The Authors 2023

1. Introduction

Nebular attenuation¹ due to dust grains is ubiquitous in astronomical observations. Besides the Galactic extinction, observations of emission lines from extragalactic H II regions inevitably suffer from different degrees of intrinsic attenuation inside their host galaxies. The overall attenuation of the emission is usually described by a parameterized nebular attenuation curve, whose shape is determined by the dust composition, grain size distribution, and geometry. Understanding the attenuation curve is important for accurate measurements of fluxes of emission lines and their derived physical properties, for example, star formation rate, metallicity, nitrogen-to-oxygen ratio (N/O), ionization parameter, among others (e.g., Pagel et al. 1979; Alloin et al. 1979; Diaz et al. 1991; Thurston et al. 1996; Dopita et al. 2000; Kennicutt & Evans 2012). Its shape also provides information about the dust sizes and compositions around emission-line regions (e.g., Spitzer 1978; Draine 2003; Flaherty et al. 2007). In addition, by comparing the average nebular attenuation curve with the average stellar attenuation curve, we can gain information about the different dust geometries around young stellar populations and older stellar populations (e.g., Charlot & Fall 2000; Inoue 2005; Chevallard et al. 2013).

Various attempts have been made to measure the nebular attenuation in both low-redshift and high-redshift emission-line galaxies (e.g., Reddy & Steidel 2004; Wild et al. 2011b; Price et al. 2014; Reddy et al. 2020; Rezaee et al. 2021). Many recent works suggest that the average nebular attenuation curve found in external galaxies has a similar shape to the average Galactic extinction curve with some slight differences (Wild et al. 2011b; Reddy et al. 2020; Rezaee et al. 2021). Although considerable dispersion around the average curve depending on galaxy properties is expected (e.g., Zafar et al. 2015; Fitzpatrick & Massa 2007; Salim & Narayanan 2020), due to the limitations of the methods and samples, this potential variation has not been investigated in depth so far. Nevertheless, there is evidence that simply applying the average Galactic extinction curve to extragalactic H II regions is inaccurate as the residual line ratio between the extinction-corrected data and the predictions from photoionization models show an unexpected dependence on the Balmer decrement (see Fig. 14 in Ji & Yan 2022). Furthermore, there could be a bias in the derived attenuation due to the complicated dust geometry, which depends on the aperture size (Vale Asari et al. 2020).

Most of the current studies on the optical attenuation curve rely on Balmer lines. In H II regions, the ratios of the intrinsic fluxes of Balmer lines are roughly fixed and are relatively insensitive to variations in temperature and density (e.g., ${\mathrm{flux}}_{\mathrm{H}\alpha }/{\mathrm{flux}}_{\mathrm{H}\beta }\approx 2.{86}_{-0.11}^{+0.18}$ for an H II region with n_H ≲ 10⁶ cm⁻³ and 5000 K ≲ T ≲ 20 000 K, Osterbrock & Ferland 2006). Therefore, any significant deviation from the theoretical ratios for Balmer lines can be interpreted as the effect of dust attenuation. However, this method can only cover a relatively small range of wavelength in the optical (roughly 4000–6600 Å). Also, higher order Balmer lines (blueward of Hγ with upper levels at n > 5) are relatively weak in observations, both because they are intrinsically weak and because they are more attenuated. In addition, absorption features in the blue part of the optical spectra add further uncertainties to the measurements of weak Balmer lines. These make this method inapplicable for a large number of galaxies whose emission-line spectra have intermediate to low signal-to-noise ratios (S/Ns). Stacking galaxy spectra could partly circumvent this problem, but this would also miss a lot of information about the variation of the nebular attenuation curve in different galaxies or regions.

Finally, it is often implicitly assumed in extragalactic studies that all of the other nebular emission lines follow the same attenuation curve probed by Balmer lines. There are two major caveats associated with this assumption. First, the attenuation for emission lines not lying in the spectral range covered by Balmer lines is not well constrained. Second, this assumption requires that all emission lines come from the same spatial location and thus experience the same level of attenuation, whose validity is rarely carefully evaluated.

Compared to high-order Balmer lines, strong optical-to-near infrared (NIR) forbidden lines cover a much larger spectral range and are easier to measure. However, the intrinsic ratios of these lines strongly depend on nebular parameters such as the metallicity and ionization parameter. To determine the intrinsic ratios of attenuation-sensitive forbidden lines, one has to rely on theoretical or empirical calibrations of nebular parameters, which could bring large systematic uncertainties since different calibrations generally do not agree with each other (e.g., Kewley & Ellison 2008; Kewley et al. 2019; Ji & Yan 2022). The key question is whether it is possible to find an empirical way to constrain these nebular parameters without specifying their values. One potential and indirect solution is to compare the attenuated spectra from galaxies with similar physical properties including redshifts, stellar masses, star formation rates (SFRs), and so forth, but different degrees of attenuation. The nebular parameters of these galaxies are likely similar given the galaxy scaling relations, and thus the differences in the attenuation-sensitive strong line ratios are likely mainly driven by dust attenuation. Wild et al. (2011a), for example, used a similar pair-matching approach to empirically measure the stellar attenuation curve in external galaxies. Still, one needs to be cautious about the scatter in the adopted galaxy scaling relations and how they transfer into variations in the galaxy spectra. Also, the uncertainties in the controlled physical parameters (e.g., stellar masses, SFRs, etc.) require careful treatments.

In this work, we adopt a slightly different approach by taking the observed quantities rather than the derived ones as the controlled physical parameters. We start with the assumption that, within the same observed star-forming (SF) region, all emission lines follow the same attenuation curve scaled by the same reddening parameter [A_V or E(B − V)]. By using a 3D line-ratio space spanned by three attenuation-insensitive line ratios, [N II]λ6583/Hα, [S II]λλ6716, 6731/Hα, and [O III]λ5007/Hβ, we constrain the nebular parameters of observed SF spatial pixels (spaxels) and thus the intrinsic ratios of forbidden lines. We show that small regions in this 3D line-ratio space correspond to a set of nearly constant nebular parameters, including metallicity (or O/H) and ionization parameter. In comparison, the variation of the dust attenuation inside these small regions can still be significant, driving most of the variations in attenuation-sensitive line ratios. This allows us to construct relative attenuation curves using a large number of extragalactic SF regions with detectable strong forbidden lines and the two strongest Balmer lines, Hα and Hβ. However, our results show that a single attenuation curve is insufficient to describe all emission lines, meaning this fundamental assumption of many extragalactic emission line studies breaks down. Based on our findings, we discuss potential explanations and suggest ways to move forward for future emission-line studies.

The layout of this paper is as follows. In Sect. 2, we introduce the observational data we use. In Sect. 3, we detail the method we use to measure the nebular attenuation curve. In Sect. 4, we show our results and compare the attenuation curve we derive with those in the literature. In Sect. 5, we examine some physical models and observational effects that could potentially explain our results. In Sect. 6, we estimate the impact of the line-specific attenuation on nebular diagnostics. We discuss the robustness of our method in Sect. 7, and draw our conclusions in Sect. 8. We also present several tests on different linear regression algorithms and emission line measurements in Appendices A and B. In Appendix C, we apply our method to a data set with high spatial resolutions from a single galaxy, IC 342. For convenience, we only show the wavelength of any given forbidden line when first introduced, unless there are other emission lines with different wavelengths but produced by the same kind of ion.

2. Data

Our main sample is drawn from the 11th product launch of the Mapping Nearby Galaxies at Apache Point Observatory survey (MaNGA, Bundy et al. 2015; Yan et al. 2016a), which is equivalent to the MaNGA products within the 17th public data release of the Sloan Digital Sky Survey (SDSS DR17, Abdurro’uf 2022). MaNGA is one of the three core programs of SDSS-IV (Blanton et al. 2017). It collected integral field unit (IFU) data using the 2.5 m Sloan telescope (Gunn et al. 2006) for over 10 000 galaxies at 0.01 < z < 0.14 with a nearly flat stellar mass distribution in 10⁹ < M_*/M_⊙ < 10¹¹ (Wake et al. 2017). Its IFUs are hexagonal bundles of optical fibers with numbers ranging from 19 to 127 (Drory et al. 2015), and a three-point dithering scheme was adopted as the main observing strategy (Law et al. 2015). Meanwhile, MaNGA used another set of small fiber bundles to target standard stars for flux calibrations, achieving an absolute calibration uncertainty less than 5% for 89% of the wavelength range (Yan et al. 2016b). The spectra were obtained through the BOSS spectrographs with a spectral resolution of R ∼ 2000 and a wavelength coverage of 3622 Å < λ < 10 354 Å (Smee et al. 2013). The observed spectra were first reduced by the Data Reduction Pipeline (DRP, Law et al. 2016, 2021a) and then processed by the Data Analysis Pipeline (DAP, Westfall et al. 2019; Belfiore et al. 2019) to produce science products including measurements of emission line fluxes. Throughout this work, we used the Gaussian fitted emission line fluxes from the DAP (see Appendix B for the results based on summed fluxes).

We selected our sample using the optical diagnostic diagrams. Specifically, we used the 3D diagnostic diagrams introduced by Ji & Yan (2020), which combines the traditional [N II]- and [S II]-based BPT diagrams (Baldwin et al. 1981; Veilleux & Osterbrock 1987). The selection function is given by

$$\begin{array}{c}\hfill P_1<-1.57P_2^2+0.53P_2-0.48,\end{array}$$(1)

where

$$\begin{array}{c}\hfill P_1=0.63N2+0.51S2+0.59R3,\end{array}$$(2)

and

$$\begin{array}{c}\hfill P_2=-0.63N2+0.78S2.\end{array}$$(3)

Here we have used N2, S2, and R3 to denote the decimal logarithms of [N II]/Hα, [S II]/Hα, and [O III]/Hβ. Before applying this criterion, we excluded spaxels that were flagged by the DAP as problematic or have S/N in Hα, Hβ, [O III], [N II], [S II], [O II]λλ3726, 3729, or [S III]λλ9069, 9531 smaller than 3. We included [O II] and [S III] during the sample selection as they are important for determining the blue end and red end of the attenuation curve.

Our selection criterion corresponds to a roughly fixed fraction of star-formation contribution based on the theoretical models of Ji & Yan (2020)². It also allows us to include some SF spaxels in the outskirt of galaxies that would have been missed by the widely adopted demarcation lines (e.g., Kewley et al. 2001; Kauffmann et al. 2003) based on the 2D diagrams. Regardless, we note that our conclusions remain largely unchanged if we use the traditional demarcation lines for selection instead. We also tried the empirical selection of SF spaxels based on gas kinematics (Law et al. 2021b), which again produced similar results. Our final sample consists of ∼2.4 × 10⁶ spaxels. Figure 1 shows the distribution of our sample in the 3D diagnostic diagram. For now we do not apply any further cut to the data. We check the effect of additional selection criteria on our results in Sect. 7.5.

Fig. 1. Distribution of the MaNGA sample viewed in a 3D line-ratio space and in one of its 2D projections. Upper left: 3D density map of the MaNGA sample, where SF spaxels are colored from black to white and other ionized regions are colored from purple to yellow. A photoionization model computed by Ji & Yan (2020) using CLOUDY is also shown, whose color coding indicates the metallicity of the simulated H II region. The blue cube is an exaggerated illustration of the 3D bins we used in this work. Upper right: the same photoionization model is shown, but color coded according to the ionization parameter of the simulated H II region. Bottom: a 2D projection of the 3D space that corresponds to an “edge-on” view of the data (color map) and model (blue). The two axes, P1 and P2, are linear combinations of the original three axes, which are denoted as N2, S2, and R3. For clarity, we only show the part of the model that covers the middle 98% of the data along the line of sight. It is clear from this view that the model surface cuts through the center of the SF spaxel distribution.

Apart from our main sample, we also made use of the MaNGA observation of a nearby galaxy IC 342, which is one of the ancillary targets in MaNGA. IC 342 is a nearby (3.3 Mpc, Saha et al. 2002) massive (M_* = 10^9.95 M_⊙, Zibetti et al. 2009) SF spiral galaxy. At the distance of IC 342, MaNGA has a spatial resolution of ∼32 pc, which is much higher than the typical resolution of 1–2 kpc for the MaNGA main sample. With IC 342, we checked the attenuation law on scales close to the size of individual H II regions. Although we only have one such nearby target in MaNGA, it provides important evidence on the scale-dependence of the nebular attenuation law. We show the analyses of the IC 342 data in Appendix C.

3. Method

In this section we introduce the method we used to measure the nebular attenuation. For a given emission line with intrinsic flux f_λ0, the attenuated flux f_λ is given by

$$\begin{array}{c}\hfill f_{\lambda }={10}^{-0.4A_{\lambda }}{f}_{\lambda 0},\end{array}$$(4)

where A_λ is the wavelength-dependent attenuation in magnitudes. The ratio of two emission lines can then be expressed as

$$\begin{array}{c}\hfill f_{{\lambda }_1}/f_{{\lambda }_2}={10}^{-0.4(A_{{\lambda }_1}-A_{{\lambda }_2})}{f}_{{\lambda }_{1}0}/f_{{\lambda }_{2}0}.\end{array}$$(5)

If the intrinsic flux ratio, f_λ10/f_λ20, is known, one can calculate the difference in attenuation between the two emission lines using the observed flux ratio. In SF H II regions, the intrinsic flux ratios of Balmer lines are roughly constants (Osterbrock & Ferland 2006). For example, under the case B recombination, the intrinsic ratio between the first two Balmer lines is roughly 2.86. Using the observed fluxes of Balmer lines, one can then measure the relative attenuation at their wavelengths, which is one of the widely adopted methods to determine the nebular attenuation curve.

Apart from Balmer lines, there are many strong forbidden lines from optical to NIR that can be observed in H II regions, such as [O II]λλ3726, 3729, [O III]λ5007, [N II]λ6583, [S II]λλ6716, 6731, [S III]λλ9069, 9531, etc. The intrinsic ratio between any two forbidden lines, however, is not a constant and varies among H II regions. In principle, the intrinsic ratio is set by the physical parameters of H II regions. To a first approximation, the variations in the intrinsic line ratios are driven by the variation in the metallicity and ionization parameter. Here the metallicity is defined as the overall chemical abundance in the gas phase, which is usually represented by the oxygen abundance, 12+log(O/H). The ionization parameter, U, is defined as the relative strength of the ionizing radiation, or $\frac{{\mathrm{\Phi }}_0}{n_{\mathrm{H}}c}$, where Φ₀ is the flux of the ionizing photons, n_H is the volume density of hydrogen, and c is the speed of light. One way to understand this approximation is to look at the BPT diagrams, where pairs of forbidden-to-Balmer line ratios with similar wavelengths are used. The ratios between these emission lines are insensitive to attenuation and are close to their intrinsic ratios, if one assumes that they follow the same attenuation curve. In these diagrams, it is shown that photoioniztion model grids of varying metallicity and ionization parameter can well cover the data locus of SF regions (e.g., Dopita et al. 2000; Kewley et al. 2001, 2019; D’Agostino et al. 2019)³.

As shown by Ji & Yan (2020), 2D diagnostic diagrams might suffer from projection effects, making it difficult to tell whether a model grid fits the data distributions in two or more diagrams in a consistent manner. By combining more than two line ratios to form a high-dimensional diagram, one can put a much more stringent constraint on the position of the ideal model that self-consistently predicts all the line ratios. In Fig. 1, we plotted our best-fit photoionization model grid for MaNGA SF spaxels (Ji & Yan 2020) by varying the metallicity and ionization parameter of a simulated H II region using the photoionization code CLOUDY (Ferland et al. 2017). Comparing the area covered by the photoionization model and the spatial distribution of the SF spaxels, we see that the location of any given datum in this line-ratio space is, to a good approximation, determined by a specific combination of a given metallicity and ionization parameter. Variations in other nebular parameters (e.g., hardness of the ionizing spectrum, abundance patterns for secondary elements, etc.) manifest themselves as scatters around the model surface, and can be covered by a more complicated model grid in principle, if we allow changes in these parameters.

With this approximation, we can rewrite Eq. (5) as

$$\begin{array}{cc}& log\frac{f_{{\lambda }_1}}{f_{{\lambda }_2}}=log{\left(\frac{f_{{\lambda }_1}}{f_{{\lambda }_2}}\right)}_0-0.4(A_{{\lambda }_1}-A_{{\lambda }_2})\hfill \\ \hfill & =F_{{\lambda }_1,{\lambda }_2}(\mathrm{O}/\mathrm{H},U,...)+\frac{A_{{\lambda }_1}-A_{{\lambda }_2}}{A_{\mathrm{H}\alpha }-A_{\mathrm{H}\beta }}log\frac{f_{\mathrm{H}\alpha }/f_{\mathrm{H}\beta }}{2.86}\hfill \\ \hfill & =F_{{\lambda }_1,{\lambda }_2}(\mathrm{O}/\mathrm{H},U,...)+m_{{\lambda }_1,{\lambda }_2}log\frac{f_{\mathrm{H}\alpha }/f_{\mathrm{H}\beta }}{2.86},\hfill \end{array}$$(6)

where F_λ1, λ2(O/H, U, ...) describes the intrinsic ratios between selected lines as a function of the metallicity, ionization parameter, and other nebular parameters of H II regions. Also, we have used the relation between the fluxes of the first two Balmer lines

$$\begin{array}{c}\hfill log(f_{\mathrm{H}\alpha }/f_{\mathrm{H}\beta }/2.86)=-0.4(A_{\mathrm{H}\alpha }-A_{\mathrm{H}\beta }).\end{array}$$(7)

If we define

$$\begin{array}{c}\hfill l_{\lambda }\equiv \frac{A_{\lambda }}{A_V},\end{array}$$(8)

to describe the normalized attenuation curve, we have

$$\begin{array}{c}\hfill m_{{\lambda }_1,{\lambda }_2}=\frac{l_{{\lambda }_1}-l_{{\lambda }_2}}{l_{\mathrm{H}\alpha }-l_{\mathrm{H}\beta }}.\end{array}$$(9)

Here m_λ1, λ2 describes the relative attenuation difference between λ₁ and λ₂ with respect to the attenuation difference between Hα and Hβ. The key of our method is to constrain the variation in F_λ1, λ2(O/H, U, ...), making it negligible compared to the variation in attenuation. If we achieve this, we can easily derive the relative attenuation, m_λ1, λ2, by fitting a linear model to the “reddening relation” between $log\frac{f_{{\lambda }_1}}{f_{{\lambda }_2}}$ and $log\frac{f_{\mathrm{H}\alpha }}{f_{\mathrm{H}\beta }}$.

We summarize our method as follows. First, we selected three line ratios that are usually thought to be attenuation-insensitive, which are [N II]/Hα, [S II]/Hα, and [O III]/Hβ. Second, we binned our sample in the 3D space spanned by log([N II]/Hα), log([S II]/Hα), and log([O III]/Hβ) (or N2, S2, and R3 for short) as illustrated in Fig. 1. Following our argument about the variation of line ratios in the 3D diagram, there is no degeneracy between O/H and U in this space. Therefore, there exists inverse functions to describe the metallicity and ionization parameter with these line ratios, that is, $(\mathrm{O}/\mathrm{H})\approx F_0^{-1}(N2,S2,R3)$ and $U\approx F_1^{-1}(N2,S2,R3)$. Thus, any small variation in O/H or U corresponds to small variations in line ratios. For example,

$$\begin{array}{c}\hfill \delta (\mathrm{O}/\mathrm{H})\approx c_0(\delta N2)+c_1(\delta S2)+c_2(\delta R3),\end{array}$$(10)

where c₀, c₁, and c₂ are coefficients also depending on N2, S2, and R3. We can then rewrite the intrinsic variations for attenuation sensitive line ratios as

$$\begin{array}{cc}& \delta F_{{\lambda }_1,{\lambda }_2}\approx a_0\delta (\mathrm{O}/\mathrm{H})+a_1\delta U+{\displaystyle \sum _i{a}_i\delta (\mathrm{secondary}\mathrm{parameters})}\hfill \\ \hfill & \approx c_0^{\prime }(\delta N2)+c_1^{\prime }(\delta S2)+c_2^{\prime }(\delta R3).\hfill \end{array}$$(11)

Again, a_i and $c_i^{\prime }$ are coefficients. It is now clear that if we make the size of each bin small enough, the variation in the metallicity, ionization parameter, and other nebular parameters would be small as well. In addition, since we binned the data in a line-ratio space that is insensitive to attenuation, we did not constrain the variation in attenuation in each bin. As a result, the variation of the second term in the rightmost hand side of Eq. (6) dominates over that of the first term. Finally, we performed a linear regression in each bin, treating log(f_λ1/f_λ2) as a linear function of log(f_Hα/f_Hβ/2.86). The slope of the relation then gives the relative attenuation with respect to the Balmer lines, as shown in Eq. (9).

We emphasize that our method does not have any dependence on photoionization models, as the selection is completely based on the data. Our only assumption is that the three line ratios we use are able to constrain the major parameters⁴ of SF spaxels, which fix their intrinsic line ratios.

In principle, the smaller the size of the bin, the more accurate our method becomes. However, reducing the bin size also reduces the number of data points within each bin. Therefore, uncertainties associated with small number statistics would become important for bins that are too small. In practice, we tried different bin sizes and eventually used a size of 0.0167 × 0.0167 × 0.0167 dex³ for each bin. With this binning scheme, we have ∼3000 bins with more than 180 data points within each bin. Meanwhile, the ±2σ range of log(f_Hα/f_Hβ/2.86) is typically $0.{20}_{-0.02}^{+0.04}\mathrm{dex}$ in each bin. Before fitting the reddening relations, we removed 3D bins with fewer than 180 data points. We also tried other cuts that select bins with more than 100 data points and 300 data points. Since the MaNGA data are more densely populated at high metallicities in the 3D line ratio space, the cut on the number of data points in each bin influences the average metallicity of our sample. We discuss the effect of different binning scheme and bin selection criteria in detail in Sect. 7.2.

Prior to binning data in 3D, we applied an attenuation correction to line ratios that form the 3D space, assuming a Fitzpatrick (1999) extinction curve with R_V = 3.1 (F99 hereafter), which is an initial guess as to the shape of the attenuation curve. This step should in principle have little effect on the derived attenuation. However, as we show in Sect. 4, our results seem to indicate different attenuation curves for different species of lines. If this is true, both the correction we applied and the assumption of the attenuation-insensitive lines become questionable, which we further discuss in Sects. 5 and 7. There is another effect associated with our choice of the line ratios that form the 3D space. Results obtained using forbidden lines that are both involved in the 3D space could be correlated by construction. This effect is discussed in Sects. 4 and 7.

The method we used to perform the linear regression for Eq. (6) is a maximum likelihood method taking into account uncertainties in both variables as well as the intrinsic scatter in the dependent variable. The natural logarithm of the likelihood function for a given data point, i, is given by

$$\begin{array}{cc}& ln{F}_{\mathrm{likelihood},i}=-0.5{(y_i-m{x}_i-b)}^2/(m^2{\sigma }_{x_i}^2+{\sigma }_{y_i}^2+{\sigma }_0^2)\hfill \\ \hfill & =-0.5{\chi }_i^2,\hfill \end{array}$$(12)

where m and b are the slope and intercept, and σ_xi, σ_yi, and σ₀ are the measurement uncertainties in x_i, y_i, and the intrinsic scatter in y (which is the same for all data points during the fit), respectively⁵. We did not consider any intrinsic scatter in x, as the Balmer ratios are relatively insensitive to temperature and density variations (Osterbrock & Ferland 2006). Still, we discuss the effect of including it in Sect. 7.1. We estimated σ_xi and σ_yi by propagating the measurement uncertainty of individual emission line. Given the test made by Belfiore et al. (2019) on emission line measurements in MaNGA, we multiplied the uncertainties reported by DAP by a factor of 1.25. We note that this adjustment does not have any significant impact on our results. We started by adopting an initial guess for the intrinsic scatter σ₀ = 0, and minimized the total χ² using the MINIMIZE function in SCIPY to obtain the slope and intercept with the Nelder & Mead (1965) method. We then checked whether the following condition is satisfied:

$$\begin{array}{c}\hfill {\displaystyle \sum _{i=1}^N{\chi }_i^2/(N-2)\le 1,}\end{array}$$(13)

where N is the number of data points. If not, we changed the intrinsic scatter by an increment of 0.0001 and recomputed the slope and intercept. We repeated the above procedure until Eq. (13) is satisfied. The method is essentially identical to the one adopted by Tremaine et al. (2002). We note that there are other ways to estimate the slope and intercept when both heteroscedastic uncertainties and intrinsic scatter are present (see e.g., Kelly 2007). We have verified the robustness of this method and present the tests in Appendix A.

4. Results

In this section we compute the relative nebular attenuation “seen” by emission lines from different species of ions and atoms. We investigate three categories of lines according to their production mechanisms as well as the ionization energies (IE) of their corresponding ions: 1) Balmer lines, including Hα, Hβ, Hγ, and Hδ; 2) high ionization lines (IE > 13.6 eV), including [Ne III]λλ3869, 3967, [O III]λ5007, and [S III]λλ9069, 9531; 3) low ionization lines (IE ≲ 13.6 eV), including [O II]λλ3726, 3729, [N II]λ6583, [S II]λλ6716, 6731, and [O I]λ6300. We summarized our results in Table 1. For the rest of the paper, we use m_{line1, line2} to represent the slope of the reddening relation, log(f_line1/f_line2) versus log(f_Hα/f_Hβ), and use $m_{{\text{line}}_{\text{1}},{\text{line}}_{\text{2}}}^{\prime }$ to represent the measured median value of the slope.

4.1. Attenuation seen by Balmer lines

We calculated two slopes related to Balmer lines: 1) the slope of log(Hα/Hγ/2.86×0.469) vs. log(Hα/Hβ/2.86) relation (or m_Hα, Hγ); 2) the slope of log(Hα/Hδ/2.86×0.259) vs. log(Hα/Hβ/2.86) relation (or m_Hα, Hδ). Here we have assumed the intrinsic Balmer ratios to follow the Case B approximation values at n_e ∼ 100 cm⁻³ and T ∼ 10⁴ K (Osterbrock & Ferland 2006). We divided the line ratios by their intrinsic values so that the intercepts, log(line ratio_obs./line ratio_theor.)|_{Hα/Hβ = 2.86}, ideally should be zero. We note that there is a small covariance between the independent and dependent variables due to the common term of Hα, which is taken into account by introducing an extra term of −2mCov(x, y) in the likelihood function (see Appendix A for more details). We also required that the S/N of all emission lines used to be greater than 3. We calculated these slopes in two different ways. The first way is the traditional method (hereafter TM), which uses the entire sample for the calculation. The second way is our new method, where we calculated slopes in 3D bins and used the median as the representative value. For Balmer lines, these two methods should give us identical results, which serves as a sanity check for our method.

The left panel of Fig. 2 shows a fitting example in one of the 3D bins, where a clear linear relation between log(Hα/Hβ) and log(Hα/Hγ) is present. The best-fit linear model shows a small uncertainty of 0.046 on the slope, and is in good agreement with the prediction from an F99 extinction curve with R_V = 3.1. Figure 3 shows the distributions of the slopes, intercepts, and intrinsic scatters of the two reddening relations. For an F99 extinction curve with R_V = 3.1, the expected slopes are $m_{\mathrm{H}\alpha ,\mathrm{H}\gamma }^{\mathrm{F}99}=1.395$ and $m_{\mathrm{H}\alpha ,\mathrm{H}\delta }^{\mathrm{F}99}=1.563$. When R_V = 4.05, these values become $m_{\mathrm{H}\alpha ,\mathrm{H}\gamma }^{\mathrm{F}99}=1.384$ and $m_{\mathrm{H}\alpha ,\mathrm{H}\delta }^{\mathrm{F}99}=1.542$. In comparison, the slopes obtained by TM are m_Hα, Hγ = 1.3400 ± 0.0005 and m_Hα, Hδ = 1.5683 ± 0.0009, where the uncertainties are calculated with a Markov chain Monte Carlo (MCMC) method using the EMCEE package in PYTHON (Foreman-Mackey et al. 2013). If we change the S/N threshold for Hγ and Hδ from 3 to 10, the resulting slopes become m_Hα, Hγ = 1.3293 ± 0.0005 and m_Hα, Hδ = 1.503 ± 0.001. While the median slopes given by our new method are $m_{\mathrm{H}\alpha ,\mathrm{H}\gamma }^{\prime }=1.343\pm 0.002$ and $m_{\mathrm{H}\alpha ,\mathrm{H}\delta }^{\prime }=1.545\pm 0.003$, where the uncertainties correspond to the 68% confidence intervals derived by using the biweight estimator of Beers et al. (1990). As a sanity check, we also ran an MCMC analysis and confirmed the uncertainties given by the MCMC were in good agreement with those returned by the biweight estimator. The 1σ widths (or the standard deviations) of the slope distributions are 0.08 and 0.16, respectively.

Fig. 2. Fitting examples in one of the 3D bins. The dashed red lines are the best-fit lines given by the maximum-likelihood method. The orange shaded regions indicate the 1σ uncertainties of the linear models. The dash-dotted blue lines are obtained by fixing the slopes using the values from an F99 extinction curve with RV = 3.1 during the fit. The error bars indicate the measurement uncertainties in different logarithmic line ratios. Left: reddening relation between Hα/Hβ and Hα/Hγ. Middle: reddening relation between Hα/Hβ and [S III]/[O III]. Right: reddening relation between Hα/Hβ and [N II]/[O II].

Fig. 3. Relative attenuation probed by Balmer lines. The first row shows the results obtained by fitting a linear function to the log(Hα/Hγ/2.86×0.469) vs. log(Hα/Hβ/2.86) relation within each 3D bin we defined. From left to right, we plot the distributions of the slope, intercept, and the intrinsic scatter. For each distribution, we show the 68% confidence interval of the median [CI(Mdn)] and the standard deviation (SD). The dashed-red lines indicate the median values among our 3D bins. The dashed-green lines indicate the values obtained by the traditional method where the whole sample is used for a single fit. The dashed-blue lines correspond to the values given by the F99 extinction curve. The second row shows similar information, but is obtained by fitting a linear function to the log(Hα/Hδ/2.86×0.259) versus log(Hα/Hβ/2.86) relation.

The results given by TM and our new method are close to each other. Our new method inevitably produces larger uncertainties. On the one hand, the intrinsic scatter in line ratios can have a larger influence on the derived slopes for bins with relatively small number of data points. On the other hand, there could be genuine variation in the shape of the attenuation curve depending on the nebular parameters or host galaxy properties, which are not included in the calculation of the uncertainties (see Sect. 7.5). The slopes derived from TM cannot be described by any single F99 extinction curve with a fixed R_V. Given the dependence of these slopes on the S/N cut, it is possible that the sample has intrinsic variations in R_V that depends on galaxy properties (see Salim & Narayanan 2020). Regardless, the corrections based on the F99 extinction curve with R_V = 3.1 would only introduce a bias less than 2% for Hα/Hγ and Hα/Hδ at A_V = 1 mag. A similar conclusion can be drawn for the median slopes derived with the new method. These results confirm that the attenuation probed by Balmer lines in galaxies can be approximately described by an average MW extinction curve (e.g., Wild et al. 2011b; Reddy et al. 2020; Rezaee et al. 2021).

The intercepts obtained by two methods also show good consistency. It is noteworthy that the intercepts we derived are slightly larger than 0. This could be due to the fact that the average intrinsic Balmer ratios are different from the values we assumed. Regardless, this offset does not affect the slopes we derived. We also notice that the median intrinsic scatter in y-axis [i.e., log(Hα/Hγ/2.86×0.469) or log(Hα/Hδ/2.86×0.259)] in each bin is small compared to the change in y. This ensures that the variation in y within each 3D bin is indeed mostly driven by the variation in x, or log(Hα/Hβ/2.86).

4.2. Attenuation seen by high ionization lines

Using our 3D binning method, we calculated slope of the log([S III]/[O III]) versus log(Hα/Hβ) relation as well as the slope of the log([O III]/[Ne III]) versus log(Hα/Hβ) relation. Since DAP ties the attenuation-uncorrected fluxes of the [S III] doublet, making f_[SIII]λ9531/f_[SIII]λ9069 ≈ 2.439, the derived slope distributions for these two lines are identical. Considering that the attenuation at NIR is relatively weak and is less dependent on wavelength, tying the observed fluxes of these two lines should not have a significant effect on our measured nebular attenuation. For [Ne III] doublet, DAP also sets a fixed ratio for the attenuation-uncorrected fluxes, making f_{[NeIII]λ3869}/f_{[NeIII]λ3967} ≈ 3.33. Although the wavelengths of [Ne III] lines are much shorter compared to [S III], the wavelength separation of the [Ne III] doublet is much smaller. The F99 extinction curve predicts that $(m_{[\mathrm{OIII}],[\mathrm{NeIII}]\lambda 3869}^{\mathrm{F}99}-m_{[\mathrm{OIII}],[\mathrm{NeIII}]\lambda 3967}^{\mathrm{F}99})/m_{[\mathrm{OIII}],[\mathrm{NeIII}]}^{\mathrm{F}99}\approx 9\%$. This is much smaller than the width of the slope distribution we derived below. Since the measurement uncertainties of [Ne III] are relatively large, an S/N cut of 3 already removes a large number of data points. To have enough 3D bins, we lowered the minimum number of spaxels required for each bin from 180 to 100. Even with this lowered cut, we only had 336 qualified 3D bins to derive the median reddening relation for log([O III]/[Ne III]).

In the middle panel of Fig. 2, we plotted the fitting result for the log([S III]/[O III]) versus log(Hα/Hβ) relation in one of the 3D bins. The overall scatter along the y axis is larger compared to the case of Balmer lines, and a large intrinsic scatter of 0.06 is returned by the fitting function, indicating there is still remaining variations in the intrinsic ratio of these forbidden lines. Regardless, the best-fit linear model matches the prediction from the F99 curve well in this case. Figure 4 shows the distribution of the slopes we derived. For the [S III] doublet, we show the result of [S III]λ9531. For the [Ne III] doublet, we show the result of [Ne III]λ3869. The slope distributions derived from the other [S III] and [Ne III] lines are identical to these. We found that the 68% confidence intervals for the central locations of the two slopes distributions are $m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\prime }=1.660\pm 0.005$ and $m_{[\mathrm{OIII}],[\mathrm{NeIII}]}^{\prime }=0.68\pm 0.06$. Meanwhile, the standard deviations of the two slope distributions are 0.27 and 1.04, respectively. Although we applied an S/N cut to [Ne III] lines, the resulting slope distribution is still very wide. Part of the reason could be due to our inclusion of the 3D bins having fewer data points. In addition, the measurements of [Ne III] fluxes could be influenced by the nearby absorption features in observed spectra. Despite the large uncertainties, our derived slopes still lie close to the values given by the F99 extinction curve with R_V = 3.1, which are $m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\mathrm{F}99}=1.695$ and $m_{[\mathrm{OIII}],[\mathrm{NeIII}]}^{\mathrm{F}99}=0.81$, after averaging the values for the two lines in each doublet. At A_V = 1 mag, using the F99 corrections would cause a bias of 1.3% for [S III]/[O III] and a bias of 5% for [O III]/[Ne III]. Therefore, the attenuation for both Balmer lines and high ionization lines in our sample galaxies can be approximately described by the F99 extinction curve. With the inclusion of the [S III] doublet, we were able to constrain the nebular attenuation curve redward of Hα.

Fig. 4. Relative attenuation probed by high ionization lines. The first row shows the results obtained by fitting a linear function to the log([S III]/[O III]) vs. log(Hα/Hβ/2.86) relation within each 3D bin we defined. From left to right: we plot the distributions of the slope, intercept, and the intrinsic scatter. For each distribution, we show the 68% confidence interval of the median [CI(Mdn)] and the standard deviation (SD). The dashed-red lines indicate the median values among our 3D bins. The dashed-blue lines correspond to the values given by the F99 extinction curve. The second row shows similar information, but is obtained by fitting a linear function to the log([Ne III]/[O III]) versus log(Hα/Hβ/2.86) relation.

From Fig. 4, one can see that the intrinsic scatter found in the reddening relations for high ionization lines is much larger compared to those for Balmer lines, which contributes to the wider slope distributions. For log([O III]/[Ne III]), the median intrinsic scatter is 0.118 dex. The standard deviation of log(Hα/Hβ) in each 3D bin is typically 0.05 dex. Combining this value with m_{[OIII],[NeIII]} ≈ 0.8, one expects that the attenuation-driven variation in log([O III]/[Ne III]) in each 3D bin is smaller than the intrinsic scatter. This further explains why the estimation of m_{[OIII],[NeIII]} is more uncertain.

4.3. Attenuation seen by low ionization lines

We examined the attenuation seen by low ionization line ratios in this subsection. Specifically, we used line ratios involving [O II] to constrain the attenuation curve blueward of Hδ.

The right panel of Fig. 2 shows a reddening relation between log([N II]/[O II]) and log(Hα/Hβ) in one of the 3D bins. In this case, unlike the linear models for Balmer lines and high ionization lines, the best-fit model shows a clear offset from the line predicted by the F99 curve. We plotted the distributions of the linear parameters of all 3D bins in Fig. 5. Our method yielded $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }=1.525\pm 0.002$, and the standard deviation of the slope distribution, σ_std, is 0.13. In comparison, the F99 extinction curve predicts that $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\mathrm{F}99}=1.844$. Our derived slope is smaller than the corresponding F99 value by roughly 2.4σ_std or 126σ′, where σ′ is the uncertainty of the median. This difference is large enough to drive a systematic bias of 13% in the attenuation-corrected [N II]/[O II] at A_V = 1 mag, if one uses the F99 extinction curve. Furthermore, the derived slope is significantly smaller than $m_{\mathrm{H}\alpha ,\mathrm{H}\delta }^{\prime }$, which cannot be explained by any single extinction curve. Since λ_[NII] > λ_Hα and λ_[OII] < λ_Hδ, for any given extinction curve that predicts the extinction to monotonically increase with decreasing wavelength in the optical, one expects A_Hδ − A_Hα < A_[OII] − A_[NII] and thus m_Hα, Hδ < m_[NII],[OII], contrary to the observed relation. Therefore, our result implies that the F99 extinction curve correctly describes the attenuation probed by Balmer lines and high ionization lines, but overpredicts the attenuation probed by low ionization lines. If we interpret this result as [N II] and [O II] both having a different magnitude of attenuation compared to that of Balmer lines, that is, A_V, Low ≠ A_V, Balmer, then from the definition of the reddening relation, we have

$$\begin{array}{c}\hfill A_{V,\mathrm{Low}}\approx \frac{m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }}{m_{[\mathrm{NII}],[\mathrm{OII}]}^{\mathrm{F}99}}{A}_{V,\mathrm{Balmer}}\approx 0.83A_{V,\mathrm{Balmer}}.\end{array}$$(14)

Fig. 5. Relative attenuation probed by low ionization lines. We fit a linear function to the log([N II]/[O II]) versus log(Hα/Hβ/2.86) relation within each 3D bin. From left to right, we plot the distributions of the slope, intercept, and the intrinsic scatter. For each distribution, we show the 68% confidence interval of the median [CI(Mdn)] and the standard deviation (SD). The dashed-red lines indicate the median values among our 3D bins. The dashed-blue lines correspond to the values given by the F99 extinction curve.

However, as we discuss in Sect. 7.3, simply using different A_V cannot explain our results and one should not use this derived “apparent A_V” to do corrections. Another complicating factor is that different low ionization lines do not necessarily have the same A_V either (see Sects. 4.4 and 5.2).

Another way to see this discrepancy is to inspect log([N II]/[O II]) corrected by the F99 extinction curve as a function of log(Hα/Hβ) in each 3D bin. Figure 6 shows such an example. We removed the fitted intercept for the log([N II]/[O II]) versus log(Hα/Hβ/2.86) relation in each 3D bin (which effectively removed the intrinsic value of log([N II]/[O II]) at zero attenuation), and plotted all the intercept-removed data in this figure. Clearly, if [N II]/[O II] is corrected by the F99 extinction law before fitting, the overall slope becomes negative. This implies that the F99 extinction law overcorrects [N II]/[O II] in our sample. Meanwhile, the distribution of $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ is relatively tight and the intrinsic scatter is small, indicating that the intrinsic scatter plays a minor role in affecting our measurement of the slope. In Fig. 6, we only show data with high S/N in Hβ. This is because the measurement uncertainty in log(Hα/Hβ/2.86) can make the attenuation-corrected trend to appear steeper than it should be as the measured value of log(Hα/Hβ/2.86) is used to derive the correction. This effect is negligible when the S/N of Hβ is high. However, we note that if there is intrinsic scatter along the x axis, a similar effect can occur even at high S/N. We discuss the effect associated with the potential intrinsic scatter in detail in Sect. 7.1.

Fig. 6. Reddening relation between log([N II]/[O II]) and log(Hα/Hβ/2.86). Only data with S/N > 30 in Hβ are included. We obtained the intrinsic log([N II]/[O II]) value at log(Hα/Hβ/2.86) = 0 and subtracted it from the observed log([N II]/[O II]) for each 3D bin. The black contours show the distribution of the data without any attenuation correction. Whereas the blue contours show the density distribution of the data with their log([N II]/[O II]) corrected by the F99 extinction curve. The contour levels trace the number density of the data points and are equally spaced on a log scale. The outermost contour and the innermost contour enclose 90% and 10% of the data, respectively. The dashed red line and the dashed green line trace the median trends in the 3D bins for the two data sets. The dash-dotted line is a horizontal line for reference.

We made a similar measurement using another combination of low ionization lines, [S II] and [O II], and found that $m_{[\mathrm{SII}],[\mathrm{OII}]}^{\prime }$ is also significantly smaller than the prediction of the F99 curve, as shown in Table 1. One might wonder if this discrepancy is due to some hidden bias induced by our 3D binning method. Indeed, [N II] and [S II] are among the emission lines we used to construct the 3D space. As a sanity check, we computed m_Hγ, [OII] using our 3D binning method. Both Hγ and [O II] are not involved in the construction of the 3D bins. We found that $m_{\mathrm{H}\gamma ,[\mathrm{OII}]}^{\prime }=0.161\pm 0.003$ (whose σ_std is 0.14), which is much smaller than $m_{\mathrm{H}\gamma ,[\mathrm{OII}]}^{\mathrm{F}99}=0.439$. This again indicates that Balmer lines and low ionization lines (in this case the [O II] line) are not subject to the same attenuation law.

4.4. Line ratios involving different species of ions

In previous subsections, we found evidence suggesting different attenuation laws for different species of lines. To further investigate this issue, we explored slopes for reddening relations involving combinations of Balmer lines, high ionization lines, and low ionization lines (which we call hybrid line ratios hereafter).

We first checked [S III]/Hα, which is a hybrid line ratio combining a high ionization line and a Balmer line. Judging from our previous measurements, one might expect that the derived slope for this line ratio should be consistent with the prediction of the F99 curve. However, we found that $m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }=0.600\pm 0.006$ (the σ_std is 0.28), which deviates from the F99 value, $m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\mathrm{F}99}=0.811$, by 0.75σ_std or 38σ′. The difference indicates a fractional bias of 8% for the attenuation-corrected [S III]/Hα at A_V = 1 mag, if one uses the F99 extinction curve for corrections. In fact, for almost all of the hybrid line ratios we checked in Table 1, the slopes of their reddening relations appear considerably lower than the corresponding F99 values.

Another noticeable effect is related to the additivity of the slope measurements. Ideally, one expects if line ratios A, B, and C satisfy log(A) = log(B)+log(C), then $m_A^{\prime }=m_B^{\prime }+m_C^{\prime }$. All line ratios should follow the law of additivity if they are subject to the same attenuation curve. From Table 1, we see that some slopes do show additivity. For example, $m_{[\mathrm{SIII}],\mathrm{H}\beta }^{\prime }=m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }+m_{\mathrm{H}\alpha ,\mathrm{H}\beta }^{\prime }=m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }+1$, $m_{[\mathrm{SIII}],[\mathrm{OII}]}^{\prime }\approx m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\prime }+m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$, etc. However, there are also a few cases where the additivity seems violated. For instance, $m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\prime }>m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }+m_{\mathrm{H}\alpha ,\mathrm{H}\beta }^{\prime }+m_{\mathrm{H}\beta ,[\mathrm{OIII}]}^{\prime }\approx m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }+0.89$, where the difference is roughly 0.16, equivalent to 1.5σ_std or 87σ′ (σ_std and σ′ are calculated from the distribution of m_{[SIII],[OIII]} − m_[SIII],Hα − m_Hα, Hβ). Also, $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }>m_{[\mathrm{NII}],[\mathrm{OIII}]}^{\prime }+m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }\approx m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }+0.89$, where the difference is roughly 0.15 (1.3σ_std or 80σ′). Interestingly, the additivity problem seems to be related to the slopes that are apparently lower than the corresponding F99 values.

One related question is whether the additivity problem comes from the bias during the fitting process associated with incorrectly estimated measurement uncertainties. In Fig. 7, we plot distributions of m_{[SIII],[OIII]} − m_[SIII],Hα − m_Hα, Hβ using samples selected based on different S/N cuts. It is clear from this figure that the additivity for these line ratios is asymptotically recovered at high S/N. In particular, we found ${\mathrm{m}}_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }$ increases noticeably with increasing S/N, while ${\mathrm{m}}_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\prime }$ and ${\mathrm{m}}_{\mathrm{H}\alpha ,[\mathrm{OIII}]}^{\prime }$ remains roughly the same. If the input measurement uncertainties are biased, the best-fit parameters returned by our likelihood function would also be biased (see Appendix A), which could lead to violation of the additivity if the bias is different for different lines. On the other hand, the bias becomes negligible when the S/N is very high. However, this explanation seems inapplicable to the case of $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$. Even after we increased the S/N cut to 20 for [N II], [O II], [O III], and Hβ, $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ is still greater than $m_{[\mathrm{NII}],[\mathrm{OIII}]}^{\prime }+m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ by 0.12 (2.1σ_std or roughly 60σ′). Also, the bias does not originate from the data lost from the S/N cut, as we performed tests with mock data and found this effect was negligible even when the S/N limit was set to 3. In Sect. 7.5 we show the additivity bias is actually more related to the surface brightness of the lines.

Fig. 7. Distributions of the slope differences under different S/N cuts. The black solid histogram shows data with S/N of [S III], [O III], and Hβ greater than 3; the dashed green histogram shows data with S/N of the aforementioned lines greater than 10; the dotted red histogram shows data with S/N of the aforementioned lines greater than 15. The slope difference peaks toward smaller values with the increasing S/N limit. However, the difference is actually due to the dependence of the slope on the surface brightness of the lines.

The slopes of the hybrid line ratios could also be affected by our construction of the 3D bins. Since we restricted the variation of [N II]/Hα, [S II]/Hα, and [O III]/Hβ in each bin, the derived slopes for certain combinations of lines are tied together. This tying should not introduce any bias if all line ratios follow the same attenuation curve, but would become problematic if different lines are attenuated by different curves. As an example, since log([S III]/Hα) = log([S III]/[N II]) + log([N II]/Hα), and log([N II]/Hα) varies little in each bin, we have forced $m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }\approx m_{[\mathrm{SIII}],[\mathrm{NII}]}^{\prime }$. Therefore, the measurement of m_[SIII],Hα is affected by m_[SIII],[NII]. Even if [S III] indeed shares the same attenuation law with Balmer lines, the involvement of the low ionization line [N II] could contribute to the deviation from the F99 value. In other words, if the observed [N II] and Hα no longer share the same attenuation law, binning data using log([N II]/Hα) can modify the measured slope for certain hybrid line ratios. In addition, our assumption of negligible variations of nebular parameters might become invalid in such cases. Mathematically, our construction of the likelihood function becomes not precise, as the single parameter σ₀ might no longer be sufficient to describe the intrinsic variation of line ratios.

Indeed, from Table 1, one can see that if a reddening relation directly or indirectly involves a low ionization line, the resulting median slope tends to be significantly lower than the expected value from the F99 curve. Interestingly, even the low ionization lines might not share the same attenuation law. One example is provided by [S III]/[O I], whose median slope is slightly smaller than that of [S III]/[N II]. If both [O I] and [N II] are subject to the same attenuation curve, we should have $m_{[\mathrm{SIII}],[\mathrm{OI}]}^{\prime }>m_{[\mathrm{SIII}],[\mathrm{NII}]}^{\prime }$. Since [O I] is an even lower ionization line compared to [N II], this discrepancy indicates a further flattening of the attenuation curve at lower ionization. There is, however, a special case where the hybrid line ratio, [N II]/Hγ, gives the slope close to the F99 value. Following the same argument we use in the preceding paragraphs, we can understand this case as follows. Even if [N II]/Hα is biased and is not attenuation-free, Hα/Hγ would still follow the unbiased attenuation curve since the slope measurements for Balmer lines do not rely on 3D binning. Given that log([N II]/Hγ) = log(Hα/Hγ) + log([N II]/Hα), and the fact that we have forced $m_{[\mathrm{NII}],\mathrm{H}\alpha }^{\prime }\approx 0$ in each 3D bin, it is no surprise $m_{[\mathrm{NII}],\mathrm{H}\gamma }^{\prime }\approx m_{\mathrm{H}\alpha ,\mathrm{H}\gamma }^{\prime }\approx m_{\mathrm{H}\alpha ,\mathrm{H}\gamma }^{F99}$.

Finally, another explanation could be that both high ionization lines and Balmer lines follow the F99 curve, but their magnitudes of attenuation differ systematically. Thus, when measuring the reddening of their combined line ratios, the result deviates from the F99 curve.

In summary, most of the hybrid line ratios we check show slopes deviating from the F99 values and some of them violate additivity as well. There could be biases associated with measurement uncertainties, which, however, cannot explain all hybrid line ratios. Since there is already some evidence that different species of lines might follow different attenuation curves, our constructions of the 3D bins could also bias the slope measurements for certain hybrid lines. To better understand these effects, we investigate some physical models that create different attenuation laws for different lines in the next section.

5. Physical model

In this section we present some physical models to explain how different observed emission lines can have different attenuation laws. In addition, we compare the model predictions with the slopes measured in our sample.

5.1. Bias from low-resolution observations

When dealing with extragalactic observations, one usually has very limited knowledge about the geometry of the dust distribution inside the corresponding galaxy. Still, indirect evidence concerning the dust distribution can be derived by comparing the attenuation of different sources. The attenuation probed by nebular emissions from SF regions is usually larger than that probed by starlight (e.g., Fanelli et al. 1988; Calzetti et al. 1994; Kashino et al. 2013; Price et al. 2014; Pannella et al. 2015; Li et al. 2021), which is attributed to the fact that the neutral gas and dust are more concentrated and clump around individual H II region (Charlot & Fall 2000; Wild et al. 2011b; Chevallard et al. 2013). Therefore, one usually assumes that emission from H II regions is “extincted” by the foreground dust rather than “attenuated” due to a complex mix of sources and dust, which is, however, not necessarily valid for unresolved “H II regions” (Pellegrini et al. 2020a).

Within each H II region, different emission lines originate from different optical depths due to their different ionization potentials. The Balmer lines are produced by recombination throughout the ionized shell. In comparison, the high ionization lines and low ionization lines are mainly produced in the inner part and outer part of the H II region, respectively. Since dust grains also exist inside the ionized shell, one might wonder whether the internal differential attenuation within the H II region can result in different attenuation laws for different lines. However, as shown by Bottorff et al. (1998), even for a highly ionized H II region with log U ∼ −2.0 and an Orion-type dust composition, the maximum optical depth at λ ∼ 0.1 μm is still smaller than 1. For the emission lines we are concerned about, the corresponding optical depths are much smaller and the dust attenuation effect should be negligible within most H II regions.

The above scenario is applicable to radiation-bounded H II regions with ionized shells well confined by the background molecular clouds. In observations, however, people found that Lyman continuum photons (Lyc) can leak from dense clouds and create emission line regions in the surrounding lower density regions (e.g., Reynolds 1990; McKee & Williams 1997; Ferguson et al. 1996; Zurita et al. 2000; Haffner et al. 2009; Howard et al. 2018), which can extend over a few kiloparsecs (see e.g., Zurita et al. 2002; Seon 2009; Belfiore et al. 2022). These low-density regions are usually referred to as the diffuse ionized gas (DIG). The ratios between low ionization lines and Balmer lines are generally higher in DIG, and emission lines produced in DIG are also less attenuated compared to H II regions.

Apart from the leaked Lyc, there is another important contributor to the emissions in DIG. As shown by Zhang et al. (2017), the leaked Lyc from H II regions cannot explain the high [O III]/Hβ observed in DIG for massive galaxies. Thus, one needs extra ionizing sources apart from H II regions to produce the observed [O III] fluxes. As suggested by Flores-Fajardo et al. (2011), the hot low-mass evolved stars (HOLMES) are likely to play an important role in ionizing DIG. Specifically, HOLMES create high [O III]/Hβ due to their harder spectra compared to H II regions. It is important to note that the enhancement of [O III]/Hβ in DIG depends on the stellar mass (or metallicity) of the galaxy, which has been confirmed on both kiloparsec and subkiloparsec scales (Zhang et al. 2017; Belfiore et al. 2022). For low-mass galaxies, [O III]/Hβ barely changes from H II regions to DIG. Whilst for high-mass galaxies, [O III]/Hβ has a conspicuous increase toward DIG. This is because low-metallicity H II regions (which are typically found in low-mass galaxies) have much higher [O III]/Hβ ratios compared to high-metallicity H II regions (as can be seen in the BPT diagrams), while DIG in general covers a relatively narrow range in [O III]/Hβ.

The above evidence suggests that both low and high ionization lines could experience different amount of attenuation compared to Balmer lines. Since the MaNGA data we use have a typical resolution of 1–2 kpc, it is possible that a considerable amount of diffuse emission has been blended with the emission from the nearby H II regions. Due to the contribution of less-attenuated DIG, we expect the measured slopes for line ratios involving low ionization lines and high ionization lines to be different from the true slopes. This appears to be consistent with our interpretation in the previous section. In addition, this effect is likely metallicity-dependent.

Figure 8 shows how the measured slopes for Hα/Hγ, [S III]/[O III], [N II]/[O II], and [O III]/[O II] correlate with the median gas-phase metallicity as well as the median Hα luminosity surface density (corrected using the F99 curve) within each bin. We derived the metallicities following the Bayesian method of Ji & Yan (2022). In short, we first computed the likelihood of the observed [N II]/Hα, [S II]/Hα, and [O III]/Hβ given the photoionization model of Ji & Yan (2020). Combining the likelihood with a flat prior in the logarithmic space, we then obtained the posterior distribution of the metallicity for each datum. Finally, the metallicity of each datum was calculated as the weighted mean using the posterior distribution. We see that m_Hα,Hγ has no dependency on the metallicity. While m_[NII],[OII] appears to have a sharp decline at very high metallicities, there seems to be no obvious trend at lower metallicities⁶. In comparison, the slopes for line ratios involving [O III] have a noticeable change with the metallicity. For high metallicity galaxies, [O III] tends to have a larger contribution from the DIG, and we see it appears less attenuated as reflected by decreasing m_{[SIII],[OIII]} and increasing m_[OIII],[OII]. However, it is worth noting that m_{[SIII],[OIII]} and m_[OIII],[OII] deviate from the F99 values at low metallicities, where the bias brought by the DIG is supposed to be lower compared to high metallicities. The contribution of DIG to [S III] and [O II] could be the complicating factor in this explanation, which we further explore in the next subsection. In addition, Fig. 8 also shows that the metallicity strongly correlates with the median Hα surface brightness of individual 3D bins. Since the Hα surface brightness is related the level of the DIG contamination (Oey et al. 2007; Zhang et al. 2017), the trends we see in Fig. 8 could also be caused by the overall change in the DIG contamination. We further discuss the effect of selecting the sample based on the Hα surface brightness in Sect. 7.5.

Fig. 8. Derived slopes for different line ratios as a function of the gas-phase metallicity. The y axes correspond to slopes of the log(fline1/fline2) versus log(Hα/Hβ/2.86) relations. The x axes correspond to metallicities derived by using Bayesian inference with a fiducial photoionization model. Each data point corresponds to a 3D bin with number of spaxels greater than 180. The 3D bins are color-coded according to the median Hα surface brightness (which is corrected by the Balmer decrement assuming an F99 extinction curve) of their spaxels. Dashed black lines and dotted black lines correspond to the median slopes and the values predicted by an F99 extinction curve, respectively.

Besides the potential DIG contamination, there is another effect that could be important due to the limited spatial resolution. For a single H II region or an association of spatially close H II regions, it is possible that not all of its emission is attenuated by the same amount. One can imagine a dust screen covers part but not all of the H II region. For such a “partially covered H II region”, the effective integrated Balmer decrement depends on the relative contributions from the covered and uncovered parts to the observed emission line spectrum. If the covered fraction varies in different H II regions, there would be an effective variation in the observed Balmer decrement, which can also give rise to line-specific effective attenuation or differential attenuation.

In summary, observations that are unable to resolve individual H II regions could lead to differential attenuation inconsistent with a unified curve either due to the distribution of dust or due to the contamination from the DIG, even though there is an underlying extinction curve at work. In the next subsection, we try setting up a few toy models to reproduce the attenuation we measured for different emission lines in the MaNGA data.

5.2. Attenuation model

We started by assuming a simple cloud model, which consists of two separate components that could have different optical depths and emission line fluxes. The first component corresponds to high attenuation, which we denote as C_High; the second component corresponds to low attenuation, which we denote as C_Low⁷. Figure 9 shows a cartoon illustration of the two-component model. Attenuated emission lines from both C_High and C_Low are mixed in observations. In this model, C_High is a single dusty H II region or a part of an H II region that is more obscured, and C_Low is the DIG surrounding the H II region or a part of the H II region that is less obscured. By definition, A_V, High > A_V, Low. We then defined the fractional contribution to the total intrinsic Hα flux from C_High to be η ≡ f_{Hα, 0(High)}/f_{Hα, 0(total)} = f_{Hα, 0(High)}/[f_{Hα, 0(High)} + f_{Hα, 0(Low)}] (0 ≤ η ≤ 1).

Fig. 9. Cartoon illustration of a two-component attenuation model. The model consists of a more attenuated component with high AV and a less attenuated component with low AV. The emission lines from both components are mixed in observations. While the high AV component is (part of) a dusty H II region, the low AV component could be a part of the H II region that is less dusty, the DIG around the H II region, or a combination of both.

The first case we considered is where both C_High and C_Low have the same intrinsic line ratios, which corresponds to a partially covered H II region. In addition, we fixed A_V, High and A_V, Low for simplicity. The only thing we varied was η, the fractional flux contributed by C_High. As a result, the observed attenuation or Balmer decrement would follow the change in η. The observed Hα flux is

$$\begin{array}{c}\hfill f_{\mathrm{H}\alpha ,\mathrm{obs}}=f_{\mathrm{H}\alpha ,0(\mathrm{total})}[\eta {10}^{-0.4A_{V,\mathrm{High}}{l}_{\mathrm{H}\alpha }}+(1-\eta ){10}^{-0.4A_{V,\mathrm{Low}}{l}_{\mathrm{H}\alpha }}].\end{array}$$(15)

Since f_Hβ, 0 = f_Hα, 0/2.86, we have

$$\begin{array}{c}\hfill {(f_{\mathrm{H}\alpha }/f_{\mathrm{H}\beta }/2.86)}_{\mathrm{obs}}=\frac{\eta {10}^{-0.4A_{V,\mathrm{High}}{l}_{\mathrm{H}\alpha }}+(1-\eta ){10}^{-0.4A_{V,\mathrm{Low}}{l}_{\mathrm{H}\alpha }}}{\eta {10}^{-0.4A_{V,\mathrm{High}}{l}_{\mathrm{H}\beta }}+(1-\eta ){10}^{-0.4A_{V,\mathrm{Low}}{l}_{\mathrm{H}\beta }}}.\end{array}$$(16)

We note that both l_Hα and l_Hβ are given by an underlying true attenuation curve, which we assume to be the F99 curve here. It is clear from the above equation that the effective V-band attenuation derived from the observed Balmer decrement, A_V, eff, changes from A_V, Low to A_V, High as η changes from 0 to 1. In general, the observed flux ratio of two lines with wavelengths λ₁ and λ₂ is given by

$$\begin{array}{c}\hfill {(f_{{\lambda }_1}/f_{{\lambda }_2})}_{\mathrm{obs}}={(\frac{f_{{\lambda }_1}}{f_{{\lambda }_2}})}_0\frac{\eta {10}^{-0.4A_{V,\mathrm{High}}{l}_{{\lambda }_1}}+(1-\eta ){10}^{-0.4A_{V,\mathrm{Low}}{l}_{{\lambda }_1}}}{\eta {10}^{-0.4A_{V,\mathrm{High}}{l}_{{\lambda }_2}}+(1-\eta ){10}^{-0.4A_{V,\mathrm{Low}}{l}_{{\lambda }_2}}},\end{array}$$(17)

where we have used the assumption that the intrinsic line ratio is the same for C_High and C_Low.

What we measured in observations is how log(f_λ1/f_λ2)_obs changes with log(f_Hα/f_Hβ/2.86)_obs, which is supposed to give us $\frac{l_{{\lambda }_1}-l_{{\lambda }_2}}{l_{\mathrm{H}\alpha }-l_{\mathrm{H}\beta }}$. However, if the change in A_V, eff is primarily driven by the change in the covering fraction, or η, we expect a modified attenuation curve to be observed. The top panel of Fig. 10 shows such an example. Instead of specifying A_V, High and A_V, Low, we set log(f_Hα/f_Hβ/2.86)_obs to 0.4 and 0 for C_High and C_Low, respectively. We see that the final m_[NII],[OII] and m_[OIII],[OII] are lower than the expected values given by the F99 curve and are closer to the median value we measured in MaNGA. The values of the Balmer decrement for C_High and C_Low are the only parameters that determine the final slope in this model. Surprisingly, this oversimplified model can already produce some changes in the measured attenuation curve. However, the weakness of this model is apparent. The result depends on the wavelengths of the lines involved rather than the species of ions that produce them. As a consequence, the slopes for Balmer lines such as m_Hα,Hγ and m_Hα,Hδ become lower than their true values, as shown in the top leftmost panel of Fig. 10, contrary to the observed results. On the other hand, for the high ionization line ratio m_{[SIII],[OIII]}, this model yields a slope that is higher than the true value, different from the observations.

Fig. 10. Effective attenuation laws generated by different models. Row a) shows several logarithmic line ratio versus logarithmic Balmer derement relations predicted by a partially covered cloud model. For comparison, we also plotted the relation predicted by the F99 extinction curve and the median relation we found in MaNGA. Row b) shows the relations predicted by a two-component model with different intrinsic line ratios. We made the H II component to have overall lower intrinsic line ratios (relative to Hα) compared to the DIG component. Row c) shows a similar model as Row b), but the underlying true attenuation curve is replaced with the CCM89 extinction curve. All y axes have arbitrary normalization.

To solve the discrepancy associated with m_{[SIII],[OIII]}, we added new parameters to the model by setting different intrinsic line ratios for C_High and C_Low. This model simulates the situation where two kinds of ionized regions, H II regions (C_High) and DIG (C_Low) are mixed in observations. If we denote r_λ, High and r_λ, Low as the intrinsic values of (f_λ/f_Hα) for C_High and C_Low, Eq. (17) becomes

$$\begin{array}{c}\hfill {(f_{{\lambda }_1}/f_{{\lambda }_2})}_{\mathrm{obs}}=\frac{r_{{\lambda }_1,\mathrm{High}}\eta {10}^{-0.4A_{V,\mathrm{High}}{l}_{{\lambda }_1}}+r_{{\lambda }_1,\mathrm{Low}}(1-\eta ){10}^{-0.4A_{V,\mathrm{Low}}{l}_{{\lambda }_1}}}{r_{{\lambda }_2,\mathrm{High}}\eta {10}^{-0.4A_{V,\mathrm{High}}{l}_{{\lambda }_2}}+r_{{\lambda }_2,\mathrm{Low}}(1-\eta ){10}^{-0.4A_{V,\mathrm{Low}}{l}_{{\lambda }_2}}}.\end{array}$$(18)

The intrinsic flux ratio r_{λ, High(Low)} should depend on the species of the ion as well as nebular parameters such as the metallicity, but for now we chose a set of values that could roughly reproduce the slopes we measured in MaNGA, which is listed in Table 2. For each line, we set a value for log(r_λ, High/r_λ, Low). From Table 2, we see that in order to fit the observed slopes, we need to make C_High (H II component) have overall lower intrinsic forbidden line ratios compared to C_Low (DIG component), which is qualitatively consistent with observations (e.g., Zhang et al. 2017; Belfiore et al. 2022)⁸.

There is, however, still a discrepancy for m_Hα,Hγ. One might wonder whether we can achieve better agreement with the observation by changing the underlying attenuation curve. Indeed, the bottom panels of Fig. 10 shows that if we instead assume a Cardelli et al. (1989) extinction curve with R_V = 3.1 (CCM89) as the true underlying extinction curve, m_Hα, Hγ gets slightly closer to the observed values. Adjusting the Balmer decrement values for C_High and C_Low can further improve the fit. In addition, we changed log(r_[SIII],High/r_[SIII],Low) from −0.20 to −0.30 in order to fit m_{[SIII],[OIII]} in this case. To determine the proper physical models and model parameters, we need evidence from high-resolution observations. One example using the MaNGA observation of the galaxy IC 342 is discussed in Appendix C.

The above models assume that the two components, C_High and C_Low have fixed magnitudes of attenuation so that the attenuation differences among H II regions are solely driven by the change in the weights of the two components. In reality, one expects the magnitude of the attenuation of the two (if not more) components also change. We can try going to another extreme by assuming that C_High and C_Low have fixed weights in their contributions to Hα and the attenuation differences among different MaNGA spaxels are solely driven by the change in the magnitudes of the attenuation for C_High and C_Low. If, for example, we keep A_V, Low constant but vary A_V, High, one expects A_V, eff to first increase with increasing A_V, High but then start to decrease at some point due to the reduction of the observed fluxes from C_High. This makes the Balmer decrement not a monotonic function of A_V, High. Besides the Balmer decrement, other line ratios would also reach an extreme value and then turn around. Figure 11 shows such a toy model for the attenuation of [N II]/[O II]. The turn-around point depends on the value of η. When η = 0.5, the turn around happens when the observed Balmer decrement changes by ∼0.05 dex from the value under A_V, Low; for η = 0.75, the turn around happens when the observed Balmer decrement changes by ∼0.10 dex from the value under A_V, Low. We did not observe this behavior clearly in our sample, indicating that this is not a dominant effect but could contribute to the scatter in the attenuation.

Fig. 11. Two-component attenuation model for [N II]/[O II], where only the attenuation of the more attenuated component is varied. Different sets of points correspond to models with different intrinsic flux contributions from the two components as well as different AV, Low. The arrow shows the direction of increasing AV, High in the model.

On the other hand, if we vary A_V, High and A_V, Low by the same amount, there is no change in the attenuation curve. This is because changing both A_Vs by the same amount effectively changes the column density of the common foreground dust screen.

The observations likely include both changes in the intrinsic flux ratios between C_High and C_Low and variations in the magnitudes of A_V, High and A_V, Low. Within our 3D bins, the latter effect might be present as scatters rather than dominating the trend for different reddening relations. To check whether these toy models are applicable to realistic data, we examined the observational evidence in the next subsection.

5.3. Comparison with observations

In Sect. 5.1, we showed how the slopes of different reddening relations change with the gas-phase metallicity in Fig. 8. Among the slope measurements we showed, m_{[SIII],[OIII]} and m_[OIII],[OII] exhibit the most clear dependence on the metallicity, which we suspected to be related to the DIG contamination on the observed [O III]. Previous studies show the [O III]/Hβ difference between the H II region and the DIG increases with increasing galaxy masses (Zhang et al. 2017; Belfiore et al. 2022), and thus increasing metallicity. The key question is whether the two-component model can explain the trends for m_{[SIII],[OIII]} and m_[OIII],[OII] given the observational evidence on the change of the [O III]/Hβ difference.

Intuitively, when the observed [O III]/Hβ has a larger contribution from the DIG, one might expect [O III] becomes less attenuated compared to other lines as the DIG contains less dust compared to H II regions. As a consequence, m_{[SIII],[OIII]} ≡ (A_[SIII] − A_[OIII])/(A_Hα − A_Hβ) would decrease and m_[OIII],[OII] ≡ (A_[OIII] − A_[OII])/(A_Hα − A_Hβ) would increase, assuming the intrinsic line ratios, r_λ, High/r_λ, Low, of other lines remain unchanged. Whereas as we have pointed out, this does not explain why m_{[SIII],[OIII]} and m_[OIII],[OII] deviate more from the expected F99 values at lower metallicities where the DIG has more similar intrinsic line ratios as H II regions.

Using our two-component model, we find the opposite of the observed trends for m_{[SIII],[OIII]} and m_[OIII],[OII]. When we fixed r_λ, High/r_λ, Low for other line ratios while increasing r_λ, High/r_λ, Low for [O III] (which effectively reduces the difference in [O III]/Hβ between H II regions and the DIG), m_{[SIII],[OIII]} became smaller and m_[OIII],[OII] became larger. This is because in our two-component model, the observed change in the Balmer decrement is driven by the change of the fractional contributions from the H II component and the DIG component. As the Balmer decrement becomes larger, the H II component gradually dominates and the observed line ratios approach the intrinsic line ratios of the H II component. Therefore, once the relative strength of [O III] increases in the H II component (i.e., less DIG contamination), our model actually predicts a boost in the increase of total [O III] when the observed Balmer decrement increases. As a result, if only the change in [O III]/Hβ is considered, the two-component model would fail to describe the metallicity dependence of m_{[SIII],[OIII]} and m_[OIII],[OII].

However, a complicating factor here is that [O III] is not the only line whose relative strength in H II regions and DIG changes with the metallicity. In fact, the value of [O II]/Hβ in the DIG is also likely a function of the stellar mass and metallicity (Zhang et al. 2017). When there are other line ratios changing together with [O III]/Hβ, the final slope for the reddening relation would depend on the both the intrinsic shape of the true extinction curve and the relative change in different line ratios. In principle we can tune the model predictions to follow the observed trends of m_{[SIII],[OIII]} and m_[OIII],[OII], but without further evidence from resolved observations we opt for not adding extra parameters related to different line ratios so as to avoid overfitting. Furthermore, in MaNGA SF spaxels, the metallicity is found to correlate with the Hα surface brightness, which is negatively correlated with the DIG contamination (Oey et al. 2007; Zhang et al. 2017). Therefore, the overall level of the DIG contamination could also change with the metallicity of the MaNGA SF spaxels, which is not captured by our simple toy model.

In summary, in order to explain the observed dependence of the reddening relations on the metallicity, our model requires the relative strengths of various lines to change in the DIG along with [O III]/Hβ. This extra assumption needs verification from resolved observations of H II regions with a range of host galaxy masses and metallicities and is beyond the scope of the current paper.

A more direct test of the scale-dependent attenuation model is to measure the attenuation of different emission lines using IFU data with high spatial resolutions. In Appendix C, we apply our method to the high-resolution IFU data within a single galaxy, IC 342. While the effect of multicomponent attenuation appears to vanish in this case, the uncertainties of the derived slopes are large due to the small number of spaxels. It would be interesting to determine a physical scale, around which the attenuations of different emission lines start to have good agreements. We plan to further investigate this problem using a larger set of high-resolution IFU data in future work.

Thus far, we have discussed the physical picture for the multicomponent attenuation, but we have not discussed its impact on observational studies. In the next section, we detail the biases from improper attenuation corrections.

6. Impact on nebular diagnositcs

There are many widely used strong-line metallicity calibrations and ionization parameter calibrations relying on attenuation-sensitive line ratios. For example, the calibration of N/O in H II regions relies on the attenuation corrected [N II]/[O II] (Thurston et al. 1996), the R23 metallicity estimator involves the attenuation corrected [O II]/Hβ (Pagel et al. 1979), the ionization parameter calibration usually involves the attenuation corrected [O III]/[O II] or [S III]/[S II] (Diaz et al. 1991; Kewley & Dopita 2002), and so forth.

If the attenuation correction is biased, the derived nebular parameters based on the line ratios above would be biased. Our results have shown that using a single attenuation curve for correcting emission lines from different transitions is inaccurate. However, how to quantify this impact is nontrivial.

If we ignored the underlying dust model and thought different lines were simply described by different attenuation curves, the fractional bias in nebular diagnostics from using a single attenuation curve was in general around 10%–20% depending on the calibration method and the overall attenuation. For example, our results show the F99 extinction curve with R_V = 3.1 fits the attenuation of Balmer line ratios well but does not fits the attenuation of [N II]/[O II]. From Table 1, one can see using the F99 curve would underestimate the intrinsic [N II]/[O II], creating an average systematic bias of 0.05 dex in N/O derived from [N II]/[O II] using Thurston et al. (1996)’s prescription at A_V = 1, where A_V is derived from the Balmer decrement. As another example, using the F99 curve to correct [S III]/[S II] would underestimate the ionization parameter by 0.08 dex at A_V = 1 (Diaz et al. 1991). Compared to the systematic difference between different strong-line calibrations (Kewley et al. 2019), the amount of bias estimated this way is not very significant.

However, as we have shown in Sect. 5, the potential existence of unresolved multicomponent attenuation cannot be simply described by different attenuation curves. In our two-component dust model, one expect the attenuation correction bias is small when the observed fluxes are dominated by C_High or C_Low, but becomes larger when C_High and C_Low have similar contributions to the observed fluxes. Using the same model as we constructed in Sect. 6, where the intrinsic log(Hα/Hβ/2.86) is set to 0.4 for C_High and 0 for C_Low, we can calculate how the line ratios corrected based on the apparent log(Hα/Hβ/2.86) differ from the intrinsic line ratios. For a partial coverage dust model, the bias due to the F99 correction in [N II]/[O II] is 0.06 dex when the apparent A_V = 1. For a H II + DIG model assuming the parameters in Table 2, the bias becomes 0.25 dex at A_V = 1, which would significantly underestimate the intrinsic N/O by 44%. Meanwhile, there is also a bias in the corrected Hα flux, which is roughly 0.23 dex at A_V = 1 in both cases, which would bias the estimated SFR by 0.16 dex or 0.23 dex depending on whether the less attenuated component is from DIG or H II regions. We note that the bias estimation strongly relies on the dust model and could vary in galaxies. For example, the larger the attenuation difference between C_High and C_Low or the larger the line ratio difference between C_High and C_Low, the larger the resulting bias. Regardless, estimating the bias using dust models are more physically motivated than simply using different attenuation curves.

We thereby conclude that the impact from using a single attenuation curve depends on the dust model and the physical conditions in galaxies. In our two-component model, the impact would be negligible if a partial coverage scenario is assumed, but would be significant if the differently attenuated components have different intrinsic line ratios. Again, to constrain the dust model, we need high-resolution data in future work.

To proceed further from the physical picture we proposed, it is important to verify the assumptions we made in our method. In the next section, we discuss the validity of our method in detail.

7. Discussion

In this section we examine three important assumptions we made in this work. First, we discuss the form of the likelihood function and check the effect of the intrinsic scatter. Second, we discuss the 3D binning scheme and check its impact on the results. Third, we check the sample selection effect and its implications.

7.1. Impact of the intrinsic scatter in the Balmer decrement

Our likelihood function, Eq. (12), considers measurement errors on both axes, but assumes the intrinsic scatter is only in y, the dependent variable. When there is truly no intrinsic scatter in x, our function can well reproduce the true slope, intercept, and the intrinsic scatter in y, as is shown in Appendix A.

We also performed tests where intrinsic scatter was manually added to x. In such cases, ignoring the intrinsic scatter in x during the fits results in an underestimation of the true slope. Also, it is noteworthy that not a great amount of the intrinsic scatter is needed to result in considerable change in the derived slope. For example, a Gaussian intrinsic scatter with σ ≈ 0.022 in the Balmer decrement is enough to make one underestimate the m_[NII],[OII] by ∼0.3, changing it from the F99 value to the median value we measured in MaNGA.

Despite its important effect, there is much evidence against the presence of the intrinsic scatter in x, the Balmer decrement. The most likely origin of the intrinsic scatter in x is the variation in the intrinsic Hα/Hβ due to the temperature and density variation (Osterbrock & Ferland 2006). However, by our construction of the 3D bins, variations in the metallicity and ionization parameter are negligibly small within each bin, which makes the temperature variations to be small as well. In fact, according to our CLOUDY model, even when we consider a wide range of metallicity and ionization parameter, with 7.39 ≤ 12 + log(O/H) ≤ 9.19 and −4.0 ≤ log(U)≤ − 1.5, the maximum variation in Hα/Hβ is only 0.022 dex from minimum to maximum. The expected variation in each bin is much smaller and cannot explain the discrepancy between the F99 values and our measurements. Furthermore, we performed another test where we added an extra dimension using the density-sensitive log([S II]λ6716/[S II]λ6731) and constructed 4D bins, within which the density variations are well constrained. The results based on the 4D bins are nearly identical to the ones based on the 3D bins, meaning that density variations do not introduce any noticeable intrinsic scatter.

Even if we assume the intrinsic scatter in the Balmer decrement were somehow responsible for the discrepancy between the observed slopes and F99, we would have to adopt different values of intrinsic scatter for different slope measurements. For example, while making m_[NII],[OII] agree with F99 requires an intrinsic scatter of 0.022 dex for the Balmer decrement, m_[OIII],[OII] needs a significantly different intrinsic scatter of 0.034 dex. In addition, slopes of Balmer line ratios and high ionization line ratios are already in good agreement with F99 and only allow an intrinsic scatter < 0.01 dex. These results are not self-consistent, as the intrinsic scatter in x should only depend on the data distribution along the x axis, which is the same for all the slope measurements using the same set of 3D bins.

However, according to Fig. 3, our likelihood function predicts nonzero intrinsic scatters for Hα/Hγ and Hα/Hδ, which seems to contradict our assumption. We note that the origin of this intrinsic scatter might come from the imperfect estimations on the measurement errors. As shown in Fig. 12, the derived intrinsic scatters for these Balmer ratios are negatively correlated with the S/N of Balmer lines. Hα/Hδ also shows larger intrinsic scatter compared to Hα/Hγ, which could come from the overall poorer detection on Hδ. Specifically, the errors resulted from the stellar continuum subtraction could play a role, which is supposed to be anticorrelated with the S/N. It is possible that the likelihood function attempts to compensate the effects from the total errors with the intrinsic scatter. If it is true, one expects Hα/Hβ to show even smaller intrinsic scatter (< 0.01), as Hβ is has overall higher S/N than Hγ⁹.

Fig. 12. Intrinsic scatters in y for the Balmer reddening relations as functions of the S/N of the weaker Balmer lines. Each point represents a 3D bin, and its median S/N is plotted. Left: how the derived intrinsic scatter in log(Hα/Hγ) varies with the S/N of Hγ. Right: how the derived intrinsic scatter in log(Hα/Hδ) varies with the S/N of Hδ.

We used a series of tests to further check if there could be any intrinsic scatter in Hα/Hβ. First, we found that when there was significant intrinsic scatter in x that was ignored, the measured intrinsic scatter for y would be overestimated, and the measured slope would be underestimated. Second, we measured a new reddening relation between log([N II]/[O II]) and log(Hα/Hγ). The F99 value for this relation is 1.321, whereas the value based on our previous measurements, $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }/m_{\mathrm{H}\alpha ,\mathrm{H}\gamma }^{\prime }$, is 1.134. When we ignored the intrinsic scatter in Hα/Hγ, we obtained a slope of 1.08, which is indeed an underestimation. If we included the intrinsic scatter in Hα/Hγ in the likelihood function, we obtained a slope of 1.145, which is more consistent with the value based on our previous measurements. In addition, the derived intrinsic scatter for log([N II]/[O II]) is 0.033, in good agreement with our previous measurement, which is 0.032. If Hα/Hβ has significant intrinsic scatter, one expects the intrinsic scatters in both log([N II]/[O II]) and log(Hα/Hγ) were overestimated previously. As a result, using the overestimated intrinsic scatter in Hα/Hγ, one should overestimate the slope of the log([N II]/[O II]) versus log(Hα/Hγ) relation. In other words, if the intrinsic scatter in Hα/Hβ is truly responsible for making $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ deviate from the F99 value, our derived slope for the log([N II]/[O II]) versus log(Hα/Hγ) relation should actually be larger than the F99 value. From our tests, we see that the intrinsic scatter in Hα/Hβ, even if present, should be negligibly small.

Finally, we would like to comment on the general approach to linear regression problems with measurement uncertainties and intrinsic scatters. The method we adopted from Tremaine et al. (2002) cannot simultaneously estimate intrinsic scatters in both axes. Only when one knows the intrinsic scatter in one axis, can one use this method to obtain an unbiased estimation of the intrinsic scatter in the other axis. The general approach by Kelly (2007) also does not deal with the case where both axes have intrinsic scatters. However, in principle one can further generalize Kelly (2007)’s approach. Another question is whether a Gaussian component is enough to describe the intrinsic scatter, especially when it is related to measurement errors. We intend to further explore the possibility of a general solution in future work.

7.2. Impact of the binning method on the derived slopes

Most of our slope measurements are based on cubic bins in the log([N II]/Hα)-log([S II]/Hα)-log([O III]/Hβ) space with number of spaxels greater than 180. We performed a series of tests to see how the size and the number cut for the 3D bins influence our results.

On the one hand, one expects the larger the bin size, the greater the intrinsic scatter within each bin, which increases the uncertainty of the measured slope. To test this effect, we made the volume of each bin eight times as big (with a bin size of 0.034 × 0.034 × 0.034 dex³), but keep the cut on the number of spaxels. While we found $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ remains roughly unchanged, the σ_std of the slope distribution increases from 0.13 to 0.17. For reddening relations that have shallower slopes and thus are more susceptible to the intrinsic scatter, one might expect the median slopes are also affected by the bin size. For example, $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ has a significant change from 0.46 to 0.30 after we adopted larger bins. However, the change in the median metallicity could play an important role in the change as well. Since the high metallicity MaNGA spaxels are more densely populated in the line-ratio space, using larger bins while keeping the cut on the number of spaxels allows more low metallicity bins to be considered while reducing the total number of high metallicity bins. This tends to lower the median slope for [O III]/[O II] (see e.g., Fig. 8).

Our conclusions also remain valid for measurements based on smaller bins. We performed another test by making the volume of each bin 5.7 times as small (with a bin size of 0.0084 × 0.0084 × 0.0116 dex³). In this case no bin has more than 180 spaxels. By lowering the number cut to 55, we obtained 2783 bins, whose $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }=1.514\pm 0.004$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }=0.508\pm 0.004$ are again close to our results from the original binning. Since each small-size bin contains fewer spaxels, the statistical uncertainty of the linear fitting becomes larger in each bin. One can see there are different uncertainties for large bins and small bins. Regardless, we have shown that our conclusions remain valid over a wide range in bin size.

On the other hand, changing the cut on the number of spaxels also has an effect on the median and the width of the slope distribution. We expect that increasing (decreasing) the lower limit of the cut reduces (enlarges) the width of the slope distribution while biasing the median slope to the high (low) metallicity value. To check this, we kept the size of each bin unchanged and performed two additional sets of fitting with bins having more than 300 data points and 100 data points, which include 527 bins and 7287 bins, respectively. We obtained $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }=1.51$ and σ_{std;[NII],[OII]} = 0.10 for bins having more than 300 data points, and $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }=1.51$ and σ_{std;[NII],[OII]} = 0.16 for bins having more than 100 data points. Whereas $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }=0.53$ and σ_{std;[OIII],[OII]} = 0.10 for bins having more than 300 data points, and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }=0.40$ and σ_{std;[OIII],[OII]} = 0.23 for bins having more than 100 data points. These results are consistent with our expectations. Specifically, $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ is increasingly reflecting the low metallicities regions as we lowered the number cut.

In brief, both the size of the 3D bins and the spaxel number cut could affect the width and the median values of the slope distribution, with the latter effect being mainly a reflection of the metallicity dependence of certain slopes. Nevertheless, our main conclusion about the inconsistent attenuation laws between different species of lines still holds, nearly irrespective of the binning method.

7.3. Bias from an attenuation dependent line-ratio space

The next important effect associated with our binning method is the correlation between different slopes. As we have seen in Sect. 4, since we constructed bins by constraining the variation in [N II]/Hα, [S II]/Hα, and [O III]/Hβ, slopes of some line ratios are tied together. This would not be a problem if all emission lines are seeing the same attenuation curve. However, if there are different attenuation laws for different lines due to, for example, contamination from the DIG, our assumption breaks down since [N II]/Hα, [S II]/Hα, and [O III]/Hβ no longer form an attenuation-insensitive space. This could bias the slope measurements. In addition, there are problems associated with the additivity for certain line ratios, especially the hybrid line ratios, which could also be related to the attenuation dependence of the 3D space.

To understand this effect, we performed a simple test to recompute $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ with a semi-empirical data set. To generate the data set, we first extracted the observed Hα and Hβ fluxes from the spaxels that were identified as SF regions in MaNGA and set them as the attenuated fluxes with no uncertainties for the simulated data. We then treated the Hα and Hβ fluxes corrected by the F99 extinction curve as the true unattenuated fluxes. The simulated “observed” fluxes were computed by adding Gaussian errors to the attenuated fluxes using the uncertainty estimations of DAP.

Next, we set the observed [N II]/Hα, [S II]/Hα, and [O III]/Hβ ratios as the intrinsic ratios for the simulated data set. Combining these with the intrinsic fluxes of Balmer lines, we obtained the intrinsic fluxes of [N II], [S II], and [O III]. However, when computing the attenuated fluxes for these lines, we used different attenuation curves for different lines. Specifically, we made $l_{[\mathrm{NII}]}=0.827l_{[\mathrm{NII}]}^{\mathrm{F}99}$ and $l_{[\mathrm{SII}]}=0.827l_{[\mathrm{SII}]}^{\mathrm{F}99}$, and $l_{[\mathrm{OIII}]}=0.948l_{[\mathrm{OIII}]}^{\mathrm{F}99}$. Here we assumed all low ionization lines share the same attenuation curve, while all high ionization lines share a different one. This assumption is oversimplified, but it can make the expected slopes including $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ to be close to our measured values in Table 1. Measurement errors were subsequently added to these lines.

The intrinsic fluxes for [O II], on the other hand, were set by a theoretical photoionization model. With the SF-ionized model computed by Ji & Yan (2020), we estimated the metallicity and ionization parameter for every SF spaxel using the assumed true values of [N II]/Hα, [S II]/Hα, and [O III]/Hβ. Then for each spaxel, we find the model prediction for [N II]/[O II] using the estimated metallicity and ionization parameter. Combining [N II]/[O II] with the intrinsic fluxes of [N II], we obtained the intrinsic fluxes for [O II]. The attenuated [O II] was then computed by setting $l_{[\mathrm{OII}]}=0.827l_{[\mathrm{OII}]}^{\mathrm{F}99}$. After we obtained the simulated observed line ratios, we added small intrinsic scatters to these line ratios. Next, we corrected the manually attenuated [N II]/Hα, [S II]/Hα, and [O III]/Hβ using the F99 curve and created 3D bins. As the final step, we remeasured $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ based on the semi-empirical data set.

As expected, the resulting slopes deviate from the F99 values due to our modification of the attenuation law. Interestingly, they are not entirely consistent with the expected values from the modified curves either. We found $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }\approx 1.59$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }\approx 0.61$, whereas the modified curves predict $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }\approx 1.53$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }\approx 0.45$. In addition, the additivity of slopes is violated and $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }-m_{[\mathrm{NII}],[\mathrm{OIII}]}^{\prime }-m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }\approx 0.09$. We note that $m_{[\mathrm{NII}],[\mathrm{OIII}]}^{\prime }\approx 0.89\approx m_{[\mathrm{NII}],[\mathrm{OIII}]}^{\mathrm{F}99}$ even after we modified the attenuation law. In comparison, the expected modified value for $m_{[\mathrm{NII}],[\mathrm{OIII}]}^{\prime }$ is 1.07. This indicates $m_{[\mathrm{NII}],[\mathrm{OIII}]}^{\prime }$ is likely fixed by our construction of the 3D bins and the preapplied extinction correction based on the F99 curve, making it almost always follow the F99 extinction curve.

When generating the semi-empirical data, we noticed another effect associated with the measurement errors. If we made the propagated measurement errors for the logarithmic line ratios Gaussians rather than made the flux errors Gaussians, the additivity seemed to be retrieved, although the slopes were still different from the expected values. In such a case, $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }\approx 1.53$, $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }\approx 0.68$, and ${|m}_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }-m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }-m_{[\mathrm{NII}],[\mathrm{OIII}]}^{\prime }|<0.01$. This effect suggests the loss of the additivity is associated with the error distributions. If the flux errors are Gaussians and the propagated errors for logarithmic line ratios deviate from Gaussians significantly, the likelihood function seems to bias the results and violate the additivity. Once we manually set the measurement errors to be 1/2 of the current values, the additivity was satisfied to a good approximation even if we assumed Gaussian errors for line fluxes. However, as we have seen in Sect. 4.4, increasing the S/N does not improve the additivity for [N II]/[O II],[O III]/[O II], and [N II]/[O III] in the MaNGA sample. In addition, we did not find such a bias associated with the error distribution in another simulated data set in Appendix A even after we manually lowered the S/N of lines to be close to 3.

In fact, by checking the dependence of the derived slopes on various physical parameters, we find the surface brightness of Hα seems more relevant to the additivity issue compared to the S/N, which we discuss in Sect. 7.5. All the evidence implies that although our semi-empirical data set can provide hints on the effects of binning biases, it is not realistic enough to reproduce our measurements with the MaNGA main sample.

In summary, when there exists different attenuation laws for different lines, the final slopes derived from 3D bins are affected but could still be different from the expected values from the modified curves. Meanwhile, the assumption of the error distribution could lead to the loss of the additivity of slopes, but it cannot explain the results we found in our sample.

7.4. Reliability of the BPT diagrams

A fundamental assumption of our 3D binning method is that the line ratio variations in the 3D-BPT space are reflecting the intrinsic variations in the metallicity, ionization parameter, and other physical parameters (except the attenuation) in SF regions. This assumption is widely adopted since the theoretical photoionization models with varying nebular parameters can well reproduce the shapes of the SF loci in the BPT diagrams (e.g., Kewley et al. 2001, 2019; Dopita et al. 2013; D’Agostino et al. 2019). In addition, although photoionization models show significant degeneracy between the metallicity and ionization parameter in the 2D-BPT diagrams, the degeneracy disappears in the 3D-BPT space (Vogt et al. 2014; Ji & Yan 2020). However, there are still several caveats we would like to discuss.

First, the SF loci in the BPT diagrams could be affected by the DIG contamination. As shown by Sanders et al. (2020), the resolved H II regions in the CHemical Abundances Of Spirals survey (CHAOS, Berg et al. 2015; Croxall et al. 2015, 2016; Berg et al. 2020) show a systematic offset from the MaNGA SF spaxels in the [S II]-based BPT diagram. If the DIG contamination is strong enough to shift observed SF regions in the BPT diagrams, there would be important systematic uncertainties when people use these line ratios to constrain metallicities and ionization parameters. However, as noted by Mannucci et al. (2021), the aperture effect for H II regions become important for high-resolution studies. When a single H II region is not fully sampled in space, the line ratios would depend on the size of the aperture and the results cannot be directly compared with photoionization models that calculate the integrated emission line spectra for the whole H II region. After taking the aperture effect into account, Mannucci et al. (2021) argue the DIG contamination only has a secondary effect on the BPT diagnostics of MaNGA SF spaxels. Therefore, the DIG contamination is likely not strong enough to shift the whole SF locus in the line-ratio space. Whereas the DIG contamination on the observed attenuation within small regions of the line-ratio space cannot be ruled out.

Second, there could be degeneracy between the evolutionary stage and the metallicity of H II regions. By modeling the time evolution of an H II region, Pellegrini et al. (2020b) were able to cover a large area in the [N II]-based BPT diagrams with a photoionization model grid at fixed metallicity. If the observed H II regions really show a significant variation in their evolutionary stages, they would show scatter similar to the metallicity variation in BPT diagrams (see also Byler et al. 2017). Still, the level of potential degeneracy introduced by this effect is still unclear. In surveys like MaNGA, the unresolved H II regions are effectively light-weighted ionized regions. Since H II regions at different evolutionary phases are characterized by different brightnesses, it is likely that the observations are biased toward some average phases, which might alleviate this degeneracy.

Finally, this work also raises a question about whether the BPT space is truly attenuation-free at low spatial resolution. Although we started our calculation with the attenuation-free assumption, the results we obtained suggest differential attenuation that depends on the species of lines, undermining the initial assumption we made. Under different attenuation curves, variations in metallicities and ionization parameters are no longer well constrained in individual bins, which is reflected by the results we show in the previous subsection. Still the amount of bias associated with this process depends on the mechanism that results in different apparent attenuation laws, which could be much more complicated than the DIG contamination model we proposed.

To conclude, despite being powerful diagnostic tools, BPT diagrams are still not fully understood in terms of their applicability and their implication. A large spectroscopic sample of well-resolved H II regions and more sophisticated photoionization models are both needed to solve the aforementioned problems.

7.5. Dependence of the derived attenuation law on galaxy properties

As summarized by Salim & Narayanan (2020), observations have shown the shape of the attenuation curve depends on a number of galaxy properties including the metallicity, stellar mass, inclination, dust column density, etc. Thus far, we have seen the relative attenuation we derived shows dependencies on the metallicity and the S/N of the lines. Some of the dependencies might come from the true variation of the attenuation law, while others might be related to the observational bias induced by differently attenuated components. In this subsection we examine these potential dependencies by checking how the sample selections based on the Hα surface brightness (Σ_Hα), inclination (b/a), and A_V influence our results. The Hα surface brightness is anticorrelated with the contamination from the DIG (Oey et al. 2007; Zhang et al. 2017), and it is also partly related to the S/N of lines and the metallicity of the galaxy. Therefore, the observational bias associated with the unresolved DIG could be reflected by the dependence of the attenuation on Σ_Hα. The inclination and A_V are closely related to the column density of the dust and could reflect the intrinsic variations of the attenuation curve due to dust geometries and compositions.

Table 3 summarizes the resulting slopes for several line ratios with different sample selection criteria. We inspected two sets of line ratios: 1) [S III]/[O III], [S III]/Hα, and Hα/[O III], 2) [N II]/[O II], [N II]/[O III], [O III]/[O II]. Ideally, slopes for each set of line ratios should satisfy additivity. When doing the sample selection, we made sure the number of spaxels selected by each selection criterion is larger than a minimum number (roughly 60 ∼ 90). Meanwhile, we used the same set of 3D bins when making cut on the same physical parameter to control the metallicity effect on the resulting slopes, though the spaxels in each bin could be somewhat different. We chose not to divide the sample according to A_V, of which the reason will be discussed later in this subsection.

For the sample selection based on Σ_Hα, we made two samples with Σ_Hα < 10³⁹ erg s⁻¹ kpc⁻² and Σ_Hα > 2 × 10³⁹ erg s⁻¹ kpc⁻², where Σ_Hα is corrected for extinction using the F99 extinction curve. One can see for the low Σ_Hα sample, $m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\prime }$ is much larger than $m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\mathrm{F}99}$ and $m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\prime }>m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }+m_{\mathrm{H}\alpha ,[\mathrm{OIII}]}^{\prime }$. On the other hand, $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ is slightly closer to $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\mathrm{F}99}$ at low Σ_Hα, but has much larger uncertainty. Meanwhile, at low Σ_Hα, $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ even becomes negative. In fact, all measured slopes except m_{Hα, [OIII]} and m_[NII],[OIII] show wide distributions at low Σ_Hα, which are reflected by the uncertainties of the medians. For $m_{\mathrm{H}\alpha ,[\mathrm{OIII}]}^{\prime }$ and $m_{[\mathrm{NII}],[\mathrm{OIII}]}^{\prime }$, their measured values are fixed by our construction of the 3D bins (see Sect. 7.3).

At high Σ_Hα, both sets of line ratios follow the additivity better, although $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ deviates more from the F99 value. The behavior of the first set of line ratios involving [S III] and [O III] seems to support the scenario of the DIG contamination, which should be reduced in brighter regions dominated by H II regions. As Σ_Hα increases, both $m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\prime }$ and $m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }$ become much closer to the expected F99 values. While for the second set of line ratios, despite the retrieval of its additivity at high Σ_Hα, $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ remain biased at high Σ_Hα, although $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ does get closer to the corresponding F99 value. One might wonder whether the slopes of Balmer lines also change with Σ_Hα. We found both $m_{\mathrm{H}\alpha ,\mathrm{H}\gamma }^{\prime }$ and $m_{\mathrm{H}\alpha ,\mathrm{H}\delta }^{\prime }$ to slightly decrease at high Σ_Hα, but $m_{\mathrm{H}\alpha ,\mathrm{H}\delta }^{\prime }$ is still greater than $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ by roughly 0.05, meaning they still cannot be fitted by the same attenuation curve.

In Sect. 4.4, we found the S/N cut was able to improve the additivity for the first set of line ratios, but failed for the second set of line ratios. Compared to the cut on S/N, the cut on Σ_Hα provides an overall better recovery of the additivity, but the physical picture behind is still unclear. In addition, the DIG contamination scenario is not able to explain the deviations of $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ from the F99 values. We also ran another test that divided the sample according to the equivalent width of Hα, which can also be used to remove the DIG (Cid Fernandes et al. 2011). However, we obtained similar results where $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ deviate from the F99 value. Still, geometric effects caused by partially covered H II regions could create a modified attenuation curve at low resolution without the DIG (see Fig. 9 and Fig. 10). Then the question becomes how the dust geometry would be different for different lines in bright H II regions. Again high-resolution data are the key to fully understanding this problem, which we will explore in future work.

For the sample selection based on b/a, we made two samples with b/a < 0.6 and b/a > 0.7, respectively. This time the differences between the slopes derived from the two samples are smaller. The high inclination sample (b/a < 0.6) shows higher $m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\prime }$ and $m_{[\mathrm{SIII}],\mathrm{H}\alpha }^{\prime }$, corresponding to a larger R_V under the parameterization of Fitzpatrick (1999). If these differences are truly driven by the change of the attenuation curve rather than the observational bias, it means the more inclined regions with larger dust column density exhibits a grayer attenuation curve, which is qualitatively consistent with the conclusion of Salim & Narayanan (2020). However, $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ also become larger at high inclination, which requires a lower R_V instead. This suggests either the intrinsic variation in the shape of the attenuation curve due to the change of the inclination is overwhelmed by the observation bias, or the intrinsic variation cannot be completely characterized by the variation in R_V.

Finally, we checked the dependence of our results on A_V. According to Salim & Narayanan (2020), A_V is a more fundamental parameter that correlates with the variation of the attenuation curve, which has been extensively investigated in works using galaxy spectral energy distribution (SED) fitting. Still, the observational determination of the A_V dependence based on galaxy SED fitting might be subject to systematic biases from correlated uncertainties (Qin et al. 2022). We calculated A_V using the Balmer decrement assuming the F99 extinction curve. However, we cannot directly construct subsamples with A_V like we did with the previous two parameters. This is because A_V is exactly the x axis in the reddening relation we defined. Making cuts on the observed independent variable with measurement errors can easily bias the linear regression results, which is similar to the “missing data problem” (see e.g., Gelman et al. 1995; Little & Rubin 2002; Kelly 2007). Indeed, we found in many cases that within a 3D bin, the slopes given by our maximum-likelihood method using the high A_V spaxels and low A_V spaxels are both larger than when using all spaxels. Therefore, instead of making subsamples, we checked how the derived slopes depend on the median A_V within each bin. We found the slope versus A_V relations are very similar to the slope versus metallicity relations. In fact, as shown by Fig. 13, the median A_V shows a strong positive correlation with the metallicity, indicating more metal-rich clouds are also dustier on average. If the changes in the slopes are mainly driven by the change in R_V due to increasing A_V, one expects in Fig. 8, $m_{[\mathrm{SIII}],[\mathrm{OIII}]}^{\prime }$ should increase with increasing metallicity (and A_V) and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$ should decrease with increasing metallicity (and A_V) as the attenuation curve become grayer. Given the opposite trends we obtained, the metallicity effect is more likely to be responsible for the changes of slopes. Due to the tight correlation between the median A_V and metallicity, we cannot find any clear dependence of the slopes on the median A_V at fixed metallicities.

Fig. 13. Correlation between the metallicity and the V-band attenuation AV for the 3D bins that contain more than 180 spaxels. The V-band attenuation is derived from the Balmer decrement. For each 3D bin, the median values of 12+log(O/H) and AV are plotted. The vertical and horizontal errorbars show the typical (median) standard deviations of AV and 12+log(O/H), respectively, in each 3D bin.

In summary, due to the observational bias possibly related to the contamination of the DIG, the relative attenuation of forbidden lines show complex dependencies on various galaxy properties. Specifically, increasing Σ_Hα seems to improve the additivity of the slopes of the reddening relations, but the deviation from the expected F99 curve remains for $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ and $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\prime }$. Changing b/a and A_V also affects the derived slopes, but the results are likely dominated by the observational bias rather than the intrinsic variation in R_V.

8. Summary and conclusions

In this work we study the nebular attenuation probed by different optical-to-NIR emission lines in MaNGA via a novel method. This method bins observed SF spaxels in a 3D line-ratio space spanned by [N II]/Hα, [S II]/Hα, and [O III]/Hβ. In theory, within a small region in this attenuation-insensitive space, variations in nebular parameters such as the metallicity and ionization parameter are well constrained, whereas the magnitude of nebular attenuation is free to vary. Therefore, variations in the attenuation-sensitive line ratios are dominated by nebular attenuation in each small 3D bin. By fitting linear models to reddening relations between logarithmic attenuation-sensitive line ratios and logarithmic Balmer decrements, we measured the nebular attenuation curve using different species of lines. We summarize our results as follows.

First, reddening relations of high ionization line ratios, including [S III]/[O III] and [O III]/[Ne III], shows median slopes consistent with the predictions of a Fitzpatrick (1999) extinction with R_V = 3.1. The same attenuation curve can well describe the attenuation of Balmer lines as well. In contrast, line ratios involving low ionization lines appear to follow a grayer attenuation curve inconsistent with high ionization line ratios and Balmer line ratios. As a result, a single attenuation curve cannot simultaneously explain the reddening of high ionization lines, Balmer lines, and low ionization lines.

Second, we suspect the low spatial resolution of MaNGA causes this discrepancy and propose a toy model that can produce different observed attenuation laws for different species of lines. In the toy model, two distinct line-emitting components with different levels of attenuation and different intrinsic line ratios are mixed and their contributions to the integrated line fluxes vary in different observations. The DIG contamination could fit into this scenario and explain the variations of the attenuation curves as a result of the different intrinsic line ratios in the DIG and H II regions, but the model parameters are not well constrained and high-resolution data are needed for verification. We found evidence in the MaNGA ancillary observations of the nearby galaxy IC 342 that at a higher spatial resolution, all emission lines follow an attenuation curve more consistent with the Fitzpatrick (1999) one. However, the uncertainties of these measurements are too large and we need a larger volume of data to confirm these results.

Third, the impact of the multicomponent dust attenuation we discovered on nebular diagnostics has a nontrivial dependence on the dust model and the physical conditions in galaxies. If one uses a single attenuation curve determined by Balmer lines to correct observed forbidden line ratios, the corrected line ratios and nebular parameters inferred from these line ratios could be biased by 0.06–0.25 dex when the apparent A_V = 1, depending on the assumptions in the attenuation model.

Finally, we checked the assumptions we made in our measurements including the intrinsic scatter, the binning method, and the sample selection. Although the intrinsic scatter in the Balmer decrement, if present, could change our slope measurements, our tests rule out any significant amount of intrinsic scatter in Hα/Hβ in our sample. We noticed the reddening relations of some line ratios are correlated due to our construction of the 3D bins. If our 3D space is no longer attenuation-free as indicated by the different attenuation laws, there would be biases associated with the measurements. Still, the main conclusion about the differential attenuation remains qualitatively unchanged, unless the 3D BPT diagrams we used are not reliable tools to reveal the overall variations in the metallicity and ionization parameter. In addition, we found strong dependence of the measured slopes on the Hα surface brightness. At high Hα surface brightness, the reddening of line ratios related to [S III] converges to the predictions of the Fitzpatrick (1999) extinction curve, whereas the reddening of line ratios involving [O II] remain different from the predictions of this curve. This result cannot be fully explained by the DIG contamination, implying the different dust geometries inside H II regions might also play a role.

Our results challenge the widely adopted assumption that different emission lines should follow the same attenuation law on all spatial scales, which is the basis for a large number of nebular diagnostics including the BPT diagnostics and metallicity calibrations. However, since the physical picture behind this phenomenon is still unclear, currently we do not have a reliable and well-tested way to make corrections to the observed emission lines, and we will try to address this issue in future work. The key to solving this issue is to perform similar analyses on large spectroscopic data sets with high spatial resolution and wide spectral coverage. Upcoming surveys such as AMASE (Yan et al. 2020) and LVM (Kollmeier et al. 2017) will provide valuable data for further investigations. Meanwhile, the PHANGS-MUSE (Emsellem et al. 2022) data, despite covering a narrower spectral range, could provide hints on the spatial distribution of dust. In addition, PHANGS-MUSE covers a few low ionization lines such as [O I]λ6300, which could be used to compare with the attenuation measurements at MaNGA’s resolution. We will compare the two data sets in a future paper.

Acknowledgments

The authors thank Kathryn Kreckel for providing the observational data of IC 342 in the data release of SDSS-IV MaNGA. The work described in this paper was partially supported by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China [Project No: CUHK 14302522]. RY also acknowledges support by the Hong Kong Global STEM Scholar scheme, by the Hong Kong Jockey Club Charities Trust through the JC STEM Lab of Astronomical Instrumentation, and by the Direct Grant of CUHK Faculty of Science. MB acknowledges support from FONDECYT regular grant 1211000 and by the ANID BASAL project FB210003. RR thanks to Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, Proj. 311223/2020-6, 304927/2017-1 and 400352/2016-8), Fundação de amparo à pesquisa do Rio Grande do Sul (FAPERGS, Proj. 16/2551-0000251-7 and 19/1750-2), and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES, Proj. 0001). RAR acknowledges the support from Conselho Nacional de Desenvolvimento Científico e Tecnológico and Fundação de Amparo à pesquisa do Estado do Rio Grande do Sul. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the US 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 https://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

Appendix A: Linear regression with measurement errors and intrinsic scatter

Fitting a straight line to a data set that has errors on both dependent and independent variables have been a problem discussed by many authors with different solutions (e.g., Press et al. 1986; Tremaine et al. 2002; Kelly 2007; Hogg et al. 2010). Since this process is of great significance for our analyses, we examine the robustness of our adopted method through some ideal simulations in this appendix.

For our situation, the problem can be described as follows. There is a data set with two observed parameters, x and y (i.e., logarithmic line ratios and logarithmic Balmer decrement), which we believe to be linearly correlated according to the theory. However, there are measurement errors for both parameters, of which the variances are estimated to be ${\sigma }_x^2$ and ${\sigma }_y^2$, respectively. Besides known observational uncertainties, there are also intrinsic variations in y (caused by variations in some nebular parameters) that have an unknown variance ${\sigma }_0^2$. Meanwhile, the intrinsic variations in x is negligible (since the Balmer line ratios are roughly constants in H II regions). The goal is to find a linear model, y = mx + b that best represents the underlying true relation.

The simple least squares method does not work because of the extra uncertainty in x. A widely adopted solution is to construct and numerically maximize a likelihood function of the model parameters. Our discussion in this appendix will mainly focus on methods based on likelihood functions. For a good review of other methods, see e.g., Kelly (2007, 2012).

Once we have the proper likelihood function, we can numerically maximize it and find the best-fit model. However, the exact form of the likelihood function as well as the maximization method varies in the literature. In the main paper, we followed the method of Tremaine et al. (2002), whose logarithmic likelihood function is given by

$$\begin{array}{c}\hfill ln{F}_0=-0.5{(y-mx-b)}^2/(m^2{\sigma }_x^2+{\sigma }_y^2+{\sigma }_0^2)=-0.5{\chi }^2.\end{array}$$(A.1)

One immediately notices that the above equation is monotonically increasing with increasing ${\sigma }_0^2$. Thus, if it is directly maximized, the solution always diverges with ${\sigma }_0^2\to \infty$. Tremaine et al. (2002)’s maximization method is to start by setting ${\sigma }_0^2=0$. After obtaining the first set of solution, (m, b), one checks the corresponding χ²/(N − 2), where N is the number of data points. If χ²/(N − 2) > 1, increase ${\sigma }_0^2$ by a small amount and recompute (m, b). This step is repeated until χ²/(N − 2)≤1.

There are, however, other potential forms of the likelihood function. For example,

$$\begin{array}{cc}\hfill ln{F}_1& =-0.5{(y-mx-b)}^2/(m^2{\sigma }_x^2+{\sigma }_y^2+{\sigma }_0^2)\hfill \\ \hfill & -0.5ln(m^2{\sigma }_x^2+{\sigma }_y^2+{\sigma }_0^2)\hfill \end{array}$$(A.2)

and

$$\begin{array}{cc}\hfill ln{F}_2& =-0.5{(y-mx-b)}^2/(m^2{\sigma }_x^2+{\sigma }_y^2+{\sigma }_0^2)\hfill \\ \hfill & -0.5ln[(m^2{\sigma }_x^2+{\sigma }_y^2+{\sigma }_0^2)/(1+m^2)].\hfill \end{array}$$(A.3)

Equation A.2 is equivalent to the logarithm of Equation 24 in Kelly (2007)¹⁰. Besides the χ² term, F₁ includes an extra term related to the total variance, which is not a constant as it involves m and ${\sigma }_0^2$. With the variance term, one can now directly maximize the likelihood function as it is no longer a monotonic function of ${\sigma }_0^2$. Equation A.3 is similar to Equation A.2, but it evaluates the variance perpendicular to the linear model rather than along the y axis, thus introducing the 1 + m² term (this term cancels out in the χ² term but not in the variance term, see e.g., Equation 32 of Hogg et al. 2010). Similar to Equation A.2, Equation A.3 can also be directly maximized.

To compare the above three likelihood methods, we set up a test based on the observed data in MaNGA. First, we generated a subsample by choosing a random spaxel classified as an SF region and selecting spaxels that lie within 0.05 dex to this spaxel in the 3D line-ratio space. By the end of this step, we obtained ∼7000 spaxels that are close to each other in the 3D line-ratio space. Second, we randomly selected 200 spaxels from this subsample. The logarithmic Balmer decrements of these spaxels were used as the intrinsic values for the independent variable, x₀. Their measurement uncertainties provided by DAP were also used to calculate σ_x. The final simulated values for x were obtained by using random variables drawn from the normal distribution $N(x_0,{\sigma }_x^2)$¹¹. The intrinsic values for the dependent variable was given by y₀ = m_truex₀ + b_true, where m_true was set to $m_{[\mathrm{OIII}],[\mathrm{OII}]}^{\mathrm{F}99}=0.95$ and b_true was set to −0.5. We then used the measurement uncertainties of log([O III]/[O II]) to represent σ_y, while setting the intrinsic scatter in y to σ₀ = 0.015. The final simulated values for y were drawn from the random distribution $y_0+N(0,{\sigma }_y^2)+N(0,{\sigma }_0^2)$. We then used the aforementioned three likelihood functions to derive m, b, and σ₀. The second step was repeated 5,000 times and distributions of these model parameters were compared.

Figure A.1 shows the comparisons between model parameters estimated by different likelihood functions. Our default function, F₀, best reproduces the true parameters for m, b, and σ₀. Although it underestimated σ₀ by yielding σ_{0, measured} = 0 for ∼30% of the data, this does not have any noticeable impact on estimations of m and b. In comparison, F₁ (F₂) tends to underestimate (overestimate) m and overestimate (underestimate) b. In addition, both F₁ and F₂ overestimate σ₀.

Fig. A.1. Distributions of slopes, intercepts, and intrinsic scatters estimated by different likelihood functions for the same simulated data set. In each panel, the vertical dashed line indicates the true value of the corresponding parameter for the simulated data set. Clearly F0 best reproduces the true parameters of the linear relation. It is also the default likelihood function we used in the main paper.

The above maximum-likelihood methods all require a good understanding of the measurement uncertainty. To check the impact of measurement uncertainties, we performed additional tests where these uncertainties were incorrectly estimated but were still fed into the likelihood functions. Table A.1 summarizes the impact of the biased measurement uncertainties on the derived median slopes. We found that when σ_x was overestimated or σ_y was underestimated, the resulting median m had an increase, making F₀ overestimate m_true. On the other hand, when σ_y was overestimated or σ_x was underestimated, the resulting median m had a decrease, making F₀ underestimate m_true. Interestingly, when both σ_x and σ_y were overestimated by the same fractional amount, the median m remained roughly the same. For example, in one of our tests where σ_x and σ_y were both overestimated by 25%, the median m returned by F₀ only differed from m_true by ∼3%. In contrast, when both σ_x and σ_y were underestimated by 25%, the median m estimated by F₀ differed from m_true by ∼20%. One can see that the fractional changes in the median m given by F₁ and F₂ when the measurement uncertainties are incorrectly estimated are actually smaller than those given by F₀ in some cases. Still, F₀ is the only likelihood function that recovers m_true when there is no bias in the measurement uncertainty. Overall F₀ appears slightly more robust against any bias in σ_y rather σ_x, probably because σ₀ can be adjusted to partly compensate the biased σ_y.

Apart from the above maximum-likelihood methods, there is a Bayesian method based on Gaussian mixture models introduced by Kelly (2007), which outperforms other methods according to Kelly (2007). Basically, Kelly (2007) approximates the intrinsic distribution of x, denoted as ψ, as Gaussians and incorporate them into the likelihood function. The posterior distributions of m, b, and ψ can be jointly inferred from the data.

We performed another test to see if our likelihood function produced similar results to those given by Kelly (2007)’s estimator. We set up a new simulated data set similar to the one produced by Kelly (2007) by drawing the intrinsic values of x from the probability distribution function (PDF)

$$\begin{array}{c}\hfill p(\xi )\propto e^{\xi }{(1+e^{2.75\xi })}^{-1}.\end{array}$$(A.4)

This PDF can be approximated by a Gaussian distribution, N(μ = −0.493, σ² = 1.1²). The true slope and intercept are set to m_true = 0.5 and b_true = 1, which gives the intrinsic distribution of y. The intrinsic scatter σ₀ is described by another Gaussian distribution, $N(0,{\sigma }_0^2={0.75}^2)$. Measurement uncertainties of x and y are independently drawn from a scaled inverse-χ² distribution whose degrees of freedom ν = 5 and scale parameters s_x = 0.5σ for x and s_y = 0.5σ₀ for y. We drew 50 data points each trial and performed a total of 5,000 trials to derived distributions of slopes, intercepts, and intrinsic scatter using F₀ as well as Kelly (2007)’s estimator.

Similar to Kelly (2007)’s test, when applying his estimator, we assumed the intrinsic distribution for x was known and used the Gaussian distribution N(μ = −0.493, σ² = 1.1²) as an approximation in Kelly (2007)’s Equation 16 to obtain the likelihood function

$$\begin{array}{c}\hfill ln{F}_3=-0.5[{(\mathbf{z}-\zeta )}^T{V}^{-1}(\mathbf{z}-\zeta )]-0.5|V|,\end{array}$$(A.5)

where z = (y, x), ζ = (m_trueμ + b_true, μ), and the covariance matrix $V=\left(\begin{array}{c}{m}_{\mathrm{true}}^2{\sigma }^2+{\sigma }_y^2+{\sigma }_0^2{m}_{\mathrm{true}}{\sigma }^2\\ m_{\mathrm{true}}{\sigma }^2{\sigma }^2+{\sigma }_x^2\end{array}\right)$. Following Kelly (2007), we simply performed maximum-likelihood estimations in this test rather than derived posteriors. Similar to F₁ and F₂, this likelihood function can be directly maximized to estimate m and b.

Figure A.2 shows the linear regression results we obtained with F₀ and F₃. It is clear that both methods not only recover the true value for each parameter equally well, but also yield similar dispersions for the distributions of the parameters. We also tried different values for the input parameters to simulate other data sets as Kelly (2007) did, and still found the two methods gave nearly identical distributions. In addition to the simulated data, we applied Kelly (2007)’s Bayesian method to several observed reddening relations in MaNGA, and again found results almost identical to what we have obtained¹². We thereby concluded our adopted method is equally good as the method recommended by Kelly (2007) in terms of the output.

Fig. A.2. Distributions of slopes, intercepts, and intrinsic scatters estimated by our default likelihood function F0 and Kelly (2007)’s likelihood function F3 for the same simulated data set. In each panel, the vertical dashed line indicates the true value of the corresponding parameter for the simulated data set.

Finally, in the aforementioned tests the measurement errors are not correlated, which apply to most of our results in the main paper. In cases when x and y contain a common line ratio, one needs to replace the variance with a covariance matrix in the likelihood function. For example, for some of the slope measurements in the main paper, both x and y have the log(Hα) term. When measuring, for example, m_Hα/Hγ, we simply replaced the variance term with $m^2{\sigma }_x^2+{\sigma }_y^2-2mVar[log(\text{H}\alpha )]+{\sigma }_0^2$, where Var[log(Hα)] is the variance of log(Hα). However, since Hα generally has much higher S/N compared to other lines, adding Var[log(Hα)] has little impact on the overall results.

Appendix B: Results based on nonparametric emission-line measurements

This work heavily relies on emission line measurements of the observed spectra. The MaNGA data were analyzed by DAP, which produced emission line measurements by using Gaussian models as templates for emission lines and fit them simultaneously with the stellar continuum models (Westfall et al. 2019; Belfiore et al. 2019). We note that there are other pipelines and codes that use different fitting schemes and might produce subtle differences for the emission line measurements (see e.g., Appendix of Belfiore et al. 2019). In addition, a single Gaussian component might not be optimal for emission line fitting in cases where multiple kinematic components are blended in observations. Belfiore et al. (2019) carefully assessed the quality of the emission line measurements of DAP, confirming the S/N reported by DAP is an appropriate metric and the emission lines are statistically well fitted by Gaussian models for S/N < 30. When S/N > 30, line profiles start to deviate from Gaussians, but the reported Gaussian fluxes are still consistent with the non-parametric fluxes obtained by direct summing. Therefore, significant bias in flux measurements and uncertainty estimations seem unlikely in our sample.

Regardless, we performed a test based on the nonparametric emission line measurements provided by DAP. These emission lines fluxes were obtained by first subtracting the best-fit stellar continuum models from the spectra, and then subtracting a baseline local to each emission line, and finally directly integrating the residual fluxes with predefined passbands (Yan et al. 2006; Westfall et al. 2019). The corresponding errors were estimated via propagation of pixelwise errors but was not fully tested or calibrated. Thus, Westfall et al. (2019) warned the usage of these error measurements¹³. In our test, we used the nonparametric fluxes and measurement errors to replace the Gaussian values and remeasure slopes for [N II]/[O II], [S III]/[O III], and [O III]/[O II], and we summarized the results in Table B.1. Similar to what we did in the main paper, we applied an S/N cut of 3 for every emission line involved. This selection criterion left us ∼8.1 × 10⁵ spaxels, which are considerably fewer compared to the Gaussian case (∼2.4 × 10⁶). This is because the estimated uncertainties for nonparametric fluxes are systematically higher than those for Gaussian fluxes due to the significant noise contribution from pixels far away from line centers. We included all 3D bins with number of spaxels greater than 80 and ended up with 1,589 bins. Slightly adjusting this bin-selection criterion does not change the overall conclusion.

We notice $m_{[\mathrm{NII}],[\mathrm{OII}]}^{\prime }$ given by nonparametric measurements is considerably larger than that given by Gaussian measurements (still it is significantly different from the values predicted by the F99 curve). For comparison, in Table B.1 we also list the results given by the Gaussian measurements using the same sample that was selected by the S/N cuts on nonparametric measurements. These Gaussian results also show some differences from those in the main paper, which could be due to the fact that the sample here was not selected using the S/N cuts on Gaussian measurements. While the nonparametric fluxes are in good agreement with Gaussian fluxes, the nonparametric flux uncertainties are overall higher. Recalling the test on measurement uncertainties in Appendix A, we suspected the difference is mainly driven by the systematically higher non-parametric flux uncertainties. Indeed, once we replaced the nonparametric uncertainties with the Gaussian uncertainties, but still using the nonparametric fluxes, the resulting slopes were lowered. Since the Gaussian uncertainties reported by MaNGA were well tested and calibrated while the non-parametric uncertainties were not (Westfall et al. 2019; Belfiore et al. 2019), we consider the Gaussian results to be more trustworthy.