Revised Simulations of the Planetary Nebulae Luminosity Function

If we can measure the nebular emission line fluxes of [O iii] 5007 for many planetary nebulae (PNs) in a galaxy and transform the fluxes into apparent magnitudes m(5007) (Jacoby magnitudes, defined by Jacoby 1989), then we can build the planetary nebulae luminosity function (PNLF). It indicates how many PNs have [O iii] 5007 at each apparent magnitude m(5007). Empirically, it was discovered that the brightest PNs in a galaxy have similar absolute magnitudes M(5007). This led to the suggestion of a universal bright end of the PNLF that could be used as a standard candle (Jacoby 1989). The PNLF has been used to determine extragalactic distances for almost 30 yr now (beginning with Ciardullo et al. 1989; Jacoby 1989) and has proven to be a reliable distance indicator (see e.g., Ciardullo 2012).

The PNLF not only gives insight into the distance of a galaxy, but also into properties of the stellar population. Consider, for example, the PN populations of elliptical galaxies and of M31's bulge. The bright ends of their observed PNLFs require the existence of very luminous central stars, approaching 7000 L_⊙ (see e.g., Jacoby & Ciardullo 1999; Méndez et al. 2005, 2008b). Because of the luminosity–core mass relation for post-AGB stars, this requires very massive central stars (about $0.63\,{M}_{\odot }$), in strong disagreement with the expected maximum final mass of at most $0.55\,{M}_{\odot }$, which would correspond to a turn-off initial mass of about $1\,{M}_{\odot }$ in these rather old stellar populations (Zhao et al. 2012; Cummings 2017; El-Badry et al. 2018). The problem has recently become more severe after individual values of internal dust extinction have become available for several of the brightest PNs in M31's bulge (Davis et al. 2018), suggesting some dereddened central star luminosities in excess of 10,000 L_⊙.

This alarming discrepancy has been somewhat reduced by the introduction of new post-AGB evolutionary tracks by Miller Bertolami (2016). A modified luminosity–core mass relation has decreased the required central star mass from 0.63 to 0.58 M_⊙ (Méndez 2017), in better agreement with theoretical expectations. The new tracks also show a much faster post-AGB evolution than previously calculated. This means that central stars with much lower masses than previously expected can produce visible PNs.

If we further assume that the brightest PNs are predominantly optically thick (opaque) to H-ionizing radiation from their central stars (Gesicki et al. 2018), it becomes possible to produce a bright end of the PNLF, which stays almost invariant for a variety of stellar population ages and star formation histories. This helps to understand the high quality of the PNLF as a standard candle.

However, there are good reasons to doubt that most PNs are completely optically thick. Looking around in our Milky Way, we see PNs with predominantly bipolar symmetry, which means that they will probably start leaking H-ionizing photons very soon, in directions where the nebular density is lower.

In view of all these complications, PNLF simulations created using Monte Carlo techniques can be used to take a closer look at how sample size, population age, and partial absorption of ionizing photons affect the properties of the PNLF (Méndez et al. 1993; Méndez & Soffner 1997; Méndez et al. 2008b). The recent advances in post-AGB theory have induced us to revisit PNLF simulations and verify how far we are from really understanding the PNLF.

In the present work we describe our revisions of the procedure outlined by Méndez & Soffner (1997) to generate PNLFs, using the new evolutionary tracks for low-mass stars (Miller Bertolami 2016) and a central star mass distribution derived from a modern DA white dwarf mass distribution (Kepler et al. 2017). The newly simulated PNLFs obtained in this way are compared with the observed PNLFs in a few selected galaxies to find out what the consequences are of the new tracks, regarding the shape and bright end of the PNLF, the unresolved tension between PNLF distances and surface brightness fluctuation (SBF) distances (Ciardullo 2012; Méndez 2017), and the PN formation rate.

First of all, the reader might ask why our PNLF generation procedure follows the simple one described in Méndez & Soffner (1997), instead of the more elaborate procedure described in Méndez et al. (2008b), which incorporated information derived from nebular models of Schönberner et al. (2007) as well as from observational work by Ciardullo et al. (2002). Let us justify this choice. Méndez et al. (2008b) included in its Figure 6 a comparison with Figure 2 of Ciardullo et al. (2002), showing good agreement. Furthermore, Figure 7 of Méndez et al. (2008b) compared the PNLF bright end generated using the models of Schönberner et al. (2007) with the PNLF bright end produced with Mendez & Soffner simulations, again showing good agreement. From these two comparisons we conclude that the method introduced in Méndez & Soffner (1997) is compatible with the empirical information later provided by Ciardullo et al. (2002). For that reason we avoided the approach used in the 2008 paper; it would have required a complete recalculation of the nebular model grid, using the new post-AGB evolutionary tracks of Miller Bertolami (2016), characterized by different post-AGB central star masses and evolutionary speeds. We decided that for the purpose of studying the impact of the new evolutionary models on the PNLF, it was simpler to use the 1997 simulation method.

It should be clear that a simulation based partly on random numbers cannot be perfect; the random numbers are indeed an attempt to compensate for the lack of complete information about nebular properties. For example, if we consider 5007 intensities, a full specification would require knowing the PN metallicity distribution, about which we do not know enough. So our approach is to extract empirical information from observed cases. If we compare with theoretical nebular models like Dopita et al. (1992), which predict that "at high metallicities up to 15% of the luminosity of the central star is reradiated in the 5007 line alone," then we might find that we are generating a few PNs with a higher ratio L(5007)/Lstar, may be about 18%. However, any distorting effect produced by the generation of a few too high 5007/Lstar ratios through random numbers can be compensated by a few nebulae leaking too much ionizing radiation (which is expected according to Section 2.3 in Méndez et al. 2008b), without affecting the general picture in any significant way. In summary, our simulations, although not perfect, are good enough to study the effect of introducing the new post-AGB evolutionary tracks, which is our purpose here.

Because the procedure we have selected closely follows the one described in Méndez & Soffner (1997), the reader may refer to that paper for more details. For brevity, here we will focus on the changes we have introduced. The evolutionary tracks of H-burning post-AGB stars are now taken from Miller Bertolami (2016). More specifically, we have selected five tracks corresponding to an overall metal mass fraction ${Z}_{0}=0.01$ from his Table 2, with initial masses 1.0, 1.25, 1.5, 2.5, and 3.0 M_⊙. These tracks are based on improved models of AGB stars and evolutionary calculations. The main differences with previous work are that the new post-AGB stellar luminosities are higher and post-AGB timescales are shorter than what earlier models suggested.

As before, we have chosen to define the post-AGB age to be equal to the time passed since the moment the central star had a temperature of ${T}_{\mathrm{eff}}={\rm{25,000}}\,{\rm{K}}$. We will call the post-AGB ages just "ages" in what follows. Because the central star mass changes continuously with time (see Equation (6) in Miller Bertolami 2016), each track is assigned the average mass between the masses at ages 0 and 30,000 yr. This hardly affects the outcome of the interpolated tracks because the change of mass along each track is much smaller than the central star mass itself, by an order of 10⁻⁴–10⁻⁵: mass-loss rates for central stars inferred from observed spectra using the theory of radiatively driven stellar winds are of the order of 10⁻⁷–10⁻⁹ M_⊙ yr⁻¹ (Kudritzki et al. 2006).

We generated a look-up table similar to the one described in Méndez & Soffner (1997) to determine $\mathrm{log}{T}_{\mathrm{eff}}$ and $\mathrm{log}L$ as functions of 3000 ages between 0 and 30,000 yr, and 260 masses between 0.530 and 0.706 solar masses. For the interpolations, the tracks were split up into three sections.

For the nearly horizontal heating tracks we used 40 temperatures between 25,000 and 81,000 K, at which we plotted $\mathrm{log}\,\mathrm{age}$ and $\mathrm{log}L$ as functions of mass. To interpolate between the values given by the tracks, we fitted curves to these plots and calculated the age and luminosity for the 260 masses at each of the 40 temperatures.

For the white dwarf cooling tracks we used 60 luminosities between $\mathrm{log}L=1.60$ and 2.75, at which we plotted $\mathrm{log}\,\mathrm{age}$ and $\mathrm{log}{T}_{\mathrm{eff}}$ as functions of mass. We fitted curves to these plots and calculated the age and temperature for the 260 masses at each of the 60 luminosities.

For the curved section joining the heating and cooling tracks we used 90 lines radiating at different angles from a fixed point at ${T}_{\mathrm{eff}}=81,000\,{\rm{K}}$ and $\mathrm{log}L=2.75$. The values obtained at the intersections between these lines and the given tracks were used to plot $\mathrm{log}\,\mathrm{age}$, $\mathrm{log}{T}_{\mathrm{eff}}$, and $\mathrm{log}L$ as functions of mass. We then fitted curves to these plots to calculate age, temperature, and luminosity for the 260 masses at each of these 90 angles.

Having 190 values of temperature and luminosity for each of the 260 masses along their respective tracks, we could interpolate between these to obtain $\mathrm{log}{T}_{\mathrm{eff}}$ and $\mathrm{log}L$ at 3000 ages between 0 and 30,000 yr. The result was a look-up table that can be used to determine the temperature and luminosity of a star with a given age and mass by bivariate spline interpolation.

Figure 1 shows the evolutionary tracks for the five final (central star) masses used to calculate the interpolations.

Zoom In Zoom Out Reset image size

When generating a sample of PN central stars, each star is assigned a random age and mass. The ages are uniformly distributed between 0 and 30,000 yr—this is approximately the expected timescale of a PN, given typical observed expansion velocities and the largest observed PN sizes.

Because of the faster post-AGB evolution, central stars with lower masses than previously expected can produce visible PNs; but of course there must be a limit. PNs with central star masses below $0.53\,{M}_{\odot }$ should dissipate before the central stars can become hot enough to ionize the gas (see Figure 1). Therefore, we generated masses between $0.53\,{M}_{\odot }$ and a maximum final mass that we expected to be somewhere below $0.60\,{M}_{\odot }$. We approximated this range of masses as a linear distribution, with a mass of $0.58\,{M}_{\odot }$ being 2.5 times as likely as a mass of $0.53\,{M}_{\odot }$. This approximation is based on the DA white dwarf mass distribution for ${T}_{\mathrm{eff}}\geqslant 13,000\,{\rm{K}}$ corrected by the $1/{V}_{\max }$ method in Kepler et al. (2017). We compare both mass distributions in Figure 2.

Zoom In Zoom Out Reset image size

The Hertzsprung–Russell (HR) diagram of 1500 PNs generated in this manner is shown in Figure 3. The values of central star temperature and luminosity were computed using the procedure described in Section 2.

Zoom In Zoom Out Reset image size

For a PN that is completely optically thick to the Lyman continuum, it is possible to derive the ${\rm{H}}\beta$ luminosity for given values of central star temperature and luminosity under the assumption of a blackbody energy distribution (Osterbrock & Ferland 2006). See more details in Méndez et al. (1993).

The simulated central star mass distribution was based as much as possible on empirical information. Unfortunately, for obvious reasons, the only such information available is for local white dwarfs. What would be the difference with populations without any recent star formation? We believe that the most important difference would be the presence in the local white dwarf sample of more massive white dwarfs, produced by those massive stars that have been able, despite their more recent birth, to complete their evolution and become white dwarfs. But it is precisely those massive white dwarfs that we eliminate from our sample by introducing the "maximum final mass." Therefore, we do not expect any problem from this kind of age difference. On the other hand, we will discuss possible differences in the history of the star formation rate, e.g., the possible lack of lower-mass white dwarfs, in Section 9 below, in connection with the PNLF shape.

To compensate for the changes in mass distribution and evolutionary tracks, we slightly adjusted the procedure that generates the λ5007 intensity, I(5007), relative to the intensity of ${\rm{H}}\beta$. We continued to compare our distribution to the one found in the Large Magellanic Cloud (LMC; using 118 PNs) and to the one found in the Milky Way (using 983 PNs). In both cases, we used the compiled data described in Section 4 of Méndez & Soffner (1997).

Our distribution is generated by following the previous procedure, with the following differences. We now use a Gaussian distribution centered at $I(5007)=950$ on a scale of $I({\rm{H}}\beta )=100$ with $\mathrm{FWHM}=375$. We then lower the values above 1200 by up to 200 (values right above 1200 are only slightly decreased while values around 2000 are maximally decreased) to imitate the steeper drop-off on the right part of the distribution. We then cap the values at 1850. Because of the lower masses being used, we lower I(5007) by $60 \%$ for central stars with masses below $0.55\,{M}_{\odot }$ that are on heating tracks with ${T}_{\mathrm{eff}}\gt 75,000\,{\rm{K}}$ (instead of by $50 \%$ for masses below $0.57\,{M}_{\odot };$ following the discussion in Section 4 of Méndez et al. 1993). The new values were chosen to allow our generated distribution to fit the shapes of the observed ones in the LMC and the Milky Way as well as possible (see Figure 4). We did not try to reproduce the Milky Way peak at very low I(5007) because it has no effect on the bright end of the PNLF.

Zoom In Zoom Out Reset image size

For comments about the relative insensitivity of I(5007) and its effect on the PNLF as a function of metallicity, please refer to Jacoby et al. (1992).

In Méndez et al. (1992), an absorbing factor μ was introduced to describe the optical thickness of a nebula in the H Lyman continuum. Its value is equal to the fraction of ionizing stellar photons that are absorbed by the nebula, which means that μ = 1 represents a completely optically thick or opaque PN, whereas a perfectly transparent PN has μ = 0. A range of absorbing factors was used in PNLF simulations by Méndez & Soffner (1997).

Now that we were using new post-AGB tracks and a different central star mass range, we had to make sure whether this absorbing factor was still required by the simulation to reproduce observed PNLFs. We assumed for the following that all PNs are optically thick, that is, μ = 1 for all of them. Having the temperatures, luminosities, and intensities of λ5007 relative to $I({\rm{H}}\beta )$, we could generate a PNLF following the procedure described in Méndez & Soffner (1997).

For the observed data, we used a statistically complete PNLF with 119 PNs of M31's bulge (samples A and B from Ciardullo et al. 1989). We transformed apparent into absolute λ5007 magnitudes adopting a distance modulus of $m-M=24.43$ and a foreground extinction correction for λ5007 of $0.28\mathrm{mag}$, based on $E(B-V)=0.08$ (Jacoby et al. 1992).

At this point we introduced a quantity $r$ that tells us how well a generated PNLF fits the observed one. We defined $r$ to resemble the way the eye evaluates the fit in a plot. The sum of the squared deviations of $\mathrm{log}N$ seemed to be a good indicator, providing a more mathematical procedure than we used in previous PNLF papers. This fitting procedure can be described as chi-square minimization in log-space. We used a histogram for M(5007) between −6 and −2 with a bin size of 0.2 to match typical PNLFs found in the literature. Since numbers below $\mathrm{log}N=0$ are much less relevant and because of ${\mathrm{lim}}_{N\to {0}^{+}}\mathrm{log}N=-\infty$, we used $\mathrm{log}(0.5)$ instead of $\mathrm{log}N$ for values below $\mathrm{log}(0.5)$. We used this minimum value because $N\lt 0.5$ means that the probability of finding a PN in that particular bin is less than $50 \%$.

By adjusting the maximum final mass and the sample size in our simulated PNLFs, we found that a mass range of 0.530–0.570 solar masses and a sample size of 285 PNs gives the best fit to the observed PNLF of M31 (see Figure 5).

Zoom In Zoom Out Reset image size

Specifically, when calculating $r$ for the three selected maximum final masses in Figure 5, we got 0.30 for $0.560\,{M}_{\odot }$, 0.20 for $0.570\,{M}_{\odot }$, and 0.37 for $0.580\,{M}_{\odot }$. Clearly, a maximum final mass of 0.570 solar masses gives us the best fit.

While the bright end can be fitted properly, the shape of the PNLF beyond $M(5007)=-3.5$ has a valley that does not match the observed PNLF particularly well. And, most importantly, the sample size of 285 PNs is too small if we compare it to 970, the total expected number of PNs in the region of M31's bulge sampled by Ciardullo et al. (1989); see the discussion in their Section V.

If we now increase the sample size by a factor of three, while keeping M31's PNLF at the right distance, the simulated PNLF fails to fit. We conclude, as already stated by Méndez & Soffner (1997), that it is necessary to allow for PNs that leak ionizing photons from the central star. The empirical evidence for this is taken from Table 4 in Méndez et al. (1992), where we find a variety of μ values, from 1 to less than 0.1, and a clear trend of decreasing μ with increasing central star temperature. This UV photon leaking is further discussed by Méndez et al. (2008b).

We have decided to avoid using any theoretical guidance about nebular optical depth taken from modern nebular modeling efforts (e.g., Schönberner et al. 2010). Such radiation-hydrodynamics models are 1D, in other words spherically symmetric. The observed predominance of non-spherically symmetric PNs introduces uncertainty in our ability to predict what fraction of H-ionizing radiation is being lost through the nebular poles. Therefore, as in previous work, we generated a distribution of μ values using random numbers. Further justification for this approach can be found in Section 2.4 of Méndez et al. (2008b). We implemented similar boundary conditions as Méndez et al. (1993) and Méndez & Soffner (1997); PNs at the beginning of the heating tracks tend to be more opaque and will become increasingly transparent at higher temperatures. For the slowly evolving low-mass central stars (see Figure 1), we expect the nebula to dissipate before reaching higher temperatures, so we assign low random absorbing factor values to them for higher ages. And finally, we use a similar method to Méndez & Soffner (1997) for generating decreasing μ-values for increasing ages on the cooling tracks, such that μ = 0 for an age of 30,000 yr. The resulting absorbing factor distribution in the HR diagram is shown in Figure 6.

Zoom In Zoom Out Reset image size

Using the absorbing factor distribution from the previous section, we could again generate PNLFs for different maximum final masses. The best fit to M31's observed PNLF was found at a sample size of 700 and a maximum final mass of between 0.580 and 0.590 solar masses (see Figure 7). Again, we used the fitting parameter $r$ defined in Section 5; we obtained 0.39 for $0.570\,{M}_{\odot }$, 0.18 for $0.580\,{M}_{\odot }$, and 0.14 for $0.590\,{M}_{\odot }$. This clearly gave us a better fit than when using completely opaque PNs ($r=0.20$ for $0.570\,{M}_{\odot }$; see Figure 5). We omitted masses beyond $0.590\,{M}_{\odot }$ because the slopes of their PNLFs' bright ends do not match the observed one. Because of the distribution of absorbing factors, we required a much larger sample size. This is more reasonable for M31's bulge and confirms that it is clearly preferable to allow for many optically thin PNs when simulating PNLFs.

Zoom In Zoom Out Reset image size

At this point it is relevant to mention a recent paper by Davis et al. (2018) reporting PNs in M31 with high internal extinction and very luminous and massive central stars. Ideally, we would like to know all the individual extinction values and plot a fully dereddened PNLF. As we only know a few individual extinction values, that is currently not possible. In fact, Davis et al. (2018) show that correcting just a few PNs destroys the typical PNLF shape (see their Figure 14). The only practical option is to assume an average value for the extinction. If we adopt a higher average extinction, the whole PNLF shifts toward higher luminosities, and therefore the maximum final mass needs to increase. We will come back to this point later on.

Because of the small number of PNs surveyed by Ciardullo et al. (1989) in M31's bulge, the uncertainty of the maximum final mass is high. For a better estimate of this value, we needed to consider an old stellar population with a larger sample size. For this purpose we used NGC 4697, a flattened elliptical galaxy in the Virgo southern extension. For the observed data, we used the method described in Méndez et al. (2001) to get a statistically complete bright end of the PNLF consisting of 328 PNs. The data were taken from the corresponding catalog (Méndez et al. 2008a).

On the assumption that the PNLFs of M31's bulge and NGC 4697 are identical, we adopted the distance modulus $m-M=30.1$ determined by the PNLF method (Méndez et al. 2001). The foreground extinction correction for λ5007 was $0.105\mathrm{mag}$, based on $E(B-V)=0.03$ from Tonry et al. (2001).

Again using the absorbing factor distribution from the previous section, we fitted our simulated PNLF to the observed one. We found the following best-fitting values: a mass range of 0.530–0.580 solar masses and a sample size of 3000 PNs (see Figure 8). Calculating the fit factor $r$ for the three maximum final masses gave us 0.16 for $0.570\,{M}_{\odot }$, 0.09 for $0.580\,{M}_{\odot }$, and 0.34 for $0.590\,{M}_{\odot }$. This time, we got a very clear best fit. The reason for this is the much larger sample size in NGC 4697 compared to that in M31. We conclude that our new PNLF simulations work very well with a maximum final mass of $0.58\,{M}_{\odot }$, confirming a preliminary result by Méndez (2017).

Zoom In Zoom Out Reset image size

To analyze the properties of PNs at the bright end of the simulated PNLF with a maximum final mass of 0.58, we took a closer look at the PNs with absolute λ5007 magnitudes brighter than $-4.2\,$ mag. In our simulations we found that the distribution of absorbing factors for these PNs has extreme values of 1.0 and 0.56, a mean value of 0.89, and a standard deviation of 0.10.

We would like to emphasize that our PNLF simulations do not include internal dust extinction. The empirical existence of a "maximum final mass" does not mean that more luminous and massive central stars cannot exist. We need to add the condition that if more luminous central stars do exist (as in the bulge of M31), they are always affected by internal dust extinction. This is not hard to imagine, because we can expect the more massive central stars to eject more material and to evolve more quickly, leading necessarily to a denser distribution of circumstellar material. The consequence would be that the PNs originating from the more massive central stars would pile up at or near the bright end of the PNLF, mixed with those PNs produced by central stars with maximum final mass, no extra dust extinction, and absorbing factors close to 1. Of course this interpretation will need to be tested with individual spectroscopic studies of complete samples of extragalactic PNs defining the bright end of the PNLF.

The use of the PNLF for distance determinations is based on the assumption that the PNLF's bright end is universal. Any systematic difference between PNLF distances and those determined by other methods is a challenge to that assumption. The SBF method of distance determination (Tonry & Schneider 1988; Tonry et al. 2001; Blakeslee et al. 2009) has shown a tendency to produce larger distance moduli than the PNLF method, by about 0.3–0.4 mag (Teodorescu et al. 2010). We refer the reader to Ciardullo (2012) for a discussion of the most likely cause for this systematic effect: zero-point errors arising from limited knowledge of internal extinction in the calibrator galaxies, to which the PNLF and SBF methods react in opposite ways. In what follows we would like to illustrate the consequences of assuming that the SBF distances of elliptical galaxies are correct.

An increased distance forces the PNLF to become more luminous. What maximum final mass would be required for our simulated PNLF to fit an observed PNLF if shifted to the SBF distance? Consider, for example, NGC 4697. Blakeslee et al. (2009) used the SBF method to determine a distance modulus $m-M=30.491\pm 0.065$. We adopted an SBF distance as close to the PNLF distance as allowed by the SBF error bars: $m-M=30.4$. As in Section 6, we used an extinction correction of $0.105\,$ mag. With those numbers we transformed apparent into absolute λ5007 magnitudes. In order to get a fit with our simulated PNLF we needed to increase the maximum final mass to 0.61 M_⊙. The sample size also increased to 5200 (see Figure 9). The fit parameter was $r=0.34$.

Zoom In Zoom Out Reset image size

As another example we took M60, an elliptical galaxy in the Virgo Cluster. For the observed data, we used the method described by Teodorescu et al. (2011) to get a statistically complete sample of 218 PNs. We adopted their PNLF distance modulus of $m-M=30.7$ and extinction factor of 0.09 in our calculations. After fitting the simulated PNLF to the observed one, we got a sample size of 2950 and a maximum final mass of $0.58\,{M}_{\odot }$ (see Figure 10). The fit parameter was $r=0.26$.

Zoom In Zoom Out Reset image size

Next, we adopted an SBF distance as close to the PNLF distance as allowed by the SBF error bars from Blakeslee et al. (2009): $m-M=31.0$ (they found $m-M=31.082\pm 0.079$), and we recalculated the absolute λ5007 magnitudes. Again, to force a fit with our simulated PNLF, we needed a maximum final mass of $0.61\,{M}_{\odot }$ and a sample size of 5350 (see Figure 11). The fit parameter was $r=0.28$.

Zoom In Zoom Out Reset image size

A summary of our tests is given in Table 1. In both cases (NGC 4697 and M60) the PNLF distance leads to a maximum final mass of $0.58\,{M}_{\odot }$, whereas the SBF distance leads to $0.61\,{M}_{\odot }$. The problem with accepting the SBF distances is, first, that we do not expect such massive central stars to originate from an old stellar population if we assume single post-AGB stellar evolution; and second, that the PNs defining the bright end of the PNLF in galaxies like NGC 4697 and M60 become more luminous than any PNs ever discovered in our Local Group. If SBF distances are confirmed, we will need to find explanations for these anomalies.

Now we see that Davis et al. (2018) may have provided such an explanation. Imagine that M31's bulge (perhaps every spiral bulge) has PNs that suffer from higher average values of internal dust extinction than what we find in elliptical galaxies like NGC 4697 and M60, perhaps as a consequence of a higher average metallicity. Assuming a higher average extinction for M31 than previously adopted, we need to shift its PNLF toward higher luminosities. In order to recover the fit, as we are using M31 as a calibrating galaxy, we need to shift the elliptical galaxy's PNLF by increasing its distance modulus to values closer to those obtained with the SBF method. Note that the difference between the average circumstellar extinctions for M31 and NGC 4697 reported by Davis et al. (2018) is 0.37 mags, similar to the systematic difference between SBF and PNLF distance moduli for distant elliptical galaxies.

There would be a price to pay: all the galaxies then would have very bright PNs, requiring very luminous and massive central stars. We would agree with Davis et al. (2018) that this situation would require rethinking the whole interpretation of the PNLF, probably requiring a significant contribution from other factors than only single post-AGB star evolution.

The impact would not be limited to PN research; our current understanding of stellar populations comes from population synthesis models, which assume in particular that we understand AGB stars. A careful revision of every result obtained using the old AGB and post-AGB stellar evolution calculations would then become a very prudent precaution.

Probably the best way to test SBF and PNLF distances will be to use Tip of the Red Giant Branch (TRGB) distances. The method is described by Makarov et al. (2006); it has been applied to many galaxies up to distances of about 10 Mpc. Ideally, we need at least ten TRGB distances in the Virgo and Fornax clusters to yield a statistically convincing result. This will be possible with the James Webb Space Telescope, which is expected to allow TRGB distance determinations up to at least 30 Mpc. If the SBF distances are confirmed to be closer to the truth than PNLF distances, then the explanation in terms of different amounts of average internal extinction becomes very likely.

Further investigation of this idea (different average internal dust extinction) would probably be best oriented toward obtaining additional deep spectrograms of PNs in elliptical galaxies like M60, with the purpose of measuring the individual Balmer decrements. Metallicity determinations from integrated light analysis indicate that NGC 4697 may be more metal-poor than the bulge of M31, while M60 (NGC 4649) may be not (see Trager et al. 2000). But this may not be decisive, because of the presence of a metallicity gradient (with detected PNs located predominantly in low surface brightness regions) and because of unknown complications in the details of dust formation.

A fit to the observed PNLF with our simulated PNLF provides a simultaneous determination of distance modulus and sample size (the total number of PNs in the surveyed area of the galaxy). What is the effect of the new Miller Bertolami tracks on PN formation rates? Having obtained new sample sizes, we can recalculate the specific PN formation rate $\dot{\xi }$ with

$$\begin{eqnarray}&&{n}_{\mathrm{PN}}=\dot{\xi }{L}_{T}{t}_{\mathrm{PN}},\end{eqnarray} \tag{ 1 }$$

where ${n}_{\mathrm{PN}}$ is the sample size, ${L}_{T}$ is the total luminosity of the sampled population, and ${t}_{\mathrm{PN}}={\rm{30,000}}\,\mathrm{yr}$ is the lifetime of a PN. We take ${L}_{T}$ from Méndez et al. (2001) for NGC 4697: ${L}_{T}=1.95\times {10}^{10}\,{L}_{\odot }$. With the new sample size ${n}_{\mathrm{PN}}=3000$ from Section 6 we get $\dot{\xi }=5\times {10}^{-12}\,{\mathrm{yr}}^{-1}\,{{L}_{\odot }}^{-1}$. This is slightly lower than the earlier result, $\dot{\xi }=6\times {10}^{-12}\,{\mathrm{yr}}^{-1}\,{{L}_{\odot }}^{-1}$ (Méndez et al. 2001). It should be noted that the uncertainty continues to be high at around $\pm 2\times {10}^{-12}\,{\mathrm{yr}}^{-1}\,{{{\rm{L}}}_{\odot }}^{-1}$.

We can perform the same calculation for M60. The total bolometric luminosity of the sample from Teodorescu et al. (2011), ${L}_{T}=8\times {10}^{10}\,{L}_{\odot }$, and the determined sample size ${n}_{\mathrm{PN}}=2950$ from Section 7 are used to get $\dot{\xi }=1.2\,\times {10}^{-12}\,{\mathrm{yr}}^{-1}\,{L}_{\odot }$. Again, this is lower than the earlier result, $\dot{\xi }=1.7\times {10}^{-12}\,{\mathrm{yr}}^{-1}\,{{L}_{\odot }}^{-1}$ (Teodorescu et al. 2011). The uncertainty is similarly high for M60.

As post-AGB stellar evolution is much faster along the Miller Bertolami tracks than what previous models suggested, one might have expected a higher PN formation rate in order to provide the amount of PNs we observe. However, the final mass range has also changed. Because we are using lower final masses, the central stars move more slowly along their respective tracks in the HR diagram. Ultimately, this leads to a somewhat smaller PN formation rate, just as we can see in our calculations for NGC 4697 and M60. Consequently the new tracks do not help to solve the old discrepancy between stellar death rates of about $15\times {10}^{-12}\,{\mathrm{yr}}^{-1}\,{L}_{\odot }^{-1}$ (Renzini & Buzzoni 1986) and PN formation rates a factor of three (or more) smaller. We still need to assume that most pre-white-dwarf stars in galaxies like NGC 4697 and M60 cannot produce a visible PN, or can at most produce a short-lived one (Méndez et al. 1993).

Méndez & Soffner (1997) found a change of slope between magnitudes $M(5007)=-3.5$ and −2.3. While our simulations do not entirely eliminate this feature, it is much less noticeable as we no longer have any negative slope in that interval. The PNs in this section of the PNLF are mostly on heating tracks and in the curved section joining the heating and cooling tracks. Comparing our distribution of PNs in the HR diagram (Figure 3) with the one in Figure 2 of Méndez & Soffner (1997), it is clear that we now have a much larger amount of PNs with luminosities right below $\mathrm{log}L=3.5$ in the heating tracks and the curved section. The distribution of PNs in this region of the HR diagram is also more uniform. Therefore, it seems that the lack of such PNs in Méndez & Soffner (1997) had been the reason for the change of slope in the PNLF. This can be traced to the existence of a quick drop in luminosity as the H-burning shell is extinguished and the central star enters the white dwarf cooling track.

It is important to note that all the H-burning central stars show the luminosity drop, but it is more pronounced for the more massive ones (see Figure 1). Thus, the adoption of the Miller Bertolami (2016) tracks and substantially lower central star masses leads to predicting a monotonically increasing PNLF for an old population, which seems to be confirmed by NGC 4697. However, other galaxies, like the LMC and SMC, do show a "camel-shaped" PNLF (Jacoby & De Marco 2002; Reid & Parker 2010a). This kind of shape might be explained by involving slightly more massive central stars, but also by a relative lack of lower-mass stars (see Section 5 of Méndez et al. 2008b). We can test this by applying our PNLF simulations to the LMC.

Reid & Parker (2010a) present a sample of PNs in the LMC over a large range of magnitudes. We used their sample of 577 PNs in the central 25 deg² region of the LMC (Tables 1 and 2 in the corresponding catalog, Reid & Parker 2010b). For the transformation into absolute λ5007 magnitudes, we adopted their calculated distance modulus of 18.50 and an extinction correction of 0.46 mag, based on $E(B-V)=0.13$ (Massey et al. 1995; Soffner et al. 1996).

Instead of the observed LMC PNLF, corrected with an average foreground extinction, we could have taken the PNLF with each PN corrected for individual internal extinction, also given by Reid & Parker (2010a); but both LMC PNLFs, shown by Reid & Parker in their Figures 6 and 7, show very similar camel shapes. Therefore, for a discussion of PNLF shapes it should not matter which one is used for the comparison.

Fitting the simulated PNLF to the observed one, we found that we can obtain the valley in the PNLF between $M(5007)=-3.5$ and −0.5 (see Figure 12), similar to what we can see in Reid & Parker (2010a). As expected from our earlier discussion, we had to increase both the minimum and maximum final masses to 0.55 and 0.59 to get this "camel shape." As there are complicating factors in this case (possible differences in the age of the population and in extinction corrections), we do not expect a perfect agreement, but are satisfied that the "camel shape" in the LMC can be reproduced by changing just a few parameters. We note in passing that Badenes et al. (2015) have reported that the LMC PNs show a limited range of ages. The lack of a very old stellar population would seem to agree with the lack of very low-mass central stars we need to reproduce the camel shape.

Zoom In Zoom Out Reset image size

It should be clear that we are far from knowing all the details that contribute to generate the PNLF in any stellar population. To mention just a few, we have the effects of binary evolution and the existence of a minor (but significant) fraction of He-burning central stars in addition to the H-burners. Such complications are very hard to model in Monte Carlo style. To all this we need to add the possibility of different amounts of internal dust extinction, which affects the maximum final mass required to fit the PNLF. The full potential of the PNLF for studies of stellar populations will probably come only after careful testing of PNLF distances and more individual studies of bright PNs in different galaxies, which will likely require 30 meter telescopes.

We have explored the influence of the new post-AGB evolutionary tracks by Miller Bertolami on the interpretation of the observed PNLFs in other galaxies. For that purpose we have revised an earlier procedure for numerical PNLF simulations (Méndez & Soffner 1997). Some adjustments were made to several routines without altering the basic philosophy of the method, which is to rely on empirical information about nebular and stellar properties. A comparison of our simulations with the observed PNLF of M31's bulge, whose distance is very well known, shows that we need to allow for optically thin PNs. Having established the distribution of absorbing factors μ with M31, we selected a galaxy with a larger sample size, NGC 4697, to get a better estimate for the maximum final mass. We clarify the meaning of the maximum final mass by remarking that it does not preclude the presence of more luminous and massive central stars within PNs suffering internal dust extinction. We briefly illustrate the consequences of adopting SBF distances as a way of remarking on the importance of solving the discrepancy between PNLF and SBF distances. We also report on the effect of the new evolutionary tracks on the calculation of PN formation rates. Finally, we show that by adjusting the minimum final mass that leads to a visible PN, we can reproduce the "camel shape" exhibited by PNLFs in the Magellanic Clouds.

L.M.V. acknowledges partial support from the German Academic Scholarship Foundation and would like to thank the Institute for Astronomy, University of Hawaii at Manoa, for their kind hospitality during the course of this research. We thank the anonymous referee for many useful comments.