© ESO 2022

1. Introduction

The [O III] λ5007 planetary nebula luminosity function (PNLF; Jacoby 1980) is the number density of planetary nebulae (PNe) in a stellar system as a function of their luminosity in that specific emission line. It has been measured in more than 50 galaxies, proving to be an important standard secondary candle on the extragalactic distance ladder out to ∼20 Mpc (see Ciardullo 2013, and references therein), and it is planned to be extended out to ∼50 Mpc (Chase et al. 2021) using novel detection and analysis techniques with the VLT Multi Unit Spectroscopic Explorer (MUSE) and similar instruments (Spriggs et al. 2020; Roth et al. 2021)

The method is based on the empirical evidence that the luminosity of a statistical set of PNe in the [O III] λ5007 nebular emission line reaches a maximum value that is invariant with galaxy type and only has a mild dependence on metallicity (e.g., Dopita et al. 1992; Ciardullo 2010). This emission includes up to 80% of the entire light emitted by a PN (typically L ∼ 10⁴ L_⊙), which, for the brightest objects, means ∼600 L_⊙ emitted in just a single line (Schönberner et al. 2007).

The bright end of the PNLF in Local Group galaxies was first modeled as an exponential by Ciardullo et al. (1989) as:

$$\begin{array}{c}\hfill N(m)\propto e^{0.307m}[1-e^{3(M^{\ast }-m)}],\end{array}$$

where N is the number of PNe in any given [O III] λ5007 magnitude bin, M^* is the bright cutoff magnitude, and m is the observed [O III] λ5007 magnitude of each PN, computed as

$$\begin{array}{c}\hfill m=-2.5·log(F_{5007})-13.74,\end{array}$$

where F₅₀₀₇ is the absolute flux of the nebula in [O III] λ5007 at a distance of 10 pc (Jacoby 1989). The current estimated value of the PNLF cutoff is M^* = −4.54 ± 0.05 (Ciardullo 2013).

Internal tests within galaxy clusters, and comparison with distances from Cepheids and surface brightness fluctuations, yield distance values that consistently agree within 10% (Ciardullo 2010). Furthermore, as PNe descend from progenitors present in any galactic subsystem, the PNLF method is equally effective for both spiral and elliptical galaxies, providing a crucial link between the Population I and Population II distance scales up to 20 Mpc. However, despite the PNLF-distance method being well tested empirically, it does not rely on a conclusive theoretical foundation (Ciardullo 2013): According to standard stellar evolution, the bright cutoff (M^*) should not be constant within galaxy types, and the brightest PNe should not be found in old populations.

Theoretical modeling (Schönberner et al. 2007; Méndez et al. 2008a) predicted that PNe with central star (CS) masses of ∼0.6 M_⊙ can attain the observed M^* if accompanied by a sufficiently delayed optically thick-to-thin transition of the nebular gas. Adopting the empirical initial-to-final mass relationships (IFMRs) of Kalirai et al. (2008), such a relatively high core mass would imply progenitors with an initial mass of ∼2 M_⊙. But, according to the evolution and visibility times of a PN, such relatively high-mass stars are too scarce (or completely absent) in old stellar systems such as elliptical galaxies or haloes to produce the number of observed bright PNe.

Although evidence points to a constant cutoff, contradictory conclusions have been reached over the years in trying to explain why M^* should be constant (Ciardullo et al. 1989; Jacoby 1989; Dopita et al. 1992; Méndez et al. 1993; Marigo et al. 2004). It is expected that the cutoff is sensitive to population age, but this fading is not observed. Paradigmatic examples are the galaxies NGC 4697 (Hubble type E6) and NGC 5128 (interacting S0p), both of which contain extremely luminous PNe, with maximum absolute magnitudes of −5.6 and −5.7, respectively (Méndez et al. 2008b; Peng et al. 2004). The massive CSs needed to generate these PNe, with lifetimes < 700 Ma, are not expected to be found in these old red populations. However, Gesicki et al. (2018) showed recently that less massive stars, M_i ∼ 1.1 − 2.0 M_⊙, can indeed produce bright PNe through the accelerated-evolution post-asymptotic giant branch (AGB) star models from Miller Bertolami (2016). Their model predicts bright PNe able to explain the PNLF and its bright cutoff at smaller CS masses and therefore larger lifetimes, ∼1 − 7 Ga.

Davis et al. (2018) tested the Gesicki et al. (2018) results in the M 31 bulge, showing that lifetimes of stars with masses able to reproduce the observed PN peak luminosities are about 1.5 Ga, too short to be present: More than ten bright PNe are found, but only ∼2 are predicted. Moreover, those authors showed that, when a proper extinction correction is accounted for, there are not only bright PNe where they are not expected, but they are substantially brighter than previously believed. Therefore, the question remains open regarding the role of progenitors with M_i ∼ 1 − 2 M_⊙ in old systems such as the M 31 bulge.

Alternatives to explain the formation of bright PNe and the PNLF invariance have been proposed, such as binary evolution (e.g., blue stragglers; Ciardullo et al. 2005), misidentification with other type of objects, or mass accretion in interacting binary stars (symbiotic stars; Soker 2006). Whether these alternative evolutionary channels are able to provide the observed number of M^* PNe is, however, not clear.

M 31 is the nearest large spiral; it has a fascinating assembly history and harbors the largest collection of PNe found in any galaxy, including the Milky Way. M 31 PNe represent an optimum laboratory for studying the PNLF and its still incomplete theoretical foundations (D’Souza & Bell 2018), attracting observers’ attention in recent years (cf. Kwitter et al. 2012; Balick et al. 2013; Corradi et al. 2015; Henry et al. 2018; Fang et al. 2015, 2018; Davis et al. 2018). In a pioneering paper, Merrett et al. (2006) cataloged 2615 likely PNe, enabling the determination of the PNLF down to ∼3.5 mag from the cutoff. More recently, Bhattacharya et al. (2021) carried out an imaging survey of M 31 over a huge area of 54 sq. deg. that included two disk regions and six inner halo substructures. Bhattacharya et al. (2021) increased the number of PN candidates to 5265 (only a small fraction of them being confirmed with spectroscopy) and extended the PNLF up to m₅₀₀₇ > 26 mag, enabling the determination of individual PN luminosity functions for each substructure. Given the large metallicity range covered by the substructures, different values for the PNLF cutoff, M^*, are expected and indeed measured by Bhattacharya et al. (2021), albeit with a slope steeper than theoretically predicted in the low metallicity regime (Dopita et al. 1992).

Given the above uncertainties, it is necessary to gain further insight into the nature of PNe at the PNLF tip in different stellar systems. As Davis et al. (2018) highlight, this problem can also lead to an incorrect estimation of a galaxy’s stellar mass, star formation rate, chemical history, and the initial mass function. Even though extragalactic PNe beyond the Magellanic Clouds are spatially unresolved, precluding knowledge of their morphology and the spatial distribution of physico-chemical properties, large telescopes provide the opportunity to obtain deep integrated spectra of the brightest PNe and to determine nebular and stellar parameters with good accuracy. With this goal, we present in this paper a detailed study of a sample of the brightest PNe in the disk of the galaxy M 31. In Sect. 2 the sample and the observations are presented. Section 3 covers the physical conditions of the observed PNe and Sect. 4 the analysis of their chemical properties. Section 5 presents the derivation of the CS parameters of each object using the results from Cloudy and post-AGB evolutionary models. Finally, we present a discussion and the conclusions derived from the study in Sects. 6 and 7.

2. Observations and data reduction

The nine target PNe were selected from Merrett et al. (2006). They are listed in Table 1. Their m₅₀₀₇ values range between 20.2 and 20.9. Four of the PNe in our sample (those marked with an asterisk in Table 1) are located on the bright tip of the PNLF. The five remaining objects, which are all within 0.5 mag of the bright tip and therefore still representative of the very bright PNe that are typically observed in external galaxies, are taken in this study as a control sample together with other objects from the literature.

Spectroscopy of these PNe was secured on different nights in September 2015 in service queue mode at the 10.4 m Gran Telescopio Canarias (GTC) telescope on the island of La Palma, Spain. The Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS) spectrograph (Cepa et al. 2000, 2003) was used in its long-slit mode, along with the R1000B grism and a slit width of 0.8″, resulting in spectral coverage from 3630 to 7850 Å with a reciprocal spectral dispersion of 2.1 Å pix⁻¹ and a spectral resolution R = 6.3 Å. During the observations, the seeing ranged from 0.6″ to 0.9″, and the data were obtained under photometric conditions and grey moon. The total exposure time per target is indicated in Table 1. The GTC calibration plan provided the G191-B2B (Oke 1990) spectrophotometric standard star for each observing night. To minimize the effects of atmospheric differential refraction, the long slit was oriented along the parallactic angle for all the PNe and standard stars.

In Table 2 the basic parameters adopted for M 31 and the observed PNe are listed as follows: Col. (2), m₅₀₀₇ as determined by Merrett et al. (2006) with quoted typical uncertainty of ∼0.07 mag; Cols. (3) and (4), RA (ξ) and Dec (η) offsets relative to the center of M 31; Cols. (5) and (6), total apparent angular distance from the center, d_app, in deg and kpc, respectively; Col. (7), de-projected distance in the plane of the disk, d_disk; Col. (8), radial velocities adopted from Merrett et al. (2006) with a typical uncertainty of ∼14 km s⁻¹; Col. (9), expected radial velocity according to the kinematic model of the extended disk of Ibata et al. (2005); and Col. (10), differences between the measured (Col. 8) and expected (Col. 9) velocities.

Small velocity differences in Col. (10) of Table 2 give further confidence that the observed PNe belong to the disk of M 31, large ones (≥100 km s⁻¹ such as for PNe M2538* and M1074) may indicate that these PNe are associated with other M 31 substructures, also considering their large de-projected galactocentric distances. Positions on the sky of the target PNe are shown in Fig. 1.

Fig. 1. Location of the target PNe. North is up, and east is to the left. Green squares correspond to the brightest sample targets and magenta diamonds to the control sample. More details are given in the text.

The GTC spectroscopic data were reduced following the standard procedure using IRAF¹. Spectra were bias-subtracted, flat-fielded, and wavelength-calibrated in 2D. Cosmic rays were removed using the lacos script (van Dokkum 2001), and then the sky background was subtracted. Finally, the spectra were flux calibrated and 1D spectra were extracted using twodspec through the traces of the spectrophotometric standards, given that, at the M 31 distance, PNe are spatially unresolved. As an example, two fully calibrated 1D spectra are shown in Fig. 2.

Fig. 2. GTC OSIRIS 1D spectra of two of the sampled PNe plotted on the same scale and showing the best and worst S/N achieved in key diagnostic faint lines. Top: spectrum of the bright M50* PN. Bottom: spectrum of the control M1074 PN.

3. PN physical conditions and chemical abundances

3.1. Emission line fluxes and reddening measurements

The emission line fluxes and the corresponding statistical errors were measured from the extracted 1D spectra by multi-Gaussian fit using the splot task in IRAF. The fluxes need to be corrected for interstellar extinction, which is one of the primary sources of uncertainty linked to the PNLF. Extinction comes from different contributions, which cannot always be determined empirically. In general, the interstellar extinction all along the line-of-sight up to each M 31 PN CS is composed of both galactic and local contributions due to (i) dust in our own Galaxy that may change significantly with the position on M 31, (ii) dust in M 31 itself, all along the path up to the PN, (iii) circumnebular dust (mainly from mass ejected during the early AGB phase), (iv) intranebular dust, mixed with the ionized gas inside the PN shell, and (v) circumstellar dust (i.e., that located around and close to the CS). By measuring the extinction from the logarithmic extinction parameter, c(Hβ), contributions (i) to (iv) are accounted for, providing the extinction to each PN.

The c(Hβ) value was derived for each PN by comparing the observed and theoretical Hα/Hβ ratio (Storey & Hummer 1995) adopting the intrinsic ratio I(Hα)/I(Hβ) of 2.86 corresponding to T_e = 10 000 K and n_e = 1000 cm⁻³. We used the extinction curve from Cardelli et al. (1989) with R_v = 3.1. The observed flux ratio, extinction coefficient c(Hβ), de-reddened magnitude (m₅₀₀₇), and corresponding absolute magnitude (M₅₀₀₇) for each PN in the sample, are presented in Table 3. The mean value found, ⟨c(Hβ)⟩ = 0.13, and its range, 0.00–0.23 mag, are similar to those of PNe located at the outer disk and external regions of M 31 by Kwitter et al. (2012, ⟨c(Hβ)⟩ = 0.15; range from 0.01 to 0.36), and Corradi et al. (2015, ⟨c(Hβ)⟩ = 0.16; range from 0.06 to 0.29), indicating a dominant interstellar, not local, origin.

In the cases where line fluxes are contaminated by almost coincident He I, He II or H I lines, and these contaminants account for less than 2% of the flux, they are considered negligible. However, in some cases, for example, the [Ar IV] λ4711 line, blended with the He I 4713 line, or [Ne III] 3968, blended with He Iλ3965 and H Iλ3970 lines, the flux contribution of the contaminant lines can reach more than 30% for some of the PNe. The contribution of H I and He I lines to the blends was therefore computed from theoretical ratios with well isolated H I and He I lines, using the effective recombination coefficients shown in Table 4 by adopting a first estimation of the electron temperature and density of each PN.

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. The error introduced by the extinction correction is also propagated. Possible errors coming from the decontamination of blended lines are not included. Measured and de-reddened line fluxes normalized to Hβ = 100 with their corresponding uncertainties are presented in the appendix (Table A.1). The [O III] λ5007 line is saturated in the four brightest objects. The adopted fluxes included in that Table have been obtained from the theoretical relation between the [O III] λ5007 and λ4959 lines (2.98, according to Storey & Zeippen 2000).

3.2. Physical conditions

Electron temperature and density, T_e and n_e, can be determined in all nebulae from collisionally excited lines (CELs) for both high and low-excitation regions of the nebulae. Physical conditions were computed with PYNEB (Luridiana et al. 2015), version 1.1.7. Diagnostic errors are computed via Monte Carlo simulations with 1500 random values generated for each line flux using a Gaussian distribution centered in the observed intensity. The selected number of random values was established when the errors converged.

Given the available diagnostic lines, two ionization zones, with different temperatures and densities (see, e.g., García-Rojas et al. 2018) are assumed: a low-ionization region characterized by an ionization potential IP ≤ 17 eV, and a high-ionization region with IP > 17 eV. We consider T_e derived from the [O III] λ4363/λ4959 ratio, T_e([O III]), representative of the high-ionization zone, while T_e determined using the [N II] λ5755/λ6548 ratio, T_e([N II]), would be representative of the low-ionization region. In the case of the electron density, n_e obtained from the [S II] 6716/6731 lines corresponds to the low-ionization zone, while for the high-ionization region the available diagnostic ratio is [Ar IV] 4711/4740. The resulting temperatures and densities are shown in Table 5. Diagnostic T_e versus n_e diagrams are shown in Figs. A.1 and A.2.

While the T_e([O III]) determination is accurate in all PNe, T_e([N II]) has in some cases significant errors. An extensive empirical study of the temperature structure in PNe made by Kaler (1986) showed that T_e measured for low-ionization species in an homogeneous nebula were similar or lower than the temperature for the high-ionization zone. However, this is not always obtained in practice, as it occurs in some of our PNe in Table 5, where T_e([N II]) > T_e([O III]). Previous works (e.g., Corradi et al. 2015; Fang et al. 2018) obviate this issue by adopting a constant T_e([N II]) of 10 300 K (cf. Kaler 1986; Magrini et al. 2009) when He IIλ4686 line is detected (as is true for all our sample). In a recent review, Morisset (2017) showed that high density clumps in PNe can lead to a possible overestimation of the [N II] temperature due to an incorrect treatment of the collisional de-excitation of the upper level where the auroral [N II] λ5755 line arises. In our work the T_e([O III]) has been assumed to be the representative of the electron temperature for the whole nebula if uncertainties in T_e([N II]) are higher than ∼2000 K (that is the case for M1687* and M1074). In all other objects, both T_e([O III]) and T_e([N II]) have been used for the temperatures of the corresponding zones, even if T_e([N II]) > T_e([O III]).

In the case of the density, a two-zone scheme was also adopted as the more realistic situation, and both n_e from [Ar IV] and [S II] line ratios have been used in all the cases except in PNe M2068 and M1675, where the [S II] density diagnostic is saturated and hence, n_e from [Ar IV] diagnostic is adopted.

Although density estimations using [Cl III] 5518/5538 ratio were also obtained for three objects, we opted for not using them given their large uncertainties (see Table 5). We emphasize that only PNe M1074 and M2860 seem to show well-differentiated density zones, while for the remaining PNe, both density diagnostics show similar values within the uncertainties.

3.3. Ionic and total abundances

Ionic and total abundances were computed using PYNEB from the fluxes reported in Table A.1 and the adopted T_e and n_e in Table 5. We used the combination of atomic data sets shown in Table 6 for CELs and Table 4 for recombination lines. Errors in line fluxes and physical conditions have been propagated via Monte Carlo simulations using 1500 random values.

Ionic abundances for different ions of O, He, N, S, Ne, Ar and Cl are listed in Table A.2. According to the two-zone scheme described in the previous section, ionic abundances of O⁺, N⁺ and S⁺ were computed for the low-ionization zone, and abundances of He⁺, He²⁺, O²⁺, Ne²⁺, Ar²⁺, Ar³⁺, Ar⁴⁺, S²⁺, and Cl²⁺ were computed adopting the measured physical conditions for the high-ionization zone. Table A.2 presents the adopted ionic abundances considered as the most reliable measure of the ionic abundances. They avoid auroral lines, as well as lines possibly affected by telluric absorption, He or H contamination, or cases with larger errors.

To compute total elemental abundances, the undetected ionization stages of the different elements are corrected by adopting the usual scheme of the ionization correction factors (ICFs). Those proposed by Delgado-Inglada et al. (2014) are used for all the elements with the exception of N, for which the classical relation N/O = N⁺/O⁺ is applied (Peimbert & Costero 1969; Kingsburgh & Barlow 1994), following the recommendations of Delgado-Inglada et al. (2015). Given the importance of this particular ICF(N) we made the exercise of computing it using Delgado-Inglada et al. (2014) recipe and found that differences were below 22% in all the objects except for M2860 and M1675, which are very uncertain since their high O²⁺/O ratio fall outside the validity range for this ICF. Moreover, using Machine Learning techniques (C. Morisset, priv. comm.) we estimated the ICFs using the same set of photoionization models as Delgado-Inglada et al. (2014) and considering as an additional constraint the ratio Ar²⁺/Ar³⁺. In this case, the differences with the classical ICF(N) scheme stay below ∼25% except, again, in M2860 and M1675, where differences are much larger, and are also discrepant with the values derived using Delgado-Inglada et al. (2014) prescription. Unfortunately, at this moment there is no alternative way to determine a reliable ICF(N) for objects with very large O²⁺/O, and we decided to retain the same ICF(N) from the classical relation N/O = N⁺/O⁺ for all objects.

The computed ICFs are presented in Table 7, and elemental abundances are shown in Table 8. We also include, for comparison, present Solar System abundances from Lodders (2021).

3.4. Comparison with previous results

One PN from our bright sample (M2538*) was observed by Corradi et al. (2015) using the 3.5m APO-ARC telescope, whereas four PNe in the control sample (M1596, M2471, M2860, and M1074) were observed with the 3.5 m APO-ARC and 8.1 m Gemini-N telescopes by Kwitter et al. (2012). The latter correspond to their PN10, PN6, PN15, and PN16, respectively (we note that M1596 was misidentified with M1583: they interchanged PN9 and PN10 IDs). The quality of our spectra supersedes those obtained in previous references, giving more precise physical conditions and chemical abundances. In general, differences in the de-reddened fluxes for the same object in bright and isolated nebular lines as the [O III]λ4959 are no bigger than 3%, but in the case of the faint auroral line [O III] λ4363 the differences can amount up to 15 − 30%. The differences in the determined physical conditions mainly affect chemical abundances that use ICFs with a strong dependence on the determined ionization degree; this is the case of, for example, N/H, and this is discussed in more detail in Sect. 4.3.

4. Chemical abundances in an extended sample of PNe in M 31: Correlations among abundance ratios

To increase the significance of the comparison between PNe at the PNLF tip with less bright PNe (the control sample), an extended sample of PNe in M 31 from previous studies is also included. Data for 21 PNe with high-quality GTC spectra are taken from Balick et al. (2013), Corradi et al. (2015), and Fang et al. (2015, 2018)². Extinction, electron density and temperature, and chemical abundances of helium, oxygen, and nitrogen have been recalculated from the published line fluxes, in order to have a complete set of homogeneous chemical abundances of those key elements. Abundances for the rest of the elements have been taken straight from the original papers, when available. We used for all the objects the Cardelli et al. (1989) extinction law, the density n_e([S II]), and the ICF for nitrogen given by the relation N/O = N⁺/O⁺. When available, T_e[O III] and T_e[N II] are adopted for the high- and low-ionization species, respectively. Otherwise, T_e[O III] has been used for the whole nebula.

We note that different subsamples were used when analyzing the elemental abundances. In the case of duplicated objects, we only used the results from the present work. PN M2507 from Corradi et al. (2015) and PN14 from Fang et al. (2018) are the same object, and we adopted the abundances by Fang et al. (2018) owing to the fact that the extinction law applied by these authors is the same as ours. PN6 from Fang et al. (2015) and M 31-372 from Corradi et al. (2015) are excluded due to the lack of diagnostic sulfur lines. Finally, PN5 from Fang et al. (2018) was excluded due to the lack of detection of the [O II] λ7320 and λ7330 lines, which precludes proper determination of the nitrogen ICF.

4.1. Alpha elements

All three PN samples show a strong positive correlation of neon, argon and sulfur abundances with oxygen, as shown in Fig. 3 and Table 9. No significant differences are found among the bright, control, or extended samples, except for a slightly larger dispersion for argon in the extended sample – but we note that the argon ICF was reported by Delgado-Inglada et al. (2014) as the most uncertain of their ICFs. A tight correlation is expected for these four α-elements because they should vary in lockstep after being produced by the same nucleosynthesis processes in massive stars.

Fig. 3. Neon (a), Argon (b), and Sulfur (c) vs. Oxygen abundances for M 31 PNe. Green squares correspond to our bright sample and magenta circles to our control sample. Grey circles show the results from the extended sample of M 31 PNe from Balick et al. (2013), Corradi et al. (2015), Fang et al. (2015), and Fang et al. (2018). The grey cross represents the typical uncertainty of the extended sample. Blue triangles in (c) show data for M 31 H II regions from Zurita & Bresolin (2012).

The relation between sulfur and oxygen deserves further comment. Figure 3c shows that correlation, including also a sample of M 31 H II regions from Zurita & Bresolin (2012).

As with Ne and Ar, neither S nor O are expected to be significantly produced or destroyed during the evolution of low-mass PN progenitors, at least not at the high metallicities of the M 31 disk (cf. Shingles & Karakas 2013), but the available data compiled along the years have shown inconclusive results³. Figure 3c shows clearly for the first time the predicted strong S versus O abundance correlation for any PN sample (correlation coefficient r = 0.8; Table 9). Comparison can be made, for example, with Fig. 8 from Shaw et al. (2010), Fig. 1 from Henry et al. (2012), Fig. 6 from Maciel et al. (2017), or Fig. 10 from Fang et al. (2018), where a rough tendency of increasing S versus increasing O abundances for the Milky Way and M 31 PNe is apparent, but not quantified. This result gives us confidence on the quality of the observations and the analysis done. On the other hand, the whole sample of PNe is located in Fig. 3c just below the area occupied by the H II regions of Zurita & Bresolin (2012). We recover here the so-called sulfur anomaly, that is, sulfur abundances in PNe are systematically lower, for a given oxygen abundance, than those found in other probes of interstellar abundances. Henry et al. (2012) suggested that this “anomaly” could be caused by an incorrect computation of the contribution of ionization stages higher than S²⁺ by commonly used ICF schemes. However, new ICF schemes making use of a large database of photoionization models cannot solve this problem (Delgado-Inglada et al. 2014; Delgado-Inglada 2017).

4.2. Nitrogen, helium, and oxygen

Nitrogen and helium, as well as carbon and s-elements, are substantially synthesized during the evolution of PN progenitors, whereas oxygen, as mentioned earlier, is not. Comparison of N, He, and O relative abundances has been carried out along the years to test nucleosynthesis models and study the progenitors’ evolution. Nitrogen, in particular, is a key element in that analysis, but its complex nucleosynthesis makes interpretation intricate.

The relations among N, He, and O abundances for the bright, control, and extended samples are plotted in Fig. 4. A remarkable result is that, while the nine PNe of the bright and control samples are strongly correlated in the three diagrams, (a)–(c), the extended sample of 21 PNe from literature is not. Bright and control PNe appear intertwined in the diagrams and both samples have virtually the same fitting slopes and dispersion.

Fig. 4. Nitrogen vs. Oxygen (a), Nitrogen vs. Helium (b), and Helium vs. Oxygen (c). Labels as in Fig. 3.

A mild correlation between N and O was reported early by Pottasch (1983) for Galactic PNe, but only for low (12+log(N/H) < 7.8) N abundances. Henry (1990) also finds a positive correlation (r = 0.61) for his non-Type I Galactic PN sample. Kingsburgh & Barlow (1994), however, did not find any correlation, and conclude that no significant production of O via the ON cycle occurs during the second dredge-up, even for Type I PNe. More recently, Maciel et al. (2017) claimed a positive correlation between N and O for a large heterogeneous compilation of PNe in the Milky Way and other galaxies, although the correlation is indeed weak (their Fig. 12) and not quantified. On the theoretical side, it is expected that, for low-intermediate mass progenitors with near solar abundances, the secondary production of N comes at the expenses of C via the carbon-nitrogen-oxygen cycle. An increase in N abundance with metallicity is therefore expected (cf. Magrini et al. 2018). Figure 4a shows a strong correlation for our PN sample but a mild one for the extended sample, the latter with a Pearson r coefficient similar to that found by Henry (1990).

Identical behavior is shown for He and N in Fig. 4b: a strong correlation for the bright and main PN control samples, but a much weaker one, if any, for the extended sample (Table 9). As both N and He are produced by the same nucleosynthesis processes, the observed correlation in our data is expected (while the lack of it for the extended sample is striking).

Data in Fig. 4c are even more extreme, showing disparate correlation coefficients r = 0.8 and 0.2 for our nine-PN set and for the extended sample PNe, respectively (Table 9). Neither Henry (1990) nor Kingsburgh & Barlow (1994) found any correlation of He versus O in their Galactic PN samples. In the most recent analysis, Henry et al. (2018) compare the observed He and O abundances with the MONASH (Karakas & Lugaro 2016) and LPCODE (Miller Bertolami 2016) models, which indeed predict a modest enhancement of He with metallicity. Henry et al. (2018) argue that observational uncertainties of He/H likely obscure this theoretically predicted trend. Our 12+log(He/H) error bars are typically ±0.02 dex (Table 8), or ±0.005 in He/H, just at the level suggested by Henry et al. (2018) to allow meaningful comparison of models with observational data. Although the range in He and metallicity of our PNe in M 31 is certainly limited, and the slope for the fitting line in Fig. 4c is somewhat larger than model results, they do show for the first time the positive He versus O correlation predicted by the models.

4.3. The nitrogen-to-oxygen ratio

The N/O ratio has been extensively used to confront models with observations since it minimizes (together with the C/O ratio) possible effects due to different initial composition of the progenitor stars (cf. Henry et al. 2018; García-Rojas et al. 2018). Furthermore, N and O are synthesized by progenitors with very different mass ranges that evolve with different enrichment timescales (see Vincenzo et al. 2016; Esteban et al. 2020). For this reason, interpreting N/O data has always been a challenge.

Data for log(N/O) versus 12+log(N/H), 12+log(He/H), and 12+log(O/H) relations are plotted in Figs. 5a–c, respectively. We emphasize that having adopted the ICF(N) from the relation N/O = N⁺/O⁺, the derived N/O ratios for all PNe are insensitive to the particular ICF(N) values (Table 7). This is shown in Fig. 5a where ICF(N) values are labeled for the bright and control sets. As expected, a tight correlation between log(N/O) and 12+log(N/H) is observed in Fig. 5a for all PN samples, indicating that the variation in N/O is mainly due to variations in nitrogen.

Fig. 5. log(N/O) vs. Nitrogen (a), Helium (b), and Oxygen (c) abundances. Labels as in Fig. 3. Nitrogen ICF values are labelled in (a) (see text). The dashed rectangle in (b) defines the Type I PNe region from Peimbert & Torres-Peimbert (1983), whereas the horizontal dashed line sets the Type I PNe limit from Kingsburgh & Barlow (1994).

In Fig. 5b, which displays log(N/O) versus 12+log(He/H), we note that there is no object located in the region of type I PNe following the original definition by Peimbert & Torres-Peimbert (1983) (He/H > 0.125 and log(N/O) > −0.3). One of the control sample PNe, M1675, is located at the N/O ≥ 0.8 modified Type I zone defined by Kingsburgh & Barlow (1994), but its relatively low He/H abundance casts doubts about being a true Type I nebula. Figure 5b provides a first indication about the masses of the progenitor stars: The brightest PNe in M 31 are not Type I objects, that is, according to current single-star nucleosynthesis models, they do not come from massive (M_i ≳ 2.5 M_⊙) progenitor stars (cf. Phillips 2001).

A strong positive correlation for the bright and control samples (r = 0.8) is apparent in Fig. 5b. Such a correlation is expected since N and O are synthesized independently, while N and He are enhanced during the evolution of all progenitors by the same nucleosynthesis processes. Rough correlations between N/O and He in samples of PNe from the Galaxy, M 31 and the Magellanic Clouds were also found by Kaler (1979), Jacoby & Ciardullo (1999), Shaw et al. (2010), Maciel et al. (2017), and Henry et al. (2018) (but not by Kingsburgh & Barlow 1994).

Figure 5c shows a strong positive correlation (r ∼ 0.8) between N/O and O for the bright and control samples. However, as in the previous case, it vanishes in the extended sample.

In the case of H II regions and H II galaxies, the primary/secondary production of N is reflected in the observed behavior of N/O versus O: A plateau is generally found at low metallicity followed by a linear rise beyond certain metallicity (Henry et al. 2000; Magrini et al. 2018). However, in the Milky Way and M 31 only a scattered plateau of H II regions is apparent up to the highest measured metallicity (Esteban et al. 2020; Arellano-Córdova et al. 2021). For PNe, the situation is even less clear: a poor, if any, correlation between N/O and O has been reported. Data from Henry (1990), Kingsburgh & Barlow (1994), and Henry et al. (2018) show a large dispersion with, at best, a rough increase (or decrease: cf. Jacoby & Ciardullo 1999; Maciel et al. 2009; Milingo et al. 2010) of N/O with increasing O abundance.

In PNe with solar-like metallicity and low-mass progenitors, the correlation between N/O and O abundance seen in Fig. 5c for the bright and control samples is expected since N is mostly synthesized at the expense of the O and C already present in the star. Furthermore, a slope of N/O versus O steeper for PNe than for H II regions is also expected since, in a simple closed-box model, the latter depends only on the secular N and O enrichment of the interstellar medium up to the present time. In contrast, the slope for PNe depends on (i) the secular N and O enrichment only up to the epoch of the formation of the progenitors (the O abundance, in the abscissa), plus (ii) the extra production of N during the evolution of the progenitors (increasing the N/O ratio in the ordinate). Our data indicate a slope around 70° for the bright and control PN samples, larger than typical values of around 45° from extragalactic H II regions at similar O/H abundances (Henry et al. 2018).

Finally, we compare in Fig. 6 the observed He/O and N/O ratios with model results from Karakas & Lugaro (2016) (MONASH) and Ventura et al. (2018) (ATON). Metallicities suitable for M 31 (i.e., embracing the solar value Z = 0.014; Saglia et al. 2018) have been chosen. Once more, the PNe of the bright and control samples appear arranged in a clear anticorrelation, r = −0.65, while the extended control sample does not. MONASH results for Z = 0.007 and Z = 0.014 adequately encompass the observed values of He/O for our nine bright PNe and most of the extended sample, although seven PNe from the latter lie in the zone of lower metallicity, with log(He/O) > 2.6. None of those PNe are He overabundant Type-I objects, so their positions reflect their lowest O abundance, which is consistent with the fact that five out of those seven outliers are located in areas belonging to the M 31 diffuse halo and outer structures. More interesting is the behavior of the N/O ratio. As will be discussed in the next section, the estimated progenitor masses for our nine PNe are restricted to a very narrow range, 1.1 → 1.8 M_⊙, and neither MONASH nor ATON models reach the high values of log(N/O) ≳ −0.4 observed in five of our PNe for that range of masses. A similar discrepancy has been reported by, among others, Henry et al. (2018) and points to either a serious limitation of the models (with rotation, magnetic field or extra-mixing not adequately taken into account) or to a wrong determination of the initial/final masses via evolutionary models and the Hertzsprung-Russell (HR) diagram.

Fig. 6. log(N/O) vs. log(He/O) for the M 31 PN samples together with MONASH (Karakas & Lugaro 2016) and ATON (Ventura et al. 2018) models, labeled with the initial mass of the progenitor. Symbols are the same as in Fig. 3.

In summary, a key characteristic of our data in Figs. 4a–c, 5b, c, and 6 is the higher degree of correlation of various elemental abundances for the bright and control samples with respect to the extended sample. Apart from the smaller uncertainties, our bright and control samples represent an homogeneous population of disk PNe of M 31, whereas the extended sample includes a broader range of populations, with distances ranging from d_disk ∼ 5 to more than 100 kpc and belonging to the halo or to its complex outer substructures, which span a wide range in metallicity covering roughly one order of magnitude (Bhattacharya et al. 2021). The location of PNe in these diagrams, and even the PNLF of each population (Bhattacharya et al. 2021), is sensitive to the initial chemical composition, star formation history, and subsequent interaction of each substructure with M 31, and that would explain the observed dispersion of the extended sample in the figures.

5. PN central stars and their progenitors

5.1. Cloudy models

The effective temperature, T_eff, and luminosity, L_*, of the PN central star (CSPN) locate it in the HR diagram and measure, through the available post-AGB tracks, the final mass M_f of the progenitor star. Then, an initial mass M_i can be derived using available semiempirical IFMRs.

A grid of models was run for each PN using the pyCloudy package (Morisset 2013) based on the 1D photoionization code Cloudy 17.01 version (Ferland et al. 2017). As starting values for T_eff and L_*, we adopted the values estimated with the method outlined in Dopita & Meatheringham (1991), using the PN excitation class (EC) and its Hβ absolute flux. The latter was derived by multiplying the de-reddened spectroscopic I(Hβ)/I(λ5007) flux ratios in Table A.1 by the I(λ5007)_M absolute fluxes from Merrett et al. (2006). For the four brightest PNe, whose λ5007 line is saturated, I(λ5007) was calculated to be 2.98 × I(λ4959) (Storey & Zeippen 2000). Excitation class values were obtained for the nine bright and control sample PNe from the I(λ4686)/I(Hβ) line ratio (Dopita & Meatheringham 1991). They show a remarkably narrow range of intermediate ECs, between 4.5 and 6.3.

The aim of the Cloudy modeling is to match a selection of observed nebular line ratios with those predicted by the model through a given ionization source and a number of assumptions necessitated by the lack of information (see, e.g., Jacoby & Ciardullo 1999; Magrini et al. 2004; Kwitter et al. 2012). Cloudy also requires as an input the energy distribution of the CSPN. We used the latest Rauch H-Ni atmospheres with log g = 6.5 or, alternatively, Rauch H-Ca when an effective temperature higher than 190 000 K is required (Rauch 2003). Chemical abundances in Table 8 are adopted, and all PNe have been assumed to be optically thick. The nebular geometry is assumed to be spherical with an inner radius of 10¹⁷ cm, and the density is set as constant. Even so, the effect of different density distributions has been explored, in particular for cases such as M1687* or M1074 where different values for the density when using low- or high-excitation ions have been obtained.

The model selection procedure is based mainly on a positive correlation of the He II 4686 and He I 5876 fluxes with the stellar effective temperature, T_eff, whereas the absolute fluxes of Hβ and [O III] λ5007 constrain the stellar luminosity, L_*. Additional line ratios, similar to those investigated in Kwitter et al. (2012), are also taken into consideration, and further nebular parameters, such as output radii and nebular masses, were always checked to be kept physically reasonable.

We note that the abundance of carbon, a major nebular coolant, is actually unknown for the PNe in our sample, and this may represent an important source of uncertainty in the modeling. For this reason, values for the carbon abundance corresponding to C/O ratios of 0.5, 1.0, 1.5, and 2.0 have been tested. However, changes in C abundance result in only small variations in the He II 4686, He I 5876, [O III] 5007, and Hβ diagnostic fluxes, and we finally set a C/O ratio equal to 1.0. In short, although there is no information on the C abundance for our objects, the main emission-line strengths predicted by the models are in good agreement with those observed, and we can be confident of the robustness of the procedure and the resulting stellar parameters.

The T_eff and L_* from Cloudy modeling are listed in Table 10. The most remarkable result is the very narrow ranges exhibited by both stellar parameters: log(L_*/L_⊙) = 3.57 ± 0.07, log T_eff = 5.11 ± 0.07. All nine CSPNe seem to be caught in a very similar post-AGB evolutionary stage and be produced by similar progenitor stars. Also, as expected, the luminosities L_* of the bright PNe are larger (∼30% on average) than those for the control sample, whereas temperatures T_eff do not show any obvious difference.

5.2. Location in the HR diagram, progenitor masses, and post-AGB model predictions

The positions in the HR diagram of our nine CSPNe are shown in Fig. 7, together with post-AGB evolutionary tracks from Miller Bertolami (2016) for Z = 0.02. Differences between these tracks and the classical tracks of Vassiliadis & Wood (1994) – basically, somewhat smaller final masses and substantially shorter, accelerated evolutionary times for the former – have been discussed in recent papers (cf. Miller Bertolami 2016; Fang et al. 2018; Henry et al. 2018). Interpolation of the evolutionary tracks allows us to estimate M_f for each CSPN, and derive tentative masses M_i of their progenitors. The uncertainty associated with the interpolation is expected to be ∼10 − 15% for accurately derived CSPN parameters (Serenelli et al. 2021).

Fig. 7. HR diagram showing the location of the CSs of our PN samples as derived from Cloudy modeling in the Miller Bertolami (2016) model tracks. The final and corresponding initial masses are indicated on the tracks.

The resulting range of CSPN masses for the four brightest PNe in M 31 (Table 10) is remarkably narrow, M_f = 0.568 to 0.578 M_⊙, with a mean value of ⟨M_f⟩ = 0.574 ± 0.004 M_⊙. That is in excellent agreement with the maximum final mass derived by Valenzuela et al. (2019) from their PNLF simulations, M_f = 0.58 M_⊙. The five PNe from the control sample show, in turn, a slightly broader range, M_f = 0.545 to 0.578 M_⊙ with a mean value of ⟨M_f⟩ = 0.566 ± 0.016 M_⊙.

This analysis shows that all nine PNe, including the brightest ones, originate from low-intermediate mass stars with M_i < 2 M_⊙. The four brightest PNe (and also, to some extent, the less-bright control PNe) appear in Fig. 7 tightly clustered around the turning point of the tracks, close to the “knee” where T_eff and L_* values simultaneously reach their maxima⁴.

6. Discussion

Understanding the [O III] λ5007 emission from the nebular shell as the CSPN quickly evolves across the HR diagram and, therefore, the origin of the observed PNLF is not trivial. The line emission strength strongly depends on the density, ionization, and chemical properties of the nebula, combined with the evolution of the post-AGB CS (Schönberner et al. 2007). In this section we analyze the properties of our targets trying to understand the nature of the brightest PNe as compared to the control sample, and to the predictions of PNLF theoretical models.

6.1. Physical properties of the most luminous PNe and their central stars

Summarizing our observational results, we find that the nebular electron temperatures for the bright and control samples do not show significant differences, whereas electron densities appear larger for the former (Table 5). The n_e weighted average of the bright sample is ∼40% and ∼18% higher than in the control sample when measured from high- and low-ionization diagnostic lines, respectively.

No significant differences are found for the masses and effective temperatures of the progenitor stars in the bright and control samples. They are in the range M_i = 1.1 to 1.8 M_⊙, and T_eff = 106 000 to 166 000 K. All four bright PNe are located before the theoretical isochrone for 5000 a, while two PNe of the control sample appear slightly older, between 5000 and 10 000 a.

Marigo et al. (2004) studied the evolution of CSPNe by modeling the observed PNLF and its bright cutoff. In their models, the most luminous PNe are those originated from intermediate mass stars with final masses M_f = 0.7 → 0.75 M_⊙, setting a lower limit for the inital mass M_i > 2.5 M_⊙. They also claimed that the maximum [O III] λ5007 emission would be attained by nebulae at the thick/thin transition during the CSPN post-AGB evolution.

On the other hand, Schönberner et al. (2007) presented a 1-D, spherically symmetric, radiative-hydrodynamics modeling of CSPNe with masses from M_f = 0.565 to 0.696 M_⊙. They concluded that, at odds with Marigo et al. (2004), the conversion of stellar UV radiation into line emission occurs most efficiently only for optically thick shells, attaining the observed PNLF [O III] λ5007 cutoff at M^* = −4.5 mag for a stellar mass of 0.62 M_⊙, provided the PN remains optically thick.

As discussed in Schönberner et al. (2007), and more recently in Valenzuela et al. (2019), the question of the actual optical depth (thin or thick) to Lyman continuum photons for the brightest PNe is paramount to understand the PNLF, since the efficiency for transforming stellar UV photons into nebular ones is governed by the changes in optical depth during the coupled CS and PN evolution. All nine of our PNe show relatively high electron densities and CS effective temperatures, together with moderate values of the excitation parameter EC ≲ 6, pointing to optically thick PN shells. This is also supported by the presence of the [O I] λ6300 line in all objects. Although the morphologies of our PNe are actually unknown, any departure from sphericity would also help to maintain the PNe opaque for longer during their evolution. Three-dimensional modeling of PNe comprising dense equatorial tori and thin lobes is found to reproduce well the observed line emissivities including those of [O III] (Gesicki et al. 2016).

Recently, Gesicki et al. (2018) reevaluated the PNLF using the Miller Bertolami (2016) evolutionary tracks. Seven CSPNe tracks with M_i = 1.0 − 3.0 M_⊙ (M_f = 0.532 − 0.706 M_⊙), and the Torun photoionization code for the PN nebular emission were assumed. They showed that, for the five models with M_i = 1.1 − 2.5 M_⊙, the [O III] emission reaches values of m₅₀₀₇ near the PNLF cutoff M^* during hundreds to thousands of years, while the nebulae remain opaque. Assuming a Salpeter initial mass function and four star formation histories representative of different types of galaxies, synthetic PNLFs were built where the observed bright cutoff is shown to be populated by progenitors with moderate post-AGB stellar luminosities log(L_*/L_⊙) = 3.75 ± 0.13. This value is in agreement with, although slightly larger than, log(L_*/L_⊙) = 3.63 ± 0.03 found for the four brightest sources in our sample (Table 10). Concerning the CSPN progenitor masses at the bright PNLF end, Gesicki et al. (2018) estimate a range M_i = 1.1 − 2.0 M_⊙, which is similar to that of our full sample, M_i = 1.1 − 1.8 M_⊙. But we note that our four brightest PNe display an even narrower range of masses, with a mean value ⟨M_i⟩ = 1.45 ± 0.09 M_⊙. According to Gesicki et al. (2018), such stars would remain above the bright-end cutoff M^* for approximately one thousand years.

In summary, while the Marigo et al. (2004) and Schönberner et al. (2007) models predict CS masses substantially larger than what we estimated, the modeling by Gesicki et al. (2018) using post-AGB tracks from Miller Bertolami (2016) agrees with our derived values for both luminosities and progenitor masses.

For completion, it should be mentioned that Davis et al. (2018) studied 23 PNe in the bulge of M 31 within 1 mag of the PNLF cutoff, finding substantial levels of extinction, c(Hβ) ∼ 0.6, most of which of circumstellar origin. This implies the presence of intrinsically overluminous objects, one-half to one magnitude brighter than the brightest PNe in our M 31 disk sample. They would require, even using the evolutionary tracks from Miller Bertolami (2016), CSPN masses > 0.66 M_⊙, that is, progenitors > 2.5 M_⊙, which are not expected to be found in a sufficient number in M 31’s bulge. Davis et al. (2018) suggest that objects populating the bright end of the PNLF in the bulge of M 31 might not be PNe at all, but nebulae around symbiotic stars, as previously proposed by Soker (2006). Even if this limitation in the number of stars with > 2.5 M_⊙ would not apply in the disk of M 31, the present study does not find any indication of such a high level of extinction for the target PNe. The c(Hβ) values range from 0 to 0.23 with, at most, ∼0.1 of circumstellar origin, resulting in the low CSPN masses calculated. We can also rule out in our sample indications of a symbiotic nature: all PNe are located well inside the PN area of the Gutierrez-Moreno et al. (1995) [O III] diagnostic diagram, and they show no signs of having bright and red CSs nor the high-excitation [O VI] line emission at 6825Å typical of symbiotic stars (cf. Mikołajewska et al. 2017).

6.2. Large N/O ratios found in PNe with low-mass progenitors

An intriguing question remains to be answered, namely, if progenitors of the brightest PNe are indeed M_i ∼ 1.5 M_⊙ stars, how can they produce the N overabundance found for at least five PNe in our sample (those with log(N/O) ≳ -0.4 in Fig. 6), as well as for many PNe in previous studies (cf. Henry et al. 2018; Fang et al. 2018). Our five sources show N/O values significantly larger than expected for their progenitor masses, predicted for M_i ∼ 3 − 3.5 M_⊙ stars according to models. Furthermore, data on our nine bright and control PNe in Fig. 6 suggest a smooth progression toward larger N/O values as He/O decreases, in contrast with the abrupt increase in N/O at M_i ∼ 3.5 − 4 M_⊙ predicted by the models.

The issue is related to the low-end mass limit for stars experiencing hot bottom burning (HBB) during the AGB. Some models set the enrichment threshold at M_i ≳ 4 M_⊙ (Karakas & Lugaro 2016; Cristallo et al. 2015; Marigo et al. 2017) while other models that include overshooting at the boundary core (Miller Bertolami 2016; Ventura et al. 2018), decrease it to M_i ∼ 3 M_⊙. The limit has been continuously pushed down by the observations: Davis et al. (2019) found a 3.4 M_⊙ CS of a Type I PN in an open cluster in M 31 where indeed HBB has taken place, whereas samples of PNe from the Milky Way (Henry et al. 2018) and M 31 (Fang et al. 2018) suggest the limit could be as low as M_i ∼ 2 M_⊙. Our data indicate even lower masses, M_i ≲ 1.8 M_⊙, for PNe with the larger values of N/O in our sample, and no extant model is able, to our knowledge, to reach near that value.

A possible way to circumvent this conundrum is to reevaluate the IFMR of low-intermediate mass stars, and a recent paper from Marigo et al. (2020) might represent an important step forward in that direction. They carefully determine initial and final masses for a sample of white dwarfs (WDs) in old Galactic open clusters, and find a kink at M_i = 1.65 − 2.10 M_⊙ in the IFMR, whose peak at M_i ∼ 1.8 − 1.9 M_⊙ corresponds to measured WD masses of M_f = 0.70 − 0.75 M_⊙. It is important to emphasize that those semiempirical data strongly disagree with current models of AGB evolution: CSPN with masses of M_f = 0.70 − 0.75 M_⊙ would originate from main sequence stars with much larger masses M_i ≳ 3 − 3.5 M_⊙ according to the Vassiliadis & Wood (1994) or Miller Bertolami (2016) models. Moreover, stars with masses around M_i ∼ 1.5 M_⊙, located at the lower-mass side of the IFMR kink, would evolve into relatively high-mass, M_f ∼ 0.6 M_⊙, WDs according to Marigo et al. (2020). That is larger than Miller Bertolami (2016) model tracks predict (a M_f ∼ 0.6 M_⊙ star would originate from a M_i ∼ 2.3 M_⊙ star according to his modeling), although we note that M_i = 1.5 M_⊙ progenitors do end as M_f ∼ 0.6 M_⊙ dwarfs in Vassiliadis & Wood (1994) tracks.

Determining accurate masses from WD stars and their progenitors is a difficult task, and the M_i ∼ 1.5 M_⊙ region of the IFMR from Marigo et al. (2020) is really based on data from a single open cluster, M67 (their Fig. 1). A recent reevaluation of seven WD masses belonging to that cluster by Canton et al. (2021) confirm the mean value of M_f = 0.60 ± 0.01 M_⊙ for the WDs, and M_i = 1.52 ± 0.04 M_⊙ for their progenitors (quoted uncertainties are formal errors of the mean, not dispersion of the sample).

Assuming that M_i ∼ 1.5 M_⊙ stars, like the proposed progenitors of our four brightest PNe, would end as M_f ∼ 0.6 M_⊙ CSPN, we wish to determine what effects this would have on the nucleosynthesis and the yields of elements such as He, C, N, and possibly O. We also want to know if the discrepancy between predicted and observed N/O ratios will be solved or alleviated, given the now more massive stellar cores of the AGB stars.

In the mass regime M_i ∼ 1.5 − 2.0 M_⊙ a process takes place that is key to understanding the AGB and post-AGB evolution: stars with solar metallicity and masses up to M_i ∼ 1.8 − 1.9 M_⊙ experience a core He flash at the beginning of their core He-burning phase, while those with M_i = 2.0 M_⊙ and larger do not. That defines the fundamental frontier between low- and intermediate-mass stars, and changes not only the evolutionary timescales, but also the strength of the 3rd dredge-up and the chemistry. As mentioned earlier, that mass limit at M_i ∼ 1.8 − 1.9 M_⊙ also coincides with the observed kink in the IFMR, making CSPN and WD stars that are more massive than previously thought (Marigo et al. 2020). On the other hand, according to Miller Bertolami (2016) modeling, it is just at M_i = 1.5 M_⊙ where the largest carbon enrichment is expected, almost double that of a M_i = 2.0 M_⊙ star (his Fig. 6). If such an M_i = 1.5 M_⊙ star were to end as an M_f = 0.60 M_⊙ CSPN its properties would be quite different from one with M_f = 0.576 − 0.580 M_⊙ in terms of 3rd dredge-up efficiency, initial H envelope, and H-burning rate, and it would have higher luminosity (Miller Bertolami 2016). However, in terms of the N surface enrichment, only at masses M_i > 3.0 M_⊙ does there appear a significant jump upward according to all extant models, and therefore the origin of the disagreement between modeled and observed N/O ratios remains unexplained.

Further efforts to improve low-mass star models, including extra-mixing processes beyond overshooting, such as rotation (meridional circulation and diffusion by shear turbulence), magnetic field, or thermohaline mixing, have been encouraged by better and better observational data in recent years (Delgado-Inglada et al. 2015; Henry et al. 2018; Fang et al. 2018). The substantial N overabundance measured in several of the most luminous PNe in M 31 provides motivation for these continuing efforts.

7. Conclusions

The PNLF is an astrophysical observable that appears to have a universal bright cutoff. This is thought to be caused by an unexpected constancy in the combined results from the complex physics of stellar interiors, stellar evolution all along the HR diagram, nucleosynthesis and dredge-up processes, the physics of envelope ejection, ionization, and the dynamical evolution of the PNe. In addition to having a mild dependence on metallicity, all these processes seem to “conspire” to yield an extremely simple result: All the brightest PNe in old and younger stellar systems have a unique cutoff magnitude at M^* = −4.54 ± 0.05.

While the building of the PNLF itself has been recently explained using updated stellar evolutionary tracks (Gesicki et al. 2018; but see Davis et al. 2018), the detailed nature of the PNe populating the bright end of the PNLF is still largely unknown due to distance uncertainties for Galactic PNe and the inherent faintness of PNe in external galaxies.

The goal of this work has been to analyze in detail a representative sample of bright PNe in the galaxy M 31, our closest Milky Way-like neighbor, which hosts the largest collection of PNe known in any galaxy. We carefully selected and observed two samples of “bright PNe" (four objects inside the bright cutoff bin) and “control PNe” (five objects typically 0.5 mag fainter). We supplemented them with an “extended sample” of another 21 PNe from previous works observed with the same telescope and instrument and reanalyzed here to ensure homogeneity. Relevant conclusions from this work are:

• 10.4 m GTC spectra obtained under excellent observing conditions allow us to reach the unprecedented depth of a 3σ rms signal at ∼0.2% of I(Hβ) in the blue, even in the faintest nebulae.

• These spectra are carefully and homogeneously analyzed using updated physics and the most recent ICFs for abundance calculations, minimizing the uncertainties in the derived physico-chemical properties.

• The four brightest and the five control PNe show a remarkable uniformity in nebular characteristics: Electron temperatures range between T_e([O III]) = 10 100 and 12 700 K. All nebulae show large electron densities, n_e([Ar IV]) = 5000–32 000 cm⁻³, those of the brightest nebulae being, on average, ∼40% larger than those of the control sample. The nebular ECs also span a narrow range of EC = 4.5–5.6 for the bright and control nebulae.

• Derived abundances for He, Ne, N, O, Ar, and S agree neatly with their predicted qualitative interrelations, and this is the first time such global agreement is observed in any PN sample, including those in the Milky Way. Namely, (i) three observed alpha-elements, Ne, Ar, and S, behave, as theoretically expected, in lockstep with respect to O, and the expected tight S/H versus O/H correlation is measured for the first time (while the sulfur anomaly persists). (ii) He and N appear strongly correlated with each other and with O, following their expected dependence on metallicity but, significantly, only for the bright and control, not the extended, samples. (iii) The expected He versus O, N/O versus He, and N/O versus O correlations are observed clearly for the first time, again only for the bright and control samples. (iv) The diagram N/O versus He/O is introduced and shows a tight anticorrelation (once more, only for the bright and control samples), allowing comparison with two families of models, MONASH (Karakas & Lugaro 2016) and ATON (Ventura et al. 2018). Whereas He/O ratios are consistent with model results, measured N/O values are 1.5–3 times larger than model predictions for the range of progenitor masses derived for five of our PNe. This, adding to previous cases from the literature, indicates an important limitation in the models. (v) The contrasting behavior of the bright and control samples, on the one hand, and the extended sample, on the other, is thought to be caused by the fact that the extended sample is an inhomogeneous collection of PNe that belong to different M 31 components (the extended halo or its outer substructures, the Warp, G1 Clump, Northern Spur, and Giant Stream), with different chemical compositions and star formation histories.

• A primary goal of this work was to determine the CS masses (and progenitor initial masses via an IFMR) of these bright PNe. Cloudy modeling indicates a remarkably homogeneous set of luminosities and effective temperatures, with average values and a 1σ dispersion of ⟨L_*/L_⊙⟩ = 4300 ± 310, T_eff = 122 000 ± 10 600 K for the CS of the bright sample, and ⟨L_*/L_⊙⟩ = 3300 ± 370, T_eff = 135 000 ± 26 000 K for those in the control sample.

• The CS of the four brightest PNe appear closely grouped at the “knee” of the M_i = 1.5 M_⊙ evolutionary track from Miller Bertolami (2016), whereas those of the five control PNe span a slightly broader range of initial masses, M_i = 1.1 → 1.8 M_⊙

Our overall conclusion is that the brightest PNe in the (extended) disk of M 31 are normal, moderately dense nebulae, with low-mass, M_i ∼ 1.5 M_⊙, progenitors, that, according to the modeling by Gesicki et al. (2018), are able to populate the bright end of the PNLF for ∼1000 years. We note that the stellar mass distribution in a star-forming environment such as the disk of M 31 spans a large range, including stars with higher masses that reach larger luminosities in their post-AGB evolution. We plan to extend this work to stellar systems with different star formation histories or metallicities, starting from the Large Magellanic Cloud, to gain more comprehensive observational constraints about PNe at the tip of the PNLF.

Acknowledgments

During the preparation of this article, while finishing her Ph.D. at the Instituto de Astrofísica de Canarias, Rebeca Galera-Rosillo untimely passed away at the age of 31. That had a devastating impact on her family, colleagues and friends. Her passion for astrophysics, her solidarity with the underprivileged and her never-ending smile will not be forgotten. AM, and RC heartily acknowledge the material, and above all, personal support provided by BB, KK, LM, and EV to Rebeca during her stays at their home institutions. Particular thanks are due to Christophe Morisset for his help with PyCloudy and Machine Learning programs to constrain the ICFs, and to Simone Madonna for many discussions and continued support to Rebeca. This work is based on observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. AM, RC, JG-R, and DJ acknowledge support from the State Research Agency (AEI) of the Spanish Ministry of Science, Innovation and Universities (MCIU) and the European Regional Development Fund (FEDER) under grant AYA2017-83383-P, and support under grant P/308614 financed by funds transferred from the Spanish Ministry of Science, Innovation and Universities, charged to the General State Budgets and with funds transferred from the General Budgets of the Autonomous Community of the Canary Islands by the MCIU. JG-R acknowledges support from the Severo-Ochoa excellence programs, SEV-2015-0548 and CEX2019-000920-S, DJ also acknowledges support from the Erasmus+ program of the European Union under grant number 2020-1-CZ01-KA203-078200, and EV support from Spanish grant PGC2018-101950-B-100. AM, RC and JG-R acknowledge support from the Canarian Agency for Research, Innovation and Information Society (ACI-ISI), of the Canary Islands Government, and the European Regional Development Fund (ERDF), under grant with reference ProID2021010074.

References

Appendix A: Supporting material

In this appendix we include the following material: plasma diagnostics for the four PNe in the bright sample (Fig. A.1); plasma diagnostics for the five PNe in the control sample (Fig. A.2); line identification and flux table for all nine PNe in our bright and control samples (Table A.1); and ionic abundances (Xⁱ⁺/H⁺) for the PNe in the bright and control samples (Table A.2).

Fig. A.1. Plasma diagnostics diagrams of the observed PNe bright tip sample.

Fig. A.2. Plasma diagnostics diagrams of the observed PNe control sample.

All Tables

All Figures

	Fig. 1. Location of the target PNe. North is up, and east is to the left. Green squares correspond to the brightest sample targets and magenta diamonds to the control sample. More details are given in the text.
In the text

	Fig. 2. GTC OSIRIS 1D spectra of two of the sampled PNe plotted on the same scale and showing the best and worst S/N achieved in key diagnostic faint lines. Top: spectrum of the bright M50* PN. Bottom: spectrum of the control M1074 PN.
In the text

Fig. 3. Neon (a), Argon (b), and Sulfur (c) vs. Oxygen abundances for M 31 PNe. Green squares correspond to our bright sample and magenta circles to our control sample. Grey circles show the results from the extended sample of M 31 PNe from Balick et al. (2013), Corradi et al. (2015), Fang et al. (2015), and Fang et al. (2018). The grey cross represents the typical uncertainty of the extended sample. Blue triangles in (c) show data for M 31 H II regions from Zurita & Bresolin (2012).

In the text

	Fig. 4. Nitrogen vs. Oxygen (a), Nitrogen vs. Helium (b), and Helium vs. Oxygen (c). Labels as in Fig. 3.
In the text

	Fig. 5. log(N/O) vs. Nitrogen (a), Helium (b), and Oxygen (c) abundances. Labels as in Fig. 3. Nitrogen ICF values are labelled in (a) (see text). The dashed rectangle in (b) defines the Type I PNe region from Peimbert & Torres-Peimbert (1983), whereas the horizontal dashed line sets the Type I PNe limit from Kingsburgh & Barlow (1994).
In the text

	Fig. 6. log(N/O) vs. log(He/O) for the M 31 PN samples together with MONASH (Karakas & Lugaro 2016) and ATON (Ventura et al. 2018) models, labeled with the initial mass of the progenitor. Symbols are the same as in Fig. 3.
In the text

	Fig. 7. HR diagram showing the location of the CSs of our PN samples as derived from Cloudy modeling in the Miller Bertolami (2016) model tracks. The final and corresponding initial masses are indicated on the tracks.
In the text

	Fig. A.1. Plasma diagnostics diagrams of the observed PNe bright tip sample.
In the text

	Fig. A.2. Plasma diagnostics diagrams of the observed PNe control sample.
In the text