© T. Ueta & M. Otsuka 2022

1. Introduction

In a recent article, Galera-Rosillo et al. (2022, hereafter GR22) have presented the results of plasma diagnostics of planetary nebulae (PNe) in the Andromeda Galaxy (M 31). The analyses were based on optical spectra of nine bright PNe (the brightest four, plus five others as control cases) taken with the OSIRIS instrument at the 10.4 m Gran Telescopio Canarias (GTC) supplemented by archival spectra. Their study was aimed at investigating the physico-chemical properties of these PNe and the progenitor mass of the central stars expected to be at the tip of the PN luminosity function (PNLF; Jacoby 1980).

The nine bright PNe in M 31 in question have been found to have achieved the maximum temperature as suggested by the post-AGB evolutionary tracks for stars with initial mass of ∼1.5 M_⊙, providing spectroscopic constraints for the stellar progenitors that define the PNLF cutoff for M 31 and for other similar star-forming galaxies. However, these PNe have also been found with the N/O abundance ratio that is 1.5 to 3 times larger than predicted for such stars (of ∼1.5 M_⊙ initial mass), indicating possible limitations for the existing theoretical models.

Meanwhile, our recent examination on a specific PN case has indicated that unless both extinction correction and plasma diagnostics are performed simultaneously and self-consistently, the results of the subsequent abundance analyses can be off by several tens of percent (Ueta & Otsuka 2021). For emission-line objects, the degree of extinction, c(Hβ)¹, is typically estimated by comparing the observed (attenuated) diagnostic line flux ratio (usually of a Balmer H line pair, most often the Hα-to-Hβ line flux ratio) with its theoretical (unattenuated) counterpart.

This represents a classic chicken-and-egg problem, because the theoretical H line ratio – the basis for extinction correction – is actually dependent on the electron density (n_e) and temperature (T_e) of the line emitting gas, which are the very quantities to be determined via plasma diagnostics that must be performed after extinction is properly corrected. Therefore, to achieve the maximum self-consistency, the theoretical H line ratio used in extinction correction must be updated simultaneously and iteratively with n_e and T_e obtained in plasma diagnostics until they converge to the optimum values. However, such an approach has been seldom practiced in the literature for some unknown reasons, most likely to alleviate the volume of non-linear numerical calculations when computational resources were still scarce (Ueta 2022). At any rate, it is simply incorrect to perform extinction correction just once with some assumed n_e and T_e before performing plasma diagnostics.

In the recent M 31 work by GR22, the Hα-to-Hβ line ratio was fixed to 2.86 (equivalent to assuming T_e = 10⁴ K and n_e = 10³ cm⁻³) in deriving the extinction, c(Hβ), irrespective of what the subsequent plasma diagnostics suggested for n_e and T_e. Here, we opt to demonstrate how this widely adopted initial assumption, namely, the fixed Hα-to-Hβ line ratio, would influence the results of extinction correction and plasma diagnostics, along with any subsequent inferences based on the results of these analyses.

2. Analyses

We adopted the measured (uncorrected) line fluxes of these nine PNe in M 31 as presented by GR22 (from their Table A.1). Regretfully, these fluxes are given in only three significant figures and without measurement uncertainties. As for the uncertainties, GR22 stated that “Total errors were determined from the quadratic propagation of the measured statistical errors of the spectra, which in general decrease from blue toward red, plus an estimated additional 3% of the measured flux in order to account for systematic errors, which include continuum determination, flux calibration and, although less important, wavelength calibration uncertainties”. Thus, without the exact uncertainties, we chose to assume a 5% uncertainty across the spectrum.

For the rest of the analyses, we followed the same procedure as outlined by GR22, as long as they are described. The apparent deviation is adopting as many valid Balmer line ratios (with respect to Hβ) as possible in extinction correction modulated by the updated n_e and T_e from plasma diagnostics in order to guarantee simultaneous self-consistency between extinction correction and plasma diagnostics. Both the extinction correction and plasma diagnostics were performed with PYNEB (Luridiana et al. 2015), using the same atomic parameter set (PYNEB_18_01, as summarized in Tables 4 and 6 in GR22). Uncertainties are assessed statistically by computing 1500 Monte Carlo simulations while allowing input line fluxes with Gaussian uncertainties.

Also, we assume a two-component nebula with high- and low-excitation regions characterized by the cut-off ionization potential (IP) of 17 eV. Transitions above this IP cut-off are considered of high-excitation and computed n_e and T_e derived from the [O III] λ4363/λ4959 and [Ar IV] λ4711/λ4740 diagnostic line ratios, while those below (of low-excitation) are calculated with the [N II] λ5755/λ6548 and [S II] λ6716/λ6731 diagnostic line ratios. A minor exception adopted by GR22 was that they used T_e([O III]) even for low-excitation regions when Te([N II]) resulted in uncertainties greater than ∼2000 K (for M2068 and M1675). As shown below, we do not make any such exceptions because the uncertainties are below 1200 K in our analyses.

3. Extinction correction

To verify that the results of our analyses can be compared squarely with those by GR22, we first emulated only a single run of extinction correction and plasma diagnostics, as performed by GR22. This means that T_e and n_e are fixed in extinction correction at 10⁴ K and 10³ cm⁻³, respectively (equivalent to adopting the fixed theoretical Hα-to-Hβ line ratio of 2.86) under the same extinction law (Cardelli et al. 1989) and the total-to-selective extinction ratio, R_V, of 3.1. Then, T_e and n_e are iterated only in plasma diagnostics. In other words, plasma diagnostic would find optimum T_e and n_e via iteration, but only after diagnostic lines are extinction-corrected with c(Hβ) based on the pre-fixed T_e = 10⁴ K and n_e = 10³ cm⁻³ (hence, the optimum T_e and n_e are most likely inconsistent with the presumed T_e and n_e). With this emulated procedure established, the only difference is the adopted input flux uncertainties. As shown in Table 1, our c(Hβ) values (“Emulated” in Table 1) in this verification are generally consistent with those obtained by GR22 (within 94.0 ± 6.5%). Thus, we consider that our calculations do indeed reproduce the results obtained by GR22 reasonably well and, hence, all the comparisons that follow are indeed valid.

In Table 1, we also quote the c(Hβ) values according to the NASA/IPAC Infrared Science Archive Galactic Dust Reddening and Extinction (GDRE) service (based on the SDSS spectra; Schlafly & Finkbeiner 2011)² The GDRE values are the mean c(Hβ) value within a sampling circle of 5 arcmin radius toward each of the nine M 31 PNe, accounting for the extinction from all but the circumstellar component for each PN, that is, the interstellar component in both the Milky Way and M 31 and the intragalactic component between the Milky Way and M 31. Thus, the GDRE values may overestimate the M 31 interstellar component because it accounts for extinction along the line of sight even beyond the target PN.

At any rate, our fully iterated c(Hβ) values (“Full” in Table 1) turn out generally greater than the GDRE values, corroborating the presence of non-zero circumstellar extinction component that varies a lot from object to object (from 0.46 to 0.1, corresponding to 65–20% attenuation). This means that individual circumstellar c(Hβ) values cannot just be neglected or generically assumed: it must be computed for each object self-consistently. At a minimum, extinction correction with fixed n_e and T_e appears already dubious.

On the whole, direct comparisons between the GR22 values and our fully iterative results show that c(Hβ) by GR22 was underestimated by more than 50% (52.5 ± 31.2 %), with individual variations from 12.6 to 88.3 %. Because c(λ) is the base-10 power-law index varying non-linearly across the spectrum, the impact of offsets in c(λ) is hard to gauge unless actually calculated. The extinction-corrected line fluxes based on the fully iterative c(Hβ) are listed in Tables A.1 and A.2, along with the original-to-revised ratio for comparison. We see that the GR22 fluxes were underestimated (not sufficiently extinction-corrected) by up to 20% at the shortest [O II] λ3727 line and overestimated (excessively extinction-corrected) by up to 59% at the longest [Ar III] λ7751 line.

4. Plasma diagnostics: Physical conditions

The resulting optimized n_e and T_e, after a full iteration of extinction correction and plasma diagnostics, are summarized in Table A.3. Convergence is achieved in less than or equal to four iterations in all cases. The ratio column shows comparisons between the GR22 values to ours. No PN shows more than two-σ deviations after one full iteration in both T_e([O III]) and T_e([N II]). This stability of T_e seems to stem from the fact that the T_e-diagnostic is not only insensitive to n_e (which is why it works as the T_e diagnostic), but also varies only slightly in T_e in the n_e–T_e diagnostic plane for the given diagnostic line ratio.

The n_e([S II]) values agree remarkably well (100 ± 3% agreement), neglecting two very deviant cases (substantially overestimated by 177 and 235% for M2068 and M1675, respectively). For these two deviant PNe, GR22 found very large uncertainties (not explicitly defined). We find the same, with about 50% uncertainties. Meanwhile, the n_e([Ar IV]) values were overestimated by GR22 (141.2 ± 24.8%). This seems to have been caused by the c(Hβ) underestimates in mitigating He Iλ4713 line contamination in the [Ar IV] λ4711 line flux, as will be described below.

Given the overall rough consistency among the n_e and T_e values obtained with and without fully iterative calculations, especially the rather invariant T_e results, one may be tempted to forgo seeking convergence in extinction correction. In determining n_e and T_e in plasma diagnostics, one might think of it as a “fitting,” for which the optimum values are found from a free excursion in the n_e − T_e space. In reality, however, the applicable range of n_e is also rather restricted, because the n_e-diagnostic curve varies from one limiting ratio to another within a narrow range of n_e given the choice of the n_e-diagnostic line.

In the n_e − T_e space, given the measured line ratios of the adopted diagnostic lines, where the solution exists (i.e., where diagnostic curves intersect) is already set. Hence, what the “fitting” does in the n_e − T_e space is effectively find the optimally closest (n_e, T_e) point to the intersection of the diagnostic curves within the measurement uncertainties. In this sense, it is expected that plasma diagnostics with the measured line ratios would yield more or less the same n_e and T_e values, irrespective of the level of rigor applied in the preceding extinction correction. A true exploration of the n_e − T_e space can only be done by allowing the measured diagnostic line ratios vary according to c(Hβ). This relative insensitivity of n_e and T_e against plasma diagnostics probably fostered a sense of “laissez faire” in the community, to the extent that the theoretical Hα-to-Hβ line ratio of 2.86 is referred to as “canonical,” when there is no canonicity whatsoever to this value. As we demonstrate below, omitting the iteration in the extinction correction is not at all a viable option, because the resulting c(Hβ) offset would seriously impact the rest of the plasma diagnostics.

5. Plasma diagnostics: Ionic abundances

As for calculating ionic abundances with PYNEB, again we followed GR22 as much as possible, in spite of the fact that at times the procedure was not explicitly described in that work. The derived ionic abundances for all the input collisionally excited lines and He recombination lines are summarized in Tables A.4 and A.5, even though GR22 presented only a subset of lines in their Table A.2.

Theoretically speaking, ionic abundances derived from the same ionic transitions should turn out identical, given the adopted n_e and T_e. In the literature, this has not necessarily been accomplished in practice. Discrepancies among the derived ionic abundances for a particular ionic species are often attributed to local temperature (and density) fluctuations in the line-emitting gas. This attribution is actually very odd, especially when the fixed “canonical” Hα-to-Hβ line ratio of 2.86 is adopted (i.e., the uniform n_e at 10³ cm⁻³ and T_e at 10⁴ K are imposed). If we are to believe varying n_e and T_e in the target nebula along the line of sight, assuming the “canonical” uniformity in n_e and T_e in the extinction correction is equivalent to “injecting inconsistency” into the analyses in the first place. This is even more so the case for 2-D plasma diagnostics with line emission maps, where spatial variation is surely expected by default. At any rate, when ionic abundances derived from multiple lines of an ionic species vary, some sort of averaging needs to be done to define a representative value. GR22 have taken the straight average, while we take the uncertainty-weighted mean.

Our derived ionic abundances show a good self-consistency nearly across the board: the derived abundances for most of the ionic species all agree within uncertainties. In all major transitions to be used in calculating the total elemental abundances (He⁺, He²⁺, O⁺, and O²⁺), with the help of the ionization correction factors (ICFs; Delgado-Inglada et al. 2014), the derived ionic abundances are consistent with each other – except for O⁺ from M1675, for which the [O II] λ7320 line was not measured among the trio of [O II] lines at 3727, 7320, and 7330 Å. Among this O⁺ trio, the [O II] λ3727 line³ is located at the blue-end of the detector bandwidth, while the [O II] λ7330 is at the red-end. Without knowing which is the more accurate one among the two, we simply took the uncertainty-weighted mean, as others have done.

GR22 obtained O⁺ abundances that were very discrepant (see their Table A.2). The O⁺λ3772 abundance came out to be less than half of the O⁺λ7320/7330 abundances, except for the recurring anomalous case of M2068 (for which the former was 2.5 times greater). There is a simple explanation as to why the O⁺λ3772 abundance came out much smaller: the underestimated c(Hβ). Here, readers are reminded that what counts is a ratio of given line flux “relative” to the Hβ flux. The [O II] λ3772 is on the blue side of Hβ, hence, the extinction-corrected flux is underestimated (i.e., not corrected enough by the underestimated c(Hβ)). On the other hand, the [O II] λ7320/7330 are on the red side of Hβ, and hence, the extinction-corrected flux is overestimated (i.e., not reduced enough by the underestimated c(Hβ)). As a result, the underestimated c(Hβ) leads to under- or over-estimated abundances on the opposite sides of the reference Hβ wavelength (as seen in Tables A.4 and A.5). The O⁺ abundance adopted by GR22 appears to have been biased toward an overestimate, as the [O II] λ7320/7330 lines are much farther away from Hβ than the [O II] λ3772 line.

The He⁺ and He²⁺ abundances by GR22 suffered from the same issue, because the He⁺ abundance was adopted from the He Iλ5876 (hence, tended to be overestimated) and the He²⁺ abundance from the He Iλ4686 (hence, tended to be underestimated). This trending is also seen in Tables A.4 and A.5. These observations reveal that anomalies in abundance analyses would arise not from offsets in the n_e and T_e values, but rather from the possible under/over-estimate by the underestimated c(Hβ) on either side of the reference wavelength of Hβ at 4861 Å. GR22 also used the corrected He⁺ and He²⁺ fluxes to constrain T_eff in comparison with theoretical models. Hence, the impact of an c(Hβ) offset can be seen at many different places in the analyses.

As for the rest of the observed lines, we find derived abundances somewhat discrepant for Ar²⁺ (7136/7751 Å), Cl²⁺ (5518/5538 Å), and S⁺ (4069/4076 Å). The [Ar III] λ7751 line is located at the red-end of the bandwidth and appears relatively strongly affected by the atmospheric absorption (by the O₂ band around 8000 Å, verified in the OSIRIS 2-D spectrum itself). The [Cl III] lines are intrinsically weak (about ∼1 when I(Hβ)=100). As for the [S II] lines, GR22 referred the [S II] λ4076 line as the [S II] λ4071 line, hence, the line identification might have been compromised. As these ionic abundances are not directly involved in the subsequent elemental abundance calculations, we leave them as they are. However, these ionic abundances, if incorrect, will affect the corresponding elemental abundances in the end. In any event, there always seem some reasonable explanations as to why the derived abundances are discrepant.

However, there is a more subtle but involved complications in the [Ar IV] λ4711/4740 lines, which play an important role in the present analyses as the high-excitation n_e-diagnostic lines. The [Ar IV] λ4711 line is contaminated by the neighboring He Iλ4713 line, whose strength needs to be estimated by scaling the He Iλ5876 measurement (the strongest, hence, the most reliable of all the detected He lines) at the corresponding n_e and T_e. Thus, the [Ar IV] λ4711 measurement is unavoidably underestimated when the He Iλ5876 measurement is overestimated by the underestimated c(Hβ). Hence, the [Ar IV] λ4711 line flux corrected for the He Iλ4713 contamination by GR22 was most likely doubly affected by this mechanism. This did indeed artificially reduce the [Ar IV] λ4711/λ4740 ratio for GR22, forcing the resulting n_e become larger by 141.2 ± 24.8%.

By the same token, the [O III] λ4363/λ4959 T_e-diagnostic ratio obtained by GR22 was at least mildly affected by the underestimated c(Hβ): the ratio was artificially reduced, forcing the resulting T_e become smaller by 96.5 ± 2.0% (Table A.3), which still remains an underestimation. Conversely, other n_e and T_e diagnostics did not seem to suffer from the underestimated c(Hβ), because in these cases with the [N II] λ5755/λ6548 and [S II] λ6716/λ6731 lines for the low-excitation region the lines are on the same side of Hβ and not far from each other. Thus, the effects of the underestimated c(Hβ) were marginalized.

6. Plasma diagnostics: Elemental abundances

As the final step of plasma diagnostics, the total elemental abundances are computed from the derived ionic abundances in terms of A(X)=12 + log(X/H), as shown in Table A.6. First, we emulated GR22 by adopting the same ICFs (Delgado-Inglada et al. 2014) with the derived ionic abundances of O⁺ and O²⁺, N⁺, Ne²⁺, S⁺, Ar²⁺, and Cl²⁺. However, there are unused ionic species such as S²⁺, Ar³⁺, and Ar⁴⁺ in the analyses. In fact, it is rather strange that the Ar³⁺ abundance was not used by GR22, because [Ar IV] was adopted as the n_e diagnostic. Because the Ar abundance can work as an important metallicity indicator, there is no reason not to adopt all three measured Ar ionic species to reduce the reliance on and uncertainty of the ICF. Hence, we also computed the total elemental abundances by adopting all observed ionic abundances. In this case, the adopted ICFs are still based on those suggested by Delgado-Inglada et al. (2014), except for Ar, for which we assume Ar = Ar²⁺ + Ar³⁺ + Ar⁴⁺.

In the ICF formulation, to yield the total elemental abundance, the sum of the observed ionic abundances for a specific element has to be scaled to account for the unobserved ionic species. These scaling factors are empirically defined as functions of He and O ionic abundances (He⁺, He²⁺, O⁺, and O²⁺; Delgado-Inglada et al. 2014). This means that the reliability of the observed ionic abundances depends on the derived abundances of these He and O ionic species.

The O elemental abundance is based on the observed O⁺ and O²⁺ ionic abundances. Both of the O⁺ and O²⁺ abundances are based on three lines, [O II] lines at 3727, 7320, and 7330 Å as well as [O III] lines at 4363, 4959, and 5007 Å, respectively. As these lines are distributed on both sides of Hβ at 4861 Å, the effects of the underestimated c(Hβ) is most likely marginal, even in the results by GR22. Indeed, the median O abundance among the nine PN sample obtained by GR22 was 8.63 ± 0.09, while ours is 8.59 ± 0.08: this is fairly consistent with relatively small uncertainties. However, this does not necessarily mean that the O abundance is insensitive to the c(Hβ) discrepancy. The present PN sample is of high-excitation, and hence, their O abundance is relatively less dependent on the more uncertain low-excitation O⁺ abundance. If targets are of low-excitation, the relative importance of the more uncertain O⁺ abundance would be greater, hence, the O abundance would have been compromised by the c(Hβ) discrepancy.

On the other hand, the N elemental abundance is based solely on the observed N⁺ ionic abundance. All three [N II] lines at 5755, 6548, and 6583 Å are located on the red-side of Hβ. Hence, their line fluxes and ionic abundances were most likely overestimated by GR22 due to the underestimated c(Hβ). The median N elemental abundance obtained was 8.27 ± 0.38 by GR22, while ours is 7.99 ± 0.35: this is a factor of 1.9 in difference. While the two are statistically indifferent given the relatively large deviation, the difference in the elemental abundance amounts to nearly 90%. When both of the N and O elemental abundances are combined to assess the N/O abundance ratio, the values come out to be 0.40 ± 0.30 (by GR22) and 0.28 ± 0.16 (by us). Again, these are statistically indifferent, but the absolute difference is an overestimate at nearly 30%.

It was further argued by GR22 that the derived N/O ratio for the M 31 PN sample was 1.5–3 times greater than expected for PNe of the ∼1.5 M_⊙ initial mass based on comparisons with theoretical models (Karakas & Lugaro 2016; Miller Bertolami 2016; Ventura et al. 2018), even referring to possible limitations of the existing models. However, we can simply interpret this as another consequence of the underestimated c(Hβ), which artificially reddened the whole spectra of target PNe.

Table 2 lists the best-fit luminosity (L_*), surface temperature (T_eff), and initial mass (M_i) of the central star, obtained via a CLOUDY model fitting by GR22 (with the c(Hβ) underestimate and the fixed solar metallicity of Z = 0.02 assumed for all PNe) and by us (without the c(Hβ) underestimate, and an appropriate metallicity in the range of Z = 0.003 − 0.009 informed from abundance analyses adopted for each PN; Otsuka & Ueta, in prep.). As has been discussed above, the apparent c(Hβ) underestimate imposed false reddening in the previous analyses. Hence, the revised models naturally suggest greater L_* (by a factor of two on average). Correspondingly, the expected initial progenitor mass for the nine-PN sample comes out to be greater.

In particular, the brightest four PNe in M 31 (M1687, M2068, M2538, and M50) are now appropriately found to be the most massive among the low-mass progenitors (1.9–2.4 M_⊙ with the average of 2.2 M_⊙, as opposed to 1.3–1.6 M_⊙ with the average of 1.5 M_⊙). For such relatively higher-mass progenitors, we would indeed expect comparatively enhanced N/O abundance ratios. However, they are by no means N over-abundant as in Type I PN (Peimbert & Torres-Peimbert 1983). Therefore, the present results from fully iterative self-consistent extinction correction and plasma diagnostics are in reasonable agreement with what is predicted by the existing theoretical models. Thus, there is no need for a greater amount of N for low-mass progenitors.

7. Concluding remarks

Plasma diagnostics determine the physico-chemical conditions of line-emitting objects in terms of n_e and T_e as well as ionic and elemental abundances. Because c(Hβ) is directly influenced by n_e and T_e, n_e, and T_e must be determined self-consistently via an iterative search for convergence through both extinction correction and plasma diagnostics. In the present exercise, we have demonstrated for a sample of nine bright PNe in M 31, how inconsistently assumed n_e and T_e in extinction correction can affect the results of the subsequent plasma diagnostics.

If c(Hβ) is not iteratively updated, and given the measured line fluxes, the resulting n_e and T_e values are more or less fixed as soon as the n_e- and T_e-diagnostic lines are selected. If c(Hβ) is not determined properly, observed line fluxes above and below the reference Hβ wavelength can be over- or under-corrected, depending on their wavelengths. Consequently, diagnostic line ratios can be erroneously amplified or reduced, especially when the lines involved are taken from both sides of Hβ. This practically means that the derived ionic and elemental abundances can be unreliable, even though the resulting n_e and T_e values may still appear reasonable. Moreover, the results of model calculations would be equally compromised, especially when such wrongly extinction-corrected lines are used as constraints.

For the present case of the nine bright PN sample in M 31, we have been able to attribute the differences between the previous and present results to c(Hβ) that was underestimated more than 50% in the previous analyses. It was because the Hα-to-Hβ line ratio was assumed to be fixed at 2.86 (equivalent to assuming n_e = 10³ cm⁻³ and T_e = 10⁴ K) in extinction correction even when the final n_e and T_e were different. We have also shown that the underestimated c(Hβ) inflicted inconsistencies in the derived ionic and elemental abundances as well as the subsequent photoionization model calculations. In the end, we have established that this bright M 31 PN sample represents the high-mass end of the low-mass PN progenitor stars of less than solar metallicities. Hence, no anomalous N overabundance has been found to raise any suspicion of irregularities in the existing evolutionary models, as suggested in previous analyses.

More specifically, the N/O abundance ratio can be significantly affected in plasma diagnostics based on optical spectra, if c(Hβ) is not determined properly. For the present case, the previously suspected N overabundance was simply caused by the underestimated c(Hβ). The empirical N abundance is based almost exclusively on the [N II] lines around 6000 Å (5755 and 6548/83 Å), on the much redder side of the reference wavelength at Hβ. Therefore, the [N II] line strengths, and hence, the N abundance can be artificially inflated when c(Hβ) is underestimated. In the mean time, the O abundance, based on multiple lines of multiple ionic species scattered across the optical spectrum on either side of Hβ, is relatively insensitive to the apparent c(Hβ) underestimate. For example, the N/O ratio is often used as an indicator of the initial stellar mass based on predictions made by theoretical models (e.g., Karakas & Lugaro 2016; Miller Bertolami 2016; Ventura et al. 2018). Hence, the erroneously estimated N/O ratio can lead to a variety of wrong conclusions beyond the physico-chemical conditions of the target PNe. This issue is of course not isolated in PNe and can happen in any line-emitting objects. Therefore, caution must be exercised when quoting abundances from the literature.

To summarize, the general lessons learned are;

• The quality of the input spectra is, of course, important; in particular, line fluxes should be free from any anomalies such as the atmospheric dispersion, sky emission, absorption, and so on.

• Extinction correction and plasma diagnostics should be performed simultaneously and self-consistently through full iteration, with a careful choice of the diagnostic lines.

• It is best to incorporate as many lines as possible in the analyses, especially from either side of the reference wavelength of extinction correction (typically at Hβ).

• The adverse effects of not performing extinction correction and plasma diagnostics self-consistently do not necessarily incur in the resulting n_e and T_e, but would be more likely in the resulting ionic and elemental abundances via diagnostic line ratios compromised by the wrongly derived value for c(Hβ).

• The incorrectly determined c(Hβ) value would systematically redden or blue the input spectrum and affect the outcomes of any subsequent analyses.

• Abundances in the literature need to be treated with care, unless plasma diagnostics were performed self-consistently in conjunction with extinction correction. It may be worthwhile to re-evaluate abundances via self-consistent extinction correction and plasma diagnostics, especially where anomalies have been reported in the previous results.

Acknowledgments

This research made use of the NASA/IPAC Infrared Science Archive Galactic Dust Reddening and Extinction (GDRE) service and PyNeb, a toolset dedicated to the analysis of emission lines (Luridiana et al. 2015). T.U. was supported partially by the Japan Society for the Promotion of Science (JSPS) through its invitation fellowship program (FY2020; L20505). M.O. was supported by JSPS Grants-in-Aid for Scientific Research(C) (JP19K03914 and 22K03675). Authors thank the anonymous referee, whose inputs helped to clarify some critical points in the manuscript.

References

Appendix A: Supplemental Materials

All Tables