The Critical Role of Dust On The [O III] Planetary Nebula Luminosity Function’s Bright-End Cutoff

Our study uses subsets of Large Magellanic Cloud and M31 PNe, whose nebular reddening and central star properties (bolometric luminosity and effective temperature) are relatively well-determined via photoionization models. To ensure self-consistency, we chose for analysis the largest and most uniform PN samples available. Although measurements exist from several smaller investigations, differences in survey depth, resolution, wavelength coverage, and reduction/analysis techniques would likely introduce unquantifiable systematics into our results. The samples described below therefore represent our attempt to achieve maximum uniformity over a sample sized large enough to draw robust conclusions.

Our LMC data come solely from the PN analyses of Dopita & Meatheringham (1991a, b) and Dopita et al. (1997), with the former two references obtaining their spectrophotometry (and nebular extinctions) from Meatheringham & Dopita (1991a, b). Since we are only concerned with the PNe that define the bright-end cutoff of the PNLF, we use the PNe’s observed [O III] $\lambda 5007$ magnitudes as measured by Jacoby et al. (1990) to restrict our analysis to objects within one magnitude of $M^{*}$. Fourteen PNe satisfy these criteria; these are listed in Table 1.

To explore a different PN population, we use similar measurements of M31 PNe from the spectrophotometric analyses of Kwitter et al. (2012), Balick et al. (2013), Fang et al. (2018), and Ueta & Otsuka (2022). (The latter is a reanalysis of data originally presented in Galera-Rosillo et al., 2022.) These works, which are all based on 8-m class telescope spectrophotometry, provide extinctions and central star properties for 28 PNe with [O III] luminosities within one magnitude of $M^{*}$. (Three additional [O III]-bright bulge PNe are available from the 4-m observations of Jacoby & Ciardullo (1999); however, since the errors on their extinction measurements are 3 to 5 times larger than those from the other studies, these objects have been excluded from the analysis.) Additionally, four of the Kwitter et al. (2012) PNe have also been observed and analyzed independently by Galera-Rosillo et al. (2022) and Ueta & Otsuka (2022). These repeat measurements provide some information about the uncertainties associated with the determinations of PN core mass and nebular extinction. The M31 PNe used in our analysis are also given in Table 1.

We note that the sets of PNe being analyzed do not constitute statistically complete samples. However, they are representative of the bright PNe present in both galaxies. In the LMC, the dataset contains 14 of the brightest 23 PNe in the galaxy, and an Anderson & Darling (1954) test confirms that the magnitude distribution of the objects in our sample is consistent with that of the full population of PNe within one magnitude of $M^{*}$. Similar agreement is found between the M31 PNe being analyzed and the full sample of $M_{5007}<M^{*}+1$ PNe cataloged by Merrett et al. (2006). Thus, the results found in our study should apply to the overall population of [O III]-bright PNe in the galaxies.

The $c_{{\rm H}\beta}$ values measured for both the LMC and M31 PNe reflect the total amount of extinction along the line-of-sight. This extinction comes from three components: foreground dust in the Milky Way, foreground dust within the program galaxies themselves, and circumnebular dust ejected by the PN’s progenitor.

To remove the first component, we used the extinction maps of Schlegel et al. (1998) to derive $E(B-V)=0.080$ for the LMC and $E(B-V)=0.062$ for M31. These reddenings are supported by the observed $c_{{\rm H}\beta}$ values of the galaxies’ PNe; when this component is removed, the floor of the two reddening distributions is close to zero.

Quantifying the amount of extinction caused by foreground dust within the target galaxies is more problematic. Yet it must be carefully considered: if this internal component is larger for PNe from high-mass progenitors, it could introduce a systematic error into our measurement of the extinction versus core-mass relation. To examine this possibility, we used the 3-D distribution of dust within the Milky Way, as measured by Guo et al. (2021). This model describes the vertical structure of dust in the solar neighborhood, i.e., at galactocentric radii between $\sim 6$ and $\sim 10$ kpc. If we assume that the dust distribution in the target galaxies is similar, then the model predicts that the maximum amount of foreground extinction a PNe could experience is $A_{5007}\sim 0.15\sec i$, where $i$ is the inclination of the galaxy.

Next, we assume that the vertical distributions of our program PNe are related to the ages of their progenitor stars. Mazzi et al. (2024) have used Gaia data to quantify the age-scale height relationship for Milky Way objects in the solar cylinder, and we adopt this same relation for our target galaxies.

Finally, we apply these assumptions to the target galaxies and estimate the amount of foreground dust that might attenuate PNe of a given age (or, equivalently, progenitor mass). This allows us to address the question of whether the foreground reddening of PNe observed in the top 1 mag of the PNLF can introduce a systematic term into our analysis of circumnebular extinction.

The results of our model show that the effect of galactic internal extinction can vary greatly from object to object, with a distribution that depends strongly on the age of the target population. For example, if a PN sub-population has a scale-height similar to that of the extinction, the distribution of their $A_{5007}$ values is roughly Gaussian in nature, with most objects having $A_{5007}$ values near the distribution’s mean. In contrast, if the scale height is much larger than the dust, the distribution of $A_{5007}$ values is almost flat and observations will record many more high- and low-extinction objects.

Yet, almost paradoxically, this scale-height dependence has almost no effect on the mean properties of the observed sample: unless a galaxy is seen edge-on ($i>80^{\circ}$), most of the PNe with [O III] absolute magnitudes in the top 1 mag of the PNLF will still be in the top 1 mag even after being reddened by foreground dust. Moreover, because of the symmetry of the problem, the mean extinction affecting PNe from high-mass progenitors is not significantly different from that of PNe evolved from $\sim 1M_{\odot}$ stars.

Obviously, this analysis is quite simplified, as it assumes the vertical distributions of stars and dust in our program galaxies are the same as those measured for the solar neighborhood. This will not always be true: for instance, the projected dust distribution in M31’s disk is quite complex, and the galaxy is known to have a ring of enhanced extinction ($A_{V}\sim 3$ mag) $\sim 9$ kpc from its center (Dalcanton et al., 2015; Blaña Díaz et al., 2018). However, all but four of the M31 PNe listed in Table 1 are well beyond this feature (the median projected galactocentric radius of the sample is over 17 kpc); at such large distances, M31’s observed extinction is relatively low. Additionally, Bhattacharya et al. (2019) have shown that the observed Balmer decrements of 413 of M31’s PNe do not correlate with the position of the aforementioned dust ring or any other feature in the galaxy’s line-of-sight extinction map. Again, this supports the argument that internal extinction will not bias our results.

Our analysis implies that as long as the total extinction in a sight line through a galaxy at the position of its PNe does not become greater than $A_{5007}\sim 1$ mag, the effect of foreground dust should be primarily random in nature; systematic shifts with stellar population should be below the level at which a bias is introduced into our analysis. With an inclination of $i\sim 35^{\circ}$, the LMC’s expected internal extinction of $A_{5007}\lesssim 0.2$ mag falls well below this threshold ($c_{{\rm H}\beta}\sim 0.08$). Due to its more extreme inclination ($i\sim 77^{\circ}$), the effect of internal extinction on M31’s PNe is more serious, but even then, we expect $A_{5007}\lesssim 0.7$ mag ($c_{{\rm H}\beta}\lesssim 0.28$) for our outer disk sample of objects.

Since the floor in the distribution of PN extinctions appears close to that found from Milky Way extinction alone, we do not correct our measurements for unrelated foreground extinction within the program galaxies. Instead, we simply treat the component as a source of random scatter in the data.

The errors in $c_{{\rm H}\beta}$ come from comparing the PNe’s observed Balmer decrements with expectations from low-density plasmas (e.g., Pengelly, 1964; Brocklehurst, 1971; Storey & Hummer, 1995). For the Ueta & Otsuka (2022) PNe, these errors are quoted explicitly in their Table 1; for the Balick et al. (2013) and Fang et al. (2018) objects, $\sigma(c_{{\rm H}\beta})$ can be found by propagating the quoted uncertainties on each emission line strength (and assuming the error on the intrinsic value calculated from the modeling is negligible). For the remaining objects, no information on the $c_{{\rm H}\beta}$ errors is given. For these PNe, we assume an uncertainty of 0.036 dex, which is the median of the Balick et al. (2013), Fang et al. (2018), and Ueta & Otsuka (2022) errors. A variance test shows that if we treat the highly discrepant PN Merrett 2860 as an outlier, then the scatter in $c_{{\rm H}\beta}$ for the remaining three PNe with multiple measurements is consistent with this estimate.

Next, we need to consider the random error on our core masses. As described above, these come from interpolating the star’s photoionization model-based $\log T_{\rm eff}$ and $\log L$ values in the grid of MB evolutionary tracks. None of the papers give errors for the central star parameters, but based on our experience with photoionization models and the various trade-offs between the fitted variables, we believe that $\sigma(M_{c})\sim 0.008$ is a reasonable estimate for the random component of the core mass error.²²2 Since our targets for analysis are among the brightest PNe in the galaxy, they are most likely optically thick to ionizing radiation. In these types of objects, Kwitter et al. (2012) derive central star luminosities (and hence masses) primarily via the extinction-corrected fluxes of the brightest Balmer lines, which have a quoted error of $10-15\%$. When we increment the reported stellar luminosities by this amount, we find that the mean mass of a central star increases by $\sim 0.008M_{\odot}$. Note that this error does not include the systematic component associated with our choice of post-AGB evolutionary tracks, nor does it reflect variations in derived core masses based on significantly discrepant emission-line observations.

Finally, we must consider the likely presence of a correlated error between extinction and core mass. This systematic comes about because any overestimate of extinction may cause the photoionization models to overestimate the total luminosity of line emission from the nebula and, therefore, (through energy balance) the total luminosity of the exciting star. This error can then propagate directly into the core mass values.

To estimate the amplitude of this correlated error, we combine our estimate of $\sigma(c_{{\rm H}\beta})$ with an $R_{V}=3.1$ Cardelli et al. (1989) extinction law and calculate the effect that a given error in H$\beta$ extinction has on the total amount of [O III] luminosity being emitted. Because all the objects in our sample are dominated by [O III] emission, this effect is always very close to $0.96\,c_{{\rm H}\beta}$.

At the same time, we estimate the (very small) effect an error in $c_{{\rm H}\beta}$ has on $\log T_{\rm eff}$ using the He I $\lambda 4686$/H$\beta$ ratio and the cross-over method (Kaler & Jacoby, 1989). Although this technique only works for optically thick nebulae with hot ($T>80,300$ K) central stars, we expect that most of the PNe in the top $\sim 1$ mag of the [O III] PNLF should at least come close to satisfying these conditions. To wit: if Lyman continuum photons are leaking out into interstellar space, then the efficiency of [O III] production would necessarily decrease, leading to fainter [O III] magnitudes. Similarly, PNe with cooler central stars will have a larger fraction of oxygen in the singly ionized state, again suppressing the efficiency of [O III] emission. But most importantly, even if a PN is not optically thick and/or has a relatively cool central star, the core’s location on the HR diagram should guarantee that it is in its final, near-constant-luminosity trek across the HR diagram (see Fig. 4); for these objects, any error in effective temperature has little, if any, effect on the measurement of core mass.

Finally, we use the error vectors calculated above to find new values for $\log L$ and $\log T_{\rm eff}$ and repeat Step 6 to find the expected change in core mass due to our hypothesized error in extinction. The random and systematic components are then combined to create a covariance matrix for each PNe. These matrices are represented by the error vectors shown in Figure 5.

From the figure, it is clear that our derived $c_{{\rm H}\beta}$-core mass relation contains a significant amount of scatter; this comes both from measurement error and lack of knowledge about foreground extinction affecting each PN. Nevertheless, the datasets for both galaxies follow the overall trend found by Ciardullo & Jacoby (1999) and Kwitter et al. (2012): there is a correlation between the circumnebular extinction of [O III]-bright PNe and their core mass. Interestingly, although the LMC PNe have half the metallicity of their M31 counterparts (which comes principally from the outer disk), the extinction-core mass relation for the two systems are similar. If we apply the orthogonal-regression (OR) algorithm of York et al. (2004), which fits heteroskedastic data in the presence of correlated errors, we obtain slopes of $8.9\pm 2.6$ for the LMC, $5.3\pm 1.7$ for M31, and $7.0\pm 1.8$ for the combined dataset.

We note that simpler algorithms, such as ordinary least squares (OLS), yield shallower slopes for the $c_{{\rm H}\beta}$-core mass relation, with values generally between $\sim 3$ and $\sim 5$. This is to be expected: while OLS assumes the observed slope is entirely due to the relationship between the two variables, the York et al. (2004) OR algorithm treats the slope as being due to two components, one intrinsic to the data and the other a product of the correlated errors. But regardless of the type of least squares algorithm used, the key result of the analysis is the same: the correlation between $c_{{\rm H}\beta}$-core mass is real. Even if our best estimates for the systematic component of the errors were increased by a factor of $\sim 3$, there would be no qualitative difference in the result.³³3Our LMC figure has far less scatter than the equivalent version in Ciardullo & Jacoby (1999). Unlike the previous work, the present sample is restricted to PNe within the top one mag of the [O III] luminosity function. This more appropriate selection criterion results in a more uniform PN population.

Table 2 summarizes our best-fit relations. Given the heterogeneous nature of the data and our lack of knowledge about true errors on the core masses, the formal uncertainties for the fitted slopes and intercepts may be underestimates. In fact, depending on the assumptions used for our fitting procedure, our analysis admits a wide range of slopes, from $b\sim 4$ to $b\sim 9$. Similarly, because the intercept values are also uncertain, we include an extra column in the table labeled “Intercept*”, which given the best-fit slope, is the intercept needed to ensure that the peak value of $M_{5007}$ matches the observed value of $M^{*}$. In all cases, the offset between Intercept* and the formal best-fit intercept is quite small, i.e., just a few percent of the uncertainty.

Of course, since the composition and distribution of the AGB dust ejecta is expected to vary with stellar mass, due to the changing physics of dredge-up and time since the AGB (e.g., Karakas & Lugaro, 2016; Dell’Agli et al., 2023), we might expect the actual relation to have features (e.g., non-linearities and discontinuities) that are hidden within the uncertainties of the measurements and the limited PN samples. Nevertheless, for purposes of this study, a simple linear relation will do, and the large positive slope suggests an answer to the long-standing puzzle of the constancy of the PNLF cutoff.

Because our explanation for the invariance of the PNLF cutoff hinges on the systematic behavior of circumnebular extinction, it would be useful to have an independent way of verifying that the results of Figure 5 are not simply due to correlated errors. The systematics of PN chemical abundances provides such a test.

Without circumnebular extinction, the brightest PNe within a star-forming galaxy should be those that evolve from the highest mass ($3\lesssim M/M_{\odot}\lesssim 8$) PN progenitors. Such stars have a unique signature: unlike their lower-mass counterparts, high-mass stars undergo a third dredge-up and hot-bottom burning while on the AGB (e.g., Herwig, 2005). As a result, PNe which evolve from these objects will be enhanced in nitrogen and helium. In other words, the most luminous PNe observed in a galaxy should be Type I planetaries.

As defined by Peimbert & Torres-Peimbert (1983)⁴⁴4See also Torres-Peimbert & Peimbert (1997) for a discussion of alternative definitions of Type I PNe. Under some of these definitions, several of the M31 PNe marginally satisfy one of the Type I criteria, i.e., are slightly overabundant in either helium or nitrogen., Type I PNe are nebulae with $\log{\rm N/O}\geq-0.3$ and He/H $>0.125$. None of the 28 M31 PNe analyzed by Kwitter et al. (2012), Balick et al. (2013), Fang et al. (2018), and Ueta & Otsuka (2022) satisfy these criteria. Similarly, only one of the 14 LMC PNe included in our analysis has a Type I classification (SMP 47). Thus, at best, only $\sim 3\%$ of PNe in the top magnitude of the [O III] luminosity function are formed from high-mass stars.

We can compare this number to the results of Jacoby & Ciardullo (1999), who analyzed the abundances and central star properties of PNe in the bulge and disk of M31. These authors found that over a quarter of their sample (4 out of 15) satisfy the Type I criteria, but most interestingly, all the Type I objects are two to three magnitudes fainter than $M^{*}$. This result supports the notion that the PNe which populate the bright end of the luminosity function do not evolve from high-mass progenitors.

We note that circumnebular extinction may not solely be responsible for the lack of Type I planetaries in the top magnitude of the PNLF; because of their shorter post-AGB lifetimes, the chance of observing a Type I PN is smaller than that of finding PNe from lower-mass stars. Unfortunately, at present, the relative importance of the two mechanisms is difficult to quantify.

We can calculate the relative rates of PN formation against their progenitor masses from the application of the fuel-consumption theorem (Renzini & Buzzoni, 1986). But any calculation for the number of PNe expected to be seen in a survey of the top $\sim 1$ mag of the PNLF would not only require detailed modeling of the time evolution of the PNe’s central stars, their nebulae, but also the formation/destruction of their dust and its spatial distribution. Such calculations are beyond the scope of this paper.

Instead, we compare the fraction of Type I PNe in the Galaxy to the fraction seen in our sample of LMC/M31 objects. There are $\sim 3500$ confirmed (and unconfirmed) PNe cataloged in the HASH database of Milky Way planetary nebulae (Parker et al., 2016), and roughly 2% are classified as Type I (Phillips, 2004). This is similar fraction found for the objects listed in Table 1 (one out of 42).

So, although we cannot quantify the relative importance of either of the two mechanisms, we see no clear argument that the short visibility time for Type I PNe is dominating. Furthermore, Figure 6 demonstrates that rapid evolution does not dramatically suppress the presence of stars with progenitor masses up to $\sim 3.3M_{\odot}$ in the top mag of the luminosity function.

Gesicki et al. (2018) showed that one can derive peak [O III] luminosities for PNe that meet or exceed the observed value of $M_{5007}\sim-4.5$ by simulating PNLFs for stellar populations younger than $\sim 3$ Gyr (i.e., progenitor masses greater than $\sim 1.5M_{\odot}$). Their analysis tracks stellar evolution from the main sequence through the post-AGB phase, follows the photoionization history of the post-AGB stars’ nebulae, and includes a prescription for extinction that involves the progenitor AGB star’s superwind mass-loss rate and velocity.

This paper does not attempt such comprehensive modeling; instead, it focuses narrowly on assessing the maximum [O III] luminosity achievable for post-AGB stars of various masses, using knowledge of the maximum efficiency for converting stellar luminosity into [O III] $\lambda 5007$ emission and an empirical prescription for the behavior of circumnebular extinction with core mass. With this approach, the uncertainties associated with nebular physics and galactic star formation histories are irrelevant. While our analysis prevents us from deriving an actual luminosity function, it does generate guidelines for future analyses.

Ultimately, our results are not so different from those of Gesicki et al. (2018). However, our analysis shows that not only is dust necessary for modeling the PNLF, but it is perhaps the driving force behind the invariance of the PNLF’s bright-end cutoff. Moreover, this conclusion follows directly from observations rather than any attempted mating of stellar evolution and nebular modeling. As long as a stellar system has stars with main sequence masses of $1.25M_{\odot}$ or greater (younger than $\sim 5$ Gyrs), dust ensures that the PNLF’s high-luminosity cutoff will yield a distance that is accurate to $<2.5\%$, provided the observed PN sample is sufficiently large to populate the brightest $\sim 1$ mag of the PNLF with $\gtrsim 30$ objects.

As discussed in Section 4, PNe which evolve from stars with initial masses $M\gtrsim 3M_{\odot}$ generally do not populate the bright-end of the PNLF. This is illustrated by the solid curves of Figure 6, all of which display a turn-over at high masses. While interesting in theory, this feature has no practical effect on PNLF distances. Stars with initial masses of $M\gtrsim 2M_{\odot}$ have relatively short lifetimes (less than $\sim 1$ Gyr), and there are no large galaxies in the local universe that are composed solely of stars this young. Even in the most vigorous star-forming galaxies, the PNLF cutoff will be defined by the objects that have evolved from the underlying older populations. The existence of this fading for young populations is, therefore, a moot issue for distance determination.

A more serious concern is the unobserved fading of the PNLF in old stellar populations. Elliptical galaxies are dominated by stars with ages $\gtrsim 6$ Gyr (e.g., Trager et al., 2000; Zibetti et al., 2020; Parikh et al., 2024); thus, according to Figure 6, their PNLF cutoffs should be $\gtrsim 0.5$ mag fainter than those seen in spiral galaxies.⁵⁵5 While a moderate fraction of field ellipticals show signs of recent star formation, the phenomenon is rare in groups and clusters, where most PNLF measurements have been made (Paspaliaris et al., 2023). In any case, since the amount of luminosity in this “young” component is usually small, the luminosity-specific stellar evolutionary flux (Renzini & Buzzoni, 1986; Buzzoni et al., 2006) predicts that it would not be sufficient to populate the bright-end of the PNLF. Deviations of this magnitude would produce PNLF distances that are $\sim 25\%$ greater than those found by other methods. This is not the case: if there is any systematic bias to the PNLF, it is towards too small a distance for old populations, and its amplitude is less than $\sim 10\%$ (Ferrarese et al., 2000; Feldmeier et al., 2007).

This mismatch between the predicted and observed PN luminosities for Pop II systems has been discussed several times in the literature (e.g., Jacoby et al., 1989; Marigo et al., 2004; Ciardullo et al., 2005; Gesicki et al., 2018; Yao & Quataert, 2023) along with various possible solutions: mass enhancement through binary interactions and stellar mergers (Ciardullo et al., 2005), the presence of a small population of intermediate-age stars whose signature is masked by red giants (Jacoby, 1997; Ciardullo et al., 2005), an incorrect reference PNLF (Hartke et al., 2017), contamination by Pop I emission-line sources such as Wolf-Rayet nebulae (Jacoby et al., 2024), the presence of interloping $z\sim 3.1$ Ly$\alpha$-emitting galaxies in the PN sample (Kudritzki et al., 2000; Jacoby et al., 2024), or even contamination by post-AGB objects whose nebulae are powered by accretion onto companion white dwarfs (Soker, 2006; Souropanis et al., 2023). To these, we add three additional possibilities:

As illustrated in Figure 5, the extinction-core mass relation exhibits a large amount of scatter. Some of this scatter is likely due to internal extinction in the host galaxies; the slightly larger scatter in the M31 dataset lends some support for this interpretation, although the more heterogeneous nature of the data set may also be responsible.

Another factor is undoubtedly measurement error, as the process of using a limited number of optical (and in some cases, near-IR) emission lines to fix a central star’s position on the HR diagram is quite complex. There are a large number of trade-offs and degeneracies associated with the process, and the result can be non-unique solutions for the stellar luminosity and temperature.

It is also likely that much of the dispersion is intrinsic. The distribution of dust in bipolar PNe is notoriously asymmetric, and even in roughly spherical nebulae, walls, ridges, and filaments of extinction are easily seen. (For excellent examples, see the in-depth study of NGC 7027 by Moraga Baez et al. (2023) and the set of PN extinction maps created by Pignata et al. (2024).) Thus, the value of $c_{{\rm H}\beta}$ for distant unresolved PNe is likely sight-line dependent. This effect alone would introduce scatter into Figure 5.

Finally, the effects of extinction on the [O III] luminosity of planetary nebulae are undoubtedly time-dependent, as dust may get destroyed or dissipate. In this study, we attempted to create as homogeneous a PN sample as possible by restricting our analysis to objects within 1 mag of the PNLF cutoff. Yet, as Fig. 3 shows, there is a suggestion of evolution in our data, as the objects in the top 0.6 mag of the PNLF are more efficient at producing [O III] $\lambda 5007$ than their fainter counterparts. This would also introduce intrinsic scatter into Figure 5.

If there is scatter in the extinction-core mass relation, then the shape of the brightest $\sim 0.5$ mag of the PNLF may largely be defined by stochastic effects associated with the dust distribution, such as viewing angle. A much larger sample of PNe will be needed to test this hypothesis.