THE PROGENITORS AND LIFETIMES OF PLANETARY NEBULAE

Despite 250 years of astronomical observations and decades of theoretical work, our understanding of planetary nebulae (PNe) remains rather poor. In the traditional theoretical picture, the progenitors of PNe are low- to intermediate-mass single stars (e.g., Abell & Goldreich 1966; Balick 1987), which eject a large fraction of their envelope at the end of the asymptotic giant branch (AGB) phase. Intense photoionization of the ejected material from the exposed stellar core then leads to the formation of an extended emission nebula. The widely held belief that the Sun will someday become a PN is rooted in this picture, as is the implication that all single stars that go through an AGB phase will eventually form PNe. Yet, there is no observational evidence for this. In fact, the single star hypothesis cannot account for the non-spherical morphologies of many PNe, the observed low rate of PN formation per unit stellar mass, or the occasional appearance of PNe in old stellar systems like globular clusters (Jacoby et al. 1997; Buzzoni et al. 2006; Moe & De Marco 2006; De Marco 2009). To explain these and other inconsistencies, new paradigms of PN formation have emerged, invoking rotation and magnetic fields in single stars (Garcia-Segura et al. 1999), the interaction of AGB winds with close binary companions (Soker 1997), and large mass transfer rates during common envelope episodes (Ciardullo et al. 2005; Nordhaus et al. 2007). All of these mechanisms, or any combination of them, could work in principle, but without robust connections between PNe and their stellar progenitors, and reliable measurements of mean PN lifetimes, it is difficult to determine which ones are at play, and what kinds of PNe they produce. Unfortunately, much of what we know about PNe as a class comes from either case-by-case studies of individual objects or the analysis of small, heterogenous samples that lack the statistical rigor needed to probe these issues.

Here we present the first statistical analysis of the relationship between PNe and their parent stellar populations in the LMC. At a well-determined distance of 50 kpc (Pietrzyński et al. 2013), our closest Galactic neighbor is the ideal setting for this kind of study. On the one hand, the population of PNe in the LMC has been examined in great detail, culminating in a catalog of 435 spectroscopically confirmed objects that extend more than 7 mag down the [O iii] λ5007 PN luminosity function (Reid & Parker 2010; Reid 2014; see Figure 1). This catalog of PNe is virtually free of interlopers, and is statistically complete to a flux limit of 10⁻¹⁴ erg cm⁻² s⁻¹ in [O iii] $\lambda 5007$ (i.e., a line luminosity of 10³³ erg s⁻¹).⁴ On the other hand, the LMC has the best-studied stellar population in the Local Group. In particular, the stellar age distribution (SAD) has been mapped across the entire galaxy through a combination of ground-based and Hubble Space Telescope (HST) observations of millions of individual stars, and their comparison to theoretical isochrones (Harris & Zaritsky 2009). This SAD map consists of 1376 cells covering the inner 64 deg² of the galaxy, with a spatial resolution of 12' × 12' (350 × 350 pc at 50 kpc) and a temporal resolution of 0.2 dex (16 bins spanning lookback times between 6 Myr and a Hubble time).

This SAD map can provide unique insights into the properties of the stellar progenitors of PNe. The SAD is the localized version of the global star formation history, which gives the amount of stellar mass formed as a function of lookback time for an entire galaxy. The stars found today in each individual region of the LMC have drifted there over the years, after having been formed at many different times and locations throughout the galaxy, but the SAD of the region still provides a complete census, within observational uncertainties, of its current stellar content. As such, it must include the progenitors of all astronomical transients found in said region (including PNe), provided that the lifetimes of these transients are short compared to the dynamical drift timescale.

From a statistical point of view, the relationship between PNe and their stellar progenitors is encapsulated in the Delay Time Distribution (DTD), defined as the rate of production of PNe as a function of time after a hypothetical brief burst of star formation. Essentially, the DTD is the impulse response, or Green's function, of the formation rate of PNe, and as such it reveals the evolutionary timescales associated with PN progenitors. Although the DTD can be a key diagnostic for objects with unknown or uncertain progenitors, such as Type Ia supernovae (see Maoz et al. 2014), it is usually difficult to recover from extragalactic data sets. Our work showcases the advantages of deriving DTDs in the Local Group, where we have access to resolved stellar populations and highly complete object catalogs. For the case at hand, the current ($t={{t}_{0}}$) PN formation rate, ${{R}_{i}}({{t}_{0}})$, in a certain region i ($i=1...K$) within the LMC is given by the convolution of ${\Psi }(t)$, the main-sequence turnoff rate of a single-age stellar population formed at time t (which has units of turnoff stars per year per solar mass formed), with the region's SAD, ${{\dot{m}}_{i}}(t)$ (which has units of total mass of stars formed per year):

$$\begin{eqnarray}&&{{R}_{i}}({{t}_{0}})=\int _{0}^{{{t}_{0}}}{{\dot{m}}_{i}}({{t}_{0}}-t){\Psi }(t)dt.\end{eqnarray} \tag{ 1 }$$

Recall that the SAD includes all stars that are now in region i, irrespective of where and when they were formed. The number of PNe expected to be observed in the region, ${{\lambda }_{i}}({{t}_{0}})$, is the product of the integrand in this expression and the mean lifetime of PNe, T_PN(t), expressed as a function of the main-sequence lifetime of their progenitors t:

$$\begin{eqnarray}&&{{\lambda }_{i}}({{t}_{0}})=\int _{0}^{{{t}_{0}}}{{\dot{m}}_{i}}({{t}_{0}}-t){\Psi }(t){{T}_{{\rm PN}}}(t)dt.\end{eqnarray} \tag{ 2 }$$

Note that T_PN(t) can equal zero if stars of a particular main-sequence lifetime t (i.e., formed with a particular mass) do not go through a PN phase.

Given a set of observed PNe in each region ${{N}_{{\rm PN},i}}$ of the LMC, and an SAD map of the galaxy ${{\dot{m}}_{i}}(t)$, recovering the DTD (${\Psi }(t){{T}_{{\rm PN}}}(t)$) is a typical inverse problem. To solve it, we follow the procedure described in Maoz & Badenes (2010), re-casting the convolution in Equation (2) as a discrete sum over time-binned versions of ${{\dot{m}}_{i}}({{t}_{0}}-t)$ and ${\Psi }(t){{T}_{{\rm PN}}}(t):$

$$\begin{eqnarray}&&{{\lambda }_{i}}=\mathop{\sum }\limits_{j=1}^{M}{{m}_{i,j}}{{({\Psi }{{T}_{{\rm PN}}})}_{j}}.\end{eqnarray} \tag{ 3 }$$

Here, the number of PNe expected in each spatial cell i of the SAD map, ${{\lambda }_{i}}$ ($i=1...K$) is the product of the discretized SAD in the cell (${{m}_{i,j}}$, the stellar mass formed during time interval j, with $j=1...M$) and the discretized DTD (${{{\Psi }}_{j}}$, the number of stars that turn off the main-sequence per unit formed stellar mass, per unit time during time interval j, multiplied by the mean lifetime ${{T}_{{\rm PN},j}}$ of PNe formed from stars with main-sequence lifetimes in the time interval j). The values of ${{({\Psi }{{T}_{{\rm PN}}})}_{j}}$ are adjusted with the Markov Chain Monte Carlo solver emcee (Foreman-Mackey et al. 2013), so that random realizations of a Poisson process with expectation value ${{\lambda }_{i}}$ give, on average, the best fit to the actual number of PNe observed in each cell, ${{N}_{{\rm PN},i}}$. The value of ${{N}_{{\rm PN},i}}$ ranges between 0 and 7 in the 732 cells in the SAD map of Harris & Zaritsky (2009), which overlap the area of the PN catalog (see Figure 1).

Zoom In Zoom Out Reset image size

Statistical 1 and 2σ error intervals on each ${{({\Psi }{{T}_{{\rm PN}}})}_{j}}$ are calculated by applying a highest density criterion to the posterior probability distributions generated by emcee. To calculate the additional errors due to the uncertainty in the SADs, we repeat the procedure using the upper and lower limits to the SAD of each cell and temporal bin of the map, and compute the difference between the new median values and those obtained from the best-fit solution. Statistical and SAD-related errors are then added in quadrature. The temporal resolution of the recovered DTD is determined by the size of the observational sample (i.e., the 435 PNe), as there is a fixed amount of "signal" to be spread among the M temporal bins. For the LMC PNe, we have found that a solution with M = 6 provides the best compromise between temporal resolution and significance of the detected progenitor populations, given the statistical and SAD uncertainties. As a verification test, we also solved for ${{({\Psi }{{T}_{{\rm PN}}})}_{j}}$ using the adaptive grid-based search of parameter space with uncertainties based on Monte Carlo simulations described in Maoz & Badenes (2010), and obtained results similar to those described below.

The recovered DTD for the LMC PNe is presented in Figure 2 and Table 1. Its most prominent feature is the clear detection of a population of PN progenitors with a narrow range of main-sequence lifetimes between 5 and 8 Gyr. These lifetimes correspond to zero-age main-sequence masses between 1.2 and 1.0 ${{M}_{\odot }}$ in isolated stellar models of LMC metallicity (Bertelli et al. 2008, 2009). The DTD also suggests the presence of a second population of progenitors with main-sequence lifetimes between 35 and 800 Myr (masses between 8.2 and 2.1 ${{M}_{\odot }}$), but when the SAD uncertainties are taken into account, the formal significance of this detection drops below $2\sigma$. Nevertheless, two separate lines of evidence lend support to this tentative detection. First, a few well-observed LMC PNe have central core masses above ∼0.7 ${{M}_{\odot }}$ (Villaver et al. 2007), which, through the initial- to final-mass relation of Kalirai et al. (2008), require progenitors with main-sequence lifetimes below 800 Myr (in some cases, as low as 60 Myr; Dopita et al. 1993). Second, the spatial correlation between the distribution of formed stellar mass in each DTD age bin and the locations of PNe in the LMC is strong enough to be seen by eye in both the 35–800 Myr and the 5–8 Gyr bins (see Figure 3). For all other main-sequence lifetimes, the correlation is less clear, and we obtain only 2σ upper limits to ${{({\Psi }{{T}_{{\rm PN}}})}_{j}}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Delay Time Distribution of PNe in the LMC

Main-sequence Lifetime	Stellar Massa	${{({\Psi }{{T}_{{\rm PN}}})}_{j}}$	T $_{{\rm PN},j}$
(Myr)	(${{M}_{\odot }}$)	(PNe / 106 ${{M}_{\odot }}$)	$({\rm kyr})$
$\lt 35$	$\gt 8.2$	$\lt 5.8$	$\lt 14$
35–800	$8.2-2.1$	$1.4_{-1.0}^{+0.7}$	$11_{-8}^{+6}$
800–2000	$2.1-1.6$	$\lt 0.5$	$\lt 9$
2000–5000	$1.6-1.2$	$\lt 0.2$	$\lt 5$
5000–8000	$1.2-1.0$	$0.7_{-0.2}^{+0.1}$	27 ± 6
$\gt 8000$	$\lt 1.0$	$\lt 0.1$	$\lt 7$

^aStellar masses that correspond to the main-sequence lifetimes listed in the first column for isolated stellar models of LMC metallicity (Bertelli et al. 2008, 2009).

Download table as:  ASCIITypeset image

We emphasize that these results do not depend on the details of the global star formation history of the LMC or the SAD in each individual region, which is effectively a nuisance function that is disentangled from the DTD in our analysis. Although dynamical processes will mix the LMC's stellar populations on timescales of a few disk crossing times (∼200 Myr; Bastian et al. 2009), the measured SAD for each region i gives the actual distribution of stellar ages that have produced the PNe in that specific location, and therefore the question of when, or where in the galaxy, those stars were formed is irrelevant.

The mean lifetimes of the PNe produced in each time bin of the DTD can be calculated by dividing the recovered ${{({\Psi }{{T}_{{\rm PN}}})}_{j}}$ for each time bin j by the rate at which stars from a coeval population with age t leave the main-sequence, ${\Psi }(t)$. This rate can be derived as follows. Let us assume a power-law dependence of main-sequence lifetime on initial mass, ${{t}_{{\rm ms}}}={{t}_{1}}{{(m/{{M}_{\odot }})}^{\beta }}$, where t₁ is the lifetime, in Gyr, of a 1 ${{M}_{\odot }}$ star. If the mass-normalized initial mass function is also a power law between mass m and the minimum stellar mass, ${{m}_{{\rm min} }}$ (i.e., $dN/dm=(\alpha +2){{({{m}_{{\rm min} }}/{{M}_{\odot }})}^{-(\alpha +2)}}{{(m/{{M}_{\odot }})}^{\alpha }}$), then the main-sequence turnoff rate, in stars per Gyr per unit stellar mass formed, will be

$$\begin{eqnarray}&&{\Psi }(t)=\frac{\alpha +2}{\beta }{{\left( \frac{{{m}_{{\rm min} }}}{{{M}_{\odot }}} \right)}^{-(\alpha +2)}}{{\left( \frac{{{t}_{1}}}{{\rm Gyr}} \right)}^{-\frac{\alpha +1}{\beta }}}{{\left( \frac{t}{{\rm Gyr}} \right)}^{\frac{\alpha +1}{\beta }-1}}.\end{eqnarray} \tag{ 4 }$$

From stellar evolution models, $\beta \approx -2.5$ and, for stars with LMC metallicities, ${{t}_{1}}$ = 8 Gyr (Bertelli et al. 2008, 2009). If, for consistency with the SAD maps, we adopt a Salpeter initial mass function ($\alpha =-2.35$, ${{m}_{{\rm min} }}=0.1$ ${{M}_{\odot }}$), we obtain a turnoff rate of ${\Psi }(t)=0.020{{(t/{\rm Gyr})}^{-0.46}}$ stars per Gyr per solar mass formed. Dividing ${{({\Psi }{{T}_{{\rm PN}}})}_{j}}$ by ${{{\Psi }}_{j}}$ gives the effective lifetime of PNe in each DTD bin. This lifetime is a mean over the PN population in a given age bin. If, for example, only half of the stars in an age bin (say, those in close binaries) actually produce PNe, the lifetimes of those PNe will be twice the value calculated in this way. This method yields lifetimes of 27 ± 6 kyr for the PNe produced by the older progenitors, and $11_{-8}^{+6}$ kyr for the PNe produced by the younger ones (see Figure 2 and Table 1). These numbers are in rough agreement with estimates based on a combination of radiation-hydrodynamics PN models and observed nebular expansion velocities of local objects (Jacob et al. 2013).

The DTD presented here allows us, for the first time, to examine the properties of PN progenitors in a context where the observational sample of PNe is highly pure and complete, and the entire underlying stellar population has been taken into account. This has important implications for our understanding of the stellar evolution of low- and intermediate-mass stars. We confirm the theoretical expectation that some stars do not make a detectable contribution to the LMC PN population, either because they explode as core collapse supernovae (t ≲ 35 Myr, M ≳ 8.0 ${{M}_{\odot }}$) or because their post-AGB evolutionary timescales are longer than the timescale for the dissipation of the ejected envelope (t ≳ 8 Gyr, M ≲ 1.0 ${{M}_{\odot }}$; Herwig 2005). Surprisingly, we find no contribution to the PN population from stars with main-sequence lifetimes between 800 Myr and 2 Gyr, which must evolve through an AGB phase of some kind (Herwig 2005). If these stars produce PNe, they must be either extremely faint (with fluxes below the completeness limits of most PN surveys in the Milky Way), or short-lived (with mean lifetimes below 9 kyr), or both.

Our DTD and derived PN lifetimes result in an integrated PN formation rate of ∼0.02 yr⁻¹ in the surveyed area, which includes ∼80% of the stellar mass of the LMC (Harris & Zaritsky 2009). Taking an absolute magnitude of ${{M}_{V}}=-18.4$ for the LMC (de Vaucouleurs et al. 1991) and a bolometric correction of −0.8 (Buzzoni et al. 2006), this translates into a bolometric-luminosity specific PN formation rate of $\sim 7\times {{10}^{-12}}$ PNe yr⁻¹ L $_{\odot }^{-1}$, which is lower than the typical values obtained in studies of Galactic PNe. However, Milky Way surveys must deal with heterogeneous, incomplete samples of questionable purity, and are affected by the systematic uncertainties associated with poorly known distances (Parker et al. 2006; Sabin et al. 2014, but see Frew et al. 2015). Because all LMC PNe are at a well-determined distance, and have known foreground reddenings, our analysis is largely free from these problems.

Among the PNe visible in the LMC today, $40_{-29}^{+23}\%$ are generated by the younger progenitors and 25 ± 8% by the older progenitors. Thus, roughly one third of the LMC PNe could come from stars in other age bins, for which our DTD only gives upper limits. Nevertheless, the presence of two distinct and disjoint progenitor populations strongly suggests the existence of two separate formation channels for PNe, a possibility that has been qualitatively discussed in the past (Moe & De Marco 2006; De Marco 2009; Frew & Parker 2010), but never observationally confirmed.

If we restrict our analysis to the PNe in the brightest third of the sample (L $_{[{\rm O}\,{\rm III}]}\gtrsim 4\times {{10}^{34}}$ erg s⁻¹), the younger progenitor population is not detected. This suggests, somewhat counter-intuitively, that it is the older, less-massive progenitors that produce the brighter PNe. If true, this might explain a long-standing issue with the use of the [O iii] PN luminosity function for calculating extragalactic distances. This method is well-calibrated in more than a dozen galaxies with known Cepheid distances (Ciardullo et al. 2002), yet it yields distances to the Virgo and Fornax clusters that are ∼10% smaller than those found with other methods (Jacoby et al. 1990; McMillan et al. 1993; Freedman et al. 2001). In the past, this offset has been attributed to the presence of an unidentified systematic error which only affects galaxies beyond ∼10 Mpc (Ferrarese et al. 2000), but a more physical explanation is that the brightest PNe in old elliptical galaxies are intrinsically more luminous than their spiral galaxy counterparts (Ciardullo 2012). This is consistent with our results.

A detailed interpretation of our DTD in the framework of specific PN formation mechanisms is beyond the scope of the present work. Nevertheless, we can still put together a basic picture using some additional information. The scale height of bipolar PNe in the Milky Way is significantly smaller than that of round or elliptical PNe (Phillips 2001; Parker et al. 2006), implying that they have younger progenitors (Corradi & Schwarz 1995). In the LMC, HST has imaged 68 of the 435 PNe in our sample, with 49 objects showing symmetric (round or elliptical) morphologies and 19 showing distinctly bipolar structures (Shaw et al. 2006). The DTD of the symmetric subsample of PNe shows a significant detection in the 5–8 Gyr bin (as suggested by Stanghellini 2009), while the subsample of bipolar PNe is too small to produce any statistically significant detections. We have also derived the DTD of the 124 PNe that show enhanced N in their spectra (W. Reid 2015, private communication). These "Type I" PNe (Peimbert 1978; Peimbert & Torres-Peimbert 1983) are associated with core masses larger than 0.64 ${{M}_{\odot }}$ (Kaler & Jacoby 1990), corresponding to main-sequence masses above ∼2.25 ${{M}_{\odot }}$. The DTD shows a marginal (1.7σ) detection in the 35–800 Myr bin, consistent with this minimum core mass. Taken together, these results suggest a scenario wherein the 35–800 Myr progenitors produce fainter, more asymmetric PNe, possibly of binary origin (Nordhaus et al. 2007), while the 5–8 Gyr progenitors produce brighter, more symmetric PNe, perhaps via traditional single star evolution.

This basic picture is in rough agreement with the properties of individual PNe in the LMC (Dopita et al. 1993; Villaver et al. 2007), and the Galactic Bulge (Gesicki et al. 2014), but a rigorous statistical validation would require a complete morphological census. The main features of our proposed scenario are also supported by theoretical work. On the one hand, the PN lifetimes stemming from the older progenitors are consistent with state-of-the-art radiation-hydrodynamics models for single-star PNe (Jacob et al. 2013). On the other hand, it has been shown that binary systems that enter a common envelope phase while the primary is on the AGB can form bipolar PNe, provided that the primary has a minimum mass of ∼2 M $_{\odot }$ (Soker 1998). We leave a more detailed evaluation of these and other specific PN formation mechanisms to future work.

We conclude by pointing out that the methods and techniques presented here are applicable to a wide variety of astronomical objects, from variable stars (Cepheids, RR Lyrae, δ Scuti, etc.) to interacting binaries (novae, cataclysmic variables, high- and low-mass X-ray binaries, etc.). The best available catalogs for these objects, and the most complete census of their progenitor stellar populations, are found in Local Group galaxies. The recovery of high-quality DTDs from these data sets can provide much-needed tests for specific stellar evolution scenarios.

We acknowledge useful discussions with Avishay Gal-Yam, Dan Foreman-Mackey, Dustin Lang, and Dennis Zaritsky. Letizia Stanghellini and Warren Reid kindly provided us with updated lists of PN morphologies from HST and PNe with N-rich spectra for the revised version of the paper, which greatly benefitted from the suggestions of the anonymous referee. D. M. acknowledges support by the I-CORE Center of the Planning and Budgeting Committee and the Israel Science Foundation. The Institute for Gravitation and the Cosmos is supported by the Eberly College of Science and the Office of the Senior Vice President for Research at the Pennsylvania State University.