The PICS project: II. Circumnebular extinction variations and their effect on the planetary nebula luminosity function

For decades, the theoretical understanding of planetary nebulae (PNe) has remained in tension with the observed universal bright-end cutoff of the PN luminosity function (PNLF). The brightest younger PN populations have been observed to be fainter in their [OIII] emission than expected. Recent studies have proposed that circumnebular extinction is a key ingredient in bringing their brightness down to the observed level. In this work, we use the recently introduced PICS (PNe In Cosmological Simulations) framework to investigate the impact of different circumnebular extinction treatments on the modeled PNe and their PNLF for a large range of stellar ages and metallicities. We test how different slopes in the observed relation of extinction versus central star mass modify the bright-end cutoffs of the PNLF, finding that steeper slopes lead to large changes for young stellar populations. In contrast, the differences for older PNe are much smaller. However, for individual PNe, the extinctions observed in nearby galaxies appear to be much higher than the models predict, showing that improvements on both the modeling and observational sides are needed to gain a better understanding of the brightest and strongly extincted PNe. These findings further advance the theoretical foundation for interpreting observed extragalactic PN populations coming from more complex composite stellar populations in the future.

1 Introduction

The bright-end cutoff of the planetary nebula luminosity function (PNLF) has long been used as a standard candle in the cosmic distance ladder (Jacoby, 1989; Rekola et al., 2005). This has been based on the empirical finding that the brightest planetary nebulae (PNe) in a given galaxy have the same [O III] $\lambda 5007$ emission line absolute magnitude of $M*(5007)=-4.5$, independent of galaxy morphology (Ciardullo et al., 2002). This universal PNLF bright-end cutoff appears to be independent of stellar age and metallicity above $12+[\mathrm{O}/\mathrm{H}]\gtrsim 8.5$ (Ciardullo, 2012). This universality has stood in contrast with theoretical expectations for a long time, where younger stellar populations with more massive PN central stars should feature much brighter PNe than older, less massive central stars (Ciardullo, 2022).

One proposal for bringing the brightness down for PNe in younger stellar populations has been self-extinction through the dust formed in the winds during the asymptotic giant branch (AGB) phase preceding the PN phase (Ciardullo and Jacoby, 1999; Ventura et al., 2014). Based on their findings and model, Jacoby and Ciardullo (2025) showed how this circumnebular extinction seems to exactly compensate for the increasing potential brightness of PNe with more massive central stars over a large mass range, which could, in part, explain the constant PNLF bright-end cutoff. Considering the brightest PNe in the bulge of M 31, Davis et al. (2018) found very large circumnebular extinctions, implying intrinsic brightnesses well beyond expectations for an older stellar population. Their results are also consistent with studies of other nearby galaxies (Méndez et al., 2005; Reid and Parker, 2010). Furthermore, circumnebular extinction has been suggested as a way to understand the connection between UV measurements of post-AGB stars and PN central star candidates (Sarzi et al., 2011). The findings demonstrate that understanding the effects of circumnebular extinction on the bright end of the PNLF and individual PNe is crucial to properly interpret observations in both local and distant galaxies.

Recently, theoretical modeling of PNe has made large advances through the novel PICS (PNe In Cosmological Simulations) framework (Valenzuela et al., 2025a; Valenzuela et al., 2025b), a model that traces all the properties necessary to obtain a realistic PN population for any given single stellar population (SSP). Built in a modular fashion, it allows one to study precisely how variations in the models and parameters impact the resulting PN and PNLF properties for any stellar population.

In this work, we use PICS to investigate how the chosen circumnebular extinction recipe affects the bright-end cutoff of the PNLF for different SSPs and compare the individual modeled extinctions of the bright PNe to observations. In Section 2, we provide a brief overview of the PICS model and the extinction recipes tested. We then present the results in Section 3 and discuss them in Section 4.

2 Methods

We use the fiducial model of the PICS framework presented by Valenzuela et al. (2025b), in which we test how the resulting PN properties change based on variations in the circumnebular extinction model. In the following, we provide a brief overview of PICS and the variations of the circumnebular extinction recipe we tested.

2.1 PICS model

The PICS framework is composed of several modules that, for a given SSP, predict the PN population found therein. Such an SSP is parameterized by its total mass, age, metallicity, and initial mass function (IMF). In the PICS fiducial model, these properties are first used to determine the initial mass of the central stars of the current PNe using the lifetime function (Miller Bertolami, 2016). The final stellar mass is determined using the initial-to-final mass relation (IFMR) and stellar properties derived from the post-AGB evolutionary tracks of Miller Bertolami (2016). All three of these steps are metallicity-dependent. Based on the stellar properties, the nebular model by Valenzuela et al. (2019) is applied together with a metallicity correction from Dopita et al. (1992). The nebular model also accounts for the possibility that the nebula is optically thin to ionizing photons. The primary resulting quantity is the [OIII] flux, $F(5007)$.

As shown by Jacoby and Ciardullo (2025), the [OIII] flux can be significantly extincted by circumnebular dust that was expelled through the AGB winds. The extinction is expected to decrease over time during the post-AGB phase as the dust and nebula diffuse into the interstellar medium (ISM) with age. The fiducial PICS model accounts for a post-AGB time-dependent circumnebular extinction consistent with the empirical findings from Jacoby and Ciardullo (2025), where we define the post-AGB time as the time after the central star reaches an effective temperature of $T_{\text{eff}}=25000\text{\hspace{0.17em}K}$. The implemented extinction recipe is solely a function of the post-AGB time of the central star, which, in particular, means that we assume it to be independent of stellar mass.

To derive this equation, an intermediate step has to be taken because the relations from (Jacoby and Ciardullo (2025); Table 2) between stellar final mass and extinction, $c(\mathrm{H}\beta )$, are derived only from the brightest PNe. For this, we determined the post-AGB time at which a central star (with a metallicity of $Z=0.01$) of a given final mass reaches its brightest intrinsic [OIII] magnitude. This will be needed to derive an “equivalent brightest PN final mass.” The relation by Jacoby and Ciardullo (2025) returns an extinction based on the final mass of the brightest PN. From these two relations, we derived a direct dependence of extinction on post-AGB time. From the post-AGB time, we obtain the “equivalent brightest PN final mass,” from which we finally obtain the extinction. Because we make the simplifying assumption that the extinction depends only on post-AGB time and not on stellar mass, we use the resulting relation for all PNe. For the fiducial PICS model, we use their orthonormal regression fit to the combined Large Magellanic Cloud (LMC) and M 31 sample of PNe, the Combined OR (orthogonal regression) fit (from the table of Jacoby and Ciardullo, 2025).

By running PICS for a given SSP, the result is a population of PNe, each with an [OIII] magnitude. The resulting PNe differ from each other because the nebular model by Valenzuela et al. (2019) is stochastic. Two reasons for this are the randomly drawn post-AGB lifetime and the randomness of the optical thickness of the nebula. In this work, we keep the fiducial PICS model fixed, with the exception of the circumnebular extinction treatment. The variations are described in the following sections.

2.2 Circumnebular extinction variations

While we used the Combined OR fit of Jacoby and Ciardullo (2025) for the fiducial PICS model, they also provided five additional fits to their PN data, varying both the subsample of PNe (Combined, LMC, and M 31) and the fitting method (OR; ordinary least squares, OLS). They also include a fit that optimizes the final mass range over which the bright-end cutoff magnitude remains constant at $M*(5007)\sim -4.5$ (named Optimized, see more details in Section 5 of Jacoby and Ciardullo, 2025).

For the Combined OR fit, we assumed a metallicity of $Z=0.01$ for the post-AGB evolution needed for the intermediate step (see the previous section). However, for the pure M 31 and LMC fits, it is more appropriate to assume metallicities that would more closely match those expected for their respective brightest PNe. Only for the Combined and Optimized fits, we assumed $Z=0.01$. For the LMC fits, we used $Z=0.007$ based on the upper metallicity found for the LMC and its PNe (Dopita and Meatheringham, 1991). For M 31, we used a solar-like metallicity of $Z=0.015$, based on metallicity measurements of luminous PNe in M 31 (Balick et al., 2013). For each of the metallicities, we fit a broken linear fit (above and below $M_{\text{final}}=0.60{\mathrm{M}}_{\odot }$) to the relation of final mass versus post-AGB time at which an optically thick PN reaches its brightest [OIII] luminosity, analogously to the $Z=0.01$ fit by Valenzuela et al. (2025b), Appendix D).

$M_{\text{final,bright}}^{Z=0.007}/{\mathrm{M}}_{\odot }=\max \left\{\begin{array}{c}\hfill 1.023-0.1740\log \left(t_{\text{post}-\text{AGB}}/\mathrm{y}\mathrm{r}\right)\\ \hfill 0.714-0.0432\log \left(t_{\text{post}-\text{AGB}}/\mathrm{y}\mathrm{r}\right)\end{array}\right\},(1)$

$M_{\text{final,bright}}^{Z=0.015}/{\mathrm{M}}_{\odot }=\max \left\{\begin{array}{c}\hfill 0.956-0.1477\log \left(t_{\text{post}-\text{AGB}}/\mathrm{y}\mathrm{r}\right)\\ \hfill 0.743-0.0531\log \left(t_{\text{post}-\text{AGB}}/\mathrm{y}\mathrm{r}\right)\end{array}\right\},(2)$

The derived relations for the two metallicities are given in Equations 1, 2: where the respective maximum of the two linear terms is taken. Based on the fit parameters presented in Table 2 of Jacoby and Ciardullo (2025), the extinction c(Hβ) is finally determined according to Equation 3,

$c\left(\mathrm{H}\beta \right)=\mathrm{S}\times M_{\text{final,bright}}+I,(3)$

where $S$ is the slope and $I$ is the intercept. One important note is that $M_{\text{final,bright}}$ should not be confused with the actual final mass of a PN. It is rather the “equivalent brightest PN final mass” as described in Section 2.1, based on the post-AGB time of a given PN (which can have any other final mass), and needed as input for the fits by Jacoby and Ciardullo (2025).

In addition to the relations taken directly from Jacoby and Ciardullo (2025), we also considered a shift of the OR extinction relations of $\mathrm{\Delta }c(\mathrm{H}\beta )=+0.08$. The motivation for this shift is introduced in Section 3.1 and is explained there in detail. Finally, the extinction in $\mathrm{H}\beta$ is converted to $A_{5007}$ according to Equation 13 of Valenzuela et al. (2025b).

The result is an [OIII] extinction, $A_{5007}$, as a function of post-AGB time, as shown in Figure 1. All recipes produce a broken linear relation of $A_{5007}$ with respect to the logarithmic post-AGB time, with the break occurring at around 200 years. The OR models have significantly higher extinction than their OLS counterparts at young ages but in turn have slightly lower extinction at ages beyond $\sim$1,500 years. The modified OR models with the $+0.08$ shift have a slight offset to larger extinction by $\mathrm{\Delta }{A}_{5007}=+0.18$. Between the different underlying observed PN samples, the LMC relations show the most extinction until about 1,000 years, with the LMC OLS relation being similar to the M 31 OR fit. The Combined relations, respectively, lie between the LMC and M 31 relations but align more closely with M 31 at later times. The Optimized relation leads to slightly stronger extinctions than the M 31 OR model.

Figure 1

Figure 1. Circumnebular extinction of [OIII], $A_{5007}$, as a function of post-AGB time of the central star for the different extinction recipes tested. The recipes are based on the relations by Jacoby and Ciardullo (2025) for the extinction $c(\mathrm{H}\beta )$ versus core mass for the brightest PNe in M 31 and the LMC. The colors indicate the Combined (blue), LMC (green), M 31 (orange), and Optimized (purple) fits, and the line styles denote the respective OR (solid), shifted OR (dotted), and OLS (dash-dotted) fits for all but Optimized. The dotted OR fit + 0.08 lines are based on their OR relations shifted by $\mathrm{\Delta }c(\mathrm{H}\beta )=+0.08$.

These variations of the circumnebular extinction recipe were implemented within the PICS framework as described. They are used in the following to investigate their effect on the resulting PN population properties.

3 Results

While the recipe for extinction was derived from the maximal brightness of PNe with a given final stellar mass, the nebular model by Valenzuela et al. (2019) allows for escaping ionizing fractions and different [OIII] to $\mathrm{H}\beta$ flux ratios. This leads to the possibility of dimmer PNe arising, leading to a range of possible [OIII] magnitudes for a given final mass and post-AGB time.

3.1 Extinction of $\mathrm{H}\beta$

As the recipe for circumnebular extinction was informed through the observations by Jacoby and Ciardullo (2025), this allows for a direct consistency check for the actually derived extinction values of the resulting PN population from the PICS model. To this end, we ran PICS for an LMC-like PN population of $Z=0.007$, for which we only kept those with extincted magnitudes $M(5007)\le -3.5$ for comparison with the observational sample. We ran PICS with the extinction recipes based on the LMC OR, OLS, and the shifted OR relations. The $c(\mathrm{H}\beta )$ extinctions of the bright PNe are shown in the left panel of Figure 2, where the lines indicate the median extinctions found for the bright PNe at the individual final masses, and the shaded regions indicate the $1\sigma$-scatter. The observed extinctions of the LMC PNe are indicated with black points. The same is shown for an M 31-like PN population of $Z=0.015$ in the right panel for the analogous three M 31 extinction fits to the observed M 31 PNe.

Figure 2

Figure 2. Circumnebular extinction of $\mathrm{H}\beta$, $c(\mathrm{H}\beta )$, as a function of stellar final mass for the brightest observed and modeled PNe assuming different extinction recipes. LMC-like PNe $(Z=0.007)$ are shown on the left and M 31-like PNe $(Z=0.015)$ on the right. Bright PNe are selected according to $M(5007)\le -3.5$, analogously to the observed PNe (black circles) as selected by Jacoby and Ciardullo (2025). The lines indicate the medians, and the shaded regions indicate the $1\sigma$-scatter of the extinction values at the respective final masses. The extinction recipes correspond to the input fit relations from Jacoby and Ciardullo (2025), orthonormal regression (OR), shifted OR by $\mathrm{\Delta }c(\mathrm{H}\beta )=+0.08$, and ordinary least squares (OLS).

For both PN populations, it can be observed that the OR model (light green and light orange) leads to overall lower extinctions of the PNe compared to the observations, although the slope remains the same. In contrast, the OLS model (dashed blue and red) results in too much extinction at lower masses $(\lesssim 0.56{\mathrm{M}}_{\odot })$ and underestimates the extinction at higher masses $(\gtrsim 0.59{\mathrm{M}}_{\odot })$ because of its overall shallower slope.

We find that the OR models lead to systematically too low values of $c(\mathrm{H}\beta )$ by approximately 0.08. For this reason, we implemented a shifted OR relation with $\mathrm{\Delta }c(\mathrm{H}\beta )=+0.08$. For the shifted models, the PN extinction distributions are brought into much better agreement with the observations (dark solid lines in Figure 2).

The motivation for such a shift is that the observations on which the relations are based are of the brightest PNe within one magnitude of the bright-end cutoff and not only of the very brightest PNe. In contrast, our implemented recipe assumes the extinction from the fits of Jacoby and Ciardullo (2025) to correspond to the point in time where the PNe are at their brightest. After the time of reaching its brightest [OIII] flux, a PN decreases in brightness, and the extinction declines. Such PNe could still reside in the brightest magnitude but would have a lower extinction. For this reason, it is expected that the average extinction of the brightest PNe $(M(5007)\le -3.5)$ will be underestimated through the application of the PICS model, suggesting that a correcting shift is necessary. We conclude that the shifted OR recipe leads to the most compatible extinctions with observations and that the shift is reasonable for adaptation within the PICS framework.

3.2 Extinction effects on the PNLF

Having verified the agreement of the modeled circumnebular extinction with the original observations, we next tested how much the variations of the extinction recipe actually affect the PNLF and its bright end. Running PICS for a sufficiently large number of PNe that arise from an SSP of a certain age and metallicity produces the PNLF expected for such an SSP. Figure 3 shows a grid of such SSP PNLFs for five different metallicities from $Z=0.001$ to 0.08 and ten SSP ages between 0.25 and 13 Gyr, using no circumnebular extinction (dotted black lines) and the different extinction recipes as described in Section 2 (solid lines for OR, dotted for shifted OR, and dash-dotted for OLS).

Figure 3

Figure 3. Grid of SSP PNLFs for different metallicities (columns; $Z=0.001$, 0.01, 0.02, 0.04, and 0.08) and ages (rows; $t_{\text{age}}=0.25$ to 13 Gyr for the different circumnebular extinction recipes tested (colored lines) and without extinction (dashed black lines). The observed universal bright-end cutoff of $M(5007)=-4.5$ is indicated by the dashed gray lines.

As expected, the extinction is the highest for the young SSPs with high-mass PN central stars, leading to dimmer bright-end cutoffs for the extincted PNe (colored lines) compared to the non-extincted PNe (dotted black lines) in the upper rows. Here, it becomes clear that the OLS recipes lead to the smallest amount of dimming (dash-dotted lines) for young SSPs, followed by the OR (solid lines) and shifted OR recipes (dotted colored lines). Compared to the Combined OR recipe used in the fiducial PICS model, most extinction recipes lead to brighter PNLFs, except for the two LMC OR recipes. However, these differences become noticeably smaller for ages of 1 Gyr and older, reducing the differences between the circumnebular extinction recipes to $\lesssim$0.5 mag. These differences between the tested recipes remain rather constant over the metallicity range considered ($Z=0.001$–0.08) and also when moving toward older SSPs. For the oldest SSPs of $t_{\text{age}}\gtrsim$8 Gyr, the dimming becomes zero or near-zero for all but the OLS recipes. This is related to the more shallow slope of the OLS relations (see Figure 2), resulting in a generally stronger extinction even at the lower final masses below 0.55${\mathrm{M}}_{\odot }$.

Having been optimized to recover a constant bright-end across a range of masses that is as large as possible, the Optimized PNLFs (purple lines in Figure 3) show a PNLF dimming that is right in the middle of the alternative models for younger SSPs. This leads to a slightly brighter cutoff than $M*(5007)=-4.5$ for some SSPs, but it quickly approaches a negligible dimming for older SSPs when the lower-mass PNe become intrinsically dimmer. Thus, the Optimized recipe produces the most constant PNLF cutoff across SSPs of different ages and metallicities using the PICS model.

3.3 Extinctions of de-reddened PNe

Finally, we investigated the relation of intrinsic de-reddened [OIII] magnitudes with their circumnebular extinctions. For this, we ran PICS for 100,000 PNe at an intermediate metallicity of $Z=0.01$ with different circumnebular extinction recipes, 1,000 PNe each for 100 uniformly spaced ages between 0.25 and 10 Gyr. The relation of their intrinsic [OIII] magnitudes, $M_{\text{dered}}(5007)$, with the [OIII] extinction, $A_{5007}$, is shown in Figure 4. We also included observations of bright PNe from four nearby galaxies, as compiled by Davis et al. (2018). These are from M 31 (Davis et al., 2018), the LMC (Reid and Parker, 2010), NGC 4697 (Méndez et al., 2005), and NGC 5128 (Walsh et al., 2012). The dashed line marks the location of the bright-end cutoff for PNe of different extinctions. In the left and middle panels, the Combined OR and shifted Combined OR recipes were applied, respectively. The elevation to higher extinctions can immediately be observed for the shifted OR recipe: first, the bulk of PNe reach higher values of $A_{5007}$, and second, there are fewer PNe that reside in the regime brighter than the bright-end cutoff (region below the dashed gray line). Overall, most PNe have lower extinction values of $A_{5007}\lesssim 0.5$ and 0.75 for OR and shifted OR, respectively. This stands in contrast to the observed PNe, where the extinctions of the PNe are more uniformly distributed. Only NGC 4,697 contains a large fraction of PNe with $A_{5007}\lesssim 0.5$, but still has a significant amount of them at higher extinctions (3 out of 12), unlike the modeled PNe. Additionally, there is a larger fraction of intrinsically bright observed PNe with $M(5007)\lesssim -5.0$ compared to the modeled PNe.

Figure 4

Figure 4. Circumnebular extinction of [O III], $A_{5007}$, versus the intrinsic de-reddened [OIII] magnitude, $M_{\text{dered}}(5007)$, of observed and modeled PNe for three different circumnebular extinction recipes. Left: Combined OR. Middle: Combined OR with a shift of $\mathrm{\Delta }c(\mathrm{H}\beta )=+0.08$. Right: Combined OR with the same shift of +0.08 and an added normally distributed scatter in $c(\mathrm{H}\beta )$ of $\sigma =0.1$. The values of the observed PNe are taken from Davis et al. (2018), who observed PNe in the bulge of M 31 and compiled the data from other works for three other galaxies: the LMC (Reid and Parker, 2010), NGC 4,697 (Méndez et al., 2005), and NGC 5128 (Walsh et al., 2012). The dashed gray lines indicate the bright-end cutoff of the PNLF of $M(5007)=-4.5$.

Since the observed extinctions are much greater than the modeled extinctions in both cases (left and center panel of Figure 4), we decided to add a random, normally distributed scatter to $c(\mathrm{H}\beta )$ after applying the respective linear relation from Jacoby and Ciardullo (2025). The motivation is that any asymmetric nebulae would be expected to also have varying amounts of dust along different lines of sight, thus changing the circumnebular extinction. We tested an additional Gaussian scatter with $\sigma =0.1$ applied to the shifted Combined OR relation. This chosen value leads to a slightly larger scatter than is observed in the $c(\mathrm{H}\beta )$ relation shown in Figure 2, but it still has the same median behavior as the non-scattered versions. In the right panel of Figure 4, we show how the $A_{5007}$ extinctions can reach higher values due to the increased scatter. However, even with this, the large relative amount of high-extinction PNe found in the observations is still not recovered. Finally, PNe brighter than the bright-end cutoff (those found below the dashed diagonal line) are now, on average, also further away from the cutoff, with several ending up at magnitudes close to $M(5007)\sim -4.8$. Therefore, the larger scatter, on one hand, increases the extinctions that can generally be reached, but on the other hand also leads to overly bright PNe, though there are also occasional observations of such PNe (Jacoby et al., 2024, Figure 59). We conclude that the simple additional scatter for the extinction is not sufficient in bringing the modeled and observed [OIII] extinctions into agreement.

4 Discussion

For the [OIII] circumnebular extinction, $A_{5007}$, the largest effect of the different tested extinction recipes is found at the beginning of the post-AGB phase (Figure 1). Because only the more massive central stars can come near the bright end within this short time frame due to their faster heating (Miller Bertolami, 2016), the only SSP PNLFs that experience significant extinction differences between the recipes are those of the youngest SSPs with $t_{\text{age}}\lesssim$0.5 Gyr (Figure 3). This implies that, in particular, the early behavior of dust build-up in the nebula is of great importance for the resulting amount of circumnebular extinction. However, given that the implemented recipes were based on the empirical results of extinction related to final mass from Jacoby and Ciardullo (2025), we note that the regime of very large extinction is largely beyond the originally fitted range of $c(\mathrm{H}\beta )\lesssim 0.5$ $(A_{5007}\lesssim 1.2)$. In turn, this range covers most of the PN lifetime after 200  years, where the differences are much smaller as they are more constrained by the observations.

Overall, the modeled PNe have similar extinctions to the observed PNe for the different recipes, but neither the pure OR nor OLS recipes reproduce the same median extinction-final mass relation as the observations (Figure 2). This shows that the average extinction relation of the PNe in the brightest magnitude is not the direct desired input relation for the PICS model (Valenzuela et al., 2025b), though the slope of the OR recipe is similar to that of the observed median relation. We addressed this offset by shifting the OR relations up by the respective amount of $\mathrm{\Delta }c(\mathrm{H}\beta )=0.08$. However, in the future, the relation between extinction and final mass will ideally be informed by a larger sample of PNe that lie directly at the bright end, to limit the scatter around the fitted relation (see the discussion by Jacoby and Ciardullo, 2025).

Despite the Combined OR recipes reaching the second highest extinctions (after the LMC OR recipes; see Figure 1) and generally reproducing the relation between extinction and final stellar mass, it is curious that the models using Combined OR do not reach as high extinctions for a full PN population as found in nearby galaxies (Figure 4). Even when including a scatter for the extinction that would produce higher values, the distribution still does not match the observations. There are multiple possible reasons for this discrepancy. For one, the modeled stellar population is artificially created with PNe at uniformly spaced ages and therefore does not represent the true star formation history of any of the four observed galaxies. However, from the model side, the only realistic possibility to obtain a similar amount of high-extinction PNe would be to have a much larger fraction of younger stars, which is likely not the case for the observed galaxies.

Furthermore, potential model-intrinsic issues of the recipes used by PICS could lead to PN populations that do not entirely match the observed PNe. Although the PNLF properties and many relations, such as the extinction-final mass relation, do agree with observations, the properties of individual PNe could still be systematically offset because of unknown degeneracies. Another possible explanation is that PICS is missing relevant channels for forming PNe at the bright end, for example, accreting white dwarfs (Soker, 2006; Souropanis et al., 2023), blue stragglers (Ciardullo et al., 2005), and white dwarf mergers (Yao and Quataert, 2023), which could result in different extinction properties. These additional possible pathways are supported by the evidence that a significant fraction of PNe are in binary systems (Jones and Boffin, 2017; Jacoby et al., 2021; Csukai et al., 2025). Additionally, a random scatter in the final mass through the IFMR, as proposed by Jacoby and Ciardullo (2025), could also increase the number of bright PNe with higher extinctions. Systematic biases in the extinction determination for the observations could lead to a preference for higher extinctions being predicted for the PNe. Then, the slope parallel to the bright-end cutoff line observed in Figure 4 for M 31 and NGC 4697 could also be the result of selection bias coupled with systematically higher extinctions.

Finally, it is interesting to consider the effects that metallicity would have on the extinction, which are not accounted for in the tested recipes. At higher metallicity, there should be more metals to contribute to the formation of dust for equal stellar masses. As also pointed out by Jacoby and Ciardullo (2025), curiously, their derived relations between extinction and final mass show no significant differences between M 31 and the LMC PNe, even though the PNe in the bulge of M 31 have roughly twice the metal mass fraction as those in the LMC. The LMC PNe even seem to reach higher extinctions at a given final mass (Figure 2), though this appears to stand in contrast to the overall lower measured [OIII] extinctions (Davis et al., 2018), as shown in Figure 4.

To conclude, while circumnebular extinction surely plays a central role in limiting the [OIII] magnitude of PNe in intermediate-to-younger stellar populations to the universal bright-end cutoff, the details of its dependency on other physical properties are not straightforward. While recipe variations in the PICS model affect the youngest stellar populations the most, due to being the least constrained observationally, the general discrepancies between the modeled extinctions and those determined from observations make it very difficult to disambiguate between the recipes. Further observations and advances in the understanding of circumnebular dust evolution will be necessary to improve the modeling.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

LV: Writing – original draft, Writing – review and editing, Methodology, Conceptualization, Visualization. GJ: Conceptualization, Writing – review and editing. R-SR: Conceptualization, Writing – review and editing. MM: Writing – review and editing. RM: Writing – review and editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. LV acknowledges support by the German Academic Scholarship Foundation (Studienstiftung des deutschen Volkes) and the Marianne-Plehn-Program of the Elite Network of Bavaria.

Acknowledgements

The authors thank the two referees for their comments that helped improve the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Balick, B., Kwitter, K. B., Corradi, R. L. M., and Henry, R. B. C. (2013). Metal-rich planetary nebulae in the outer reaches of M31. ApJ 774, 3. doi:10.1088/0004-637X/774/1/3

CrossRef Full Text | Google Scholar

Ciardullo, R. (2012). The planetary nebula luminosity function at the dawn of gaia. Ap&SS 341, 151–161. doi:10.1007/s10509-012-1061-2

CrossRef Full Text | Google Scholar

Ciardullo, R. (2022). The planetary nebula luminosity function in the era of precision cosmology. Front. Astronomy Space Sci. 9, 896326. doi:10.3389/fspas.2022.896326

CrossRef Full Text | Google Scholar

Ciardullo, R., and Jacoby, G. H. (1999). The circumstellar extinction of planetary nebulae. ApJ 515, 191–195. doi:10.1086/307025

CrossRef Full Text | Google Scholar

Ciardullo, R., Feldmeier, J. J., Jacoby, G. H., Kuzio de Naray, R., Laychak, M. B., and Durrell, P. R. (2002). Planetary nebulae as standard candles. XII. Connecting the population I and population II distance scales. ApJ 577, 31–50. doi:10.1086/342180

CrossRef Full Text | Google Scholar

Ciardullo, R., Sigurdsson, S., Feldmeier, J. J., and Jacoby, G. H. (2005). Close binaries as the progenitors of the brightest planetary nebulae. ApJ 629, 499–506. doi:10.1086/431353

CrossRef Full Text | Google Scholar

Csukai, A., Zijlstra, A. A., McDonald, I., and De Marco, O. (2025). Central-star extinctions towards planetary nebulae. Mon. Not. R. Astron. Soc. 543, 3035–3054. doi:10.1093/mnras/staf1552

CrossRef Full Text | Google Scholar

Davis, B. D., Ciardullo, R., Jacoby, G. H., Feldmeier, J. J., and Indahl, B. L. (2018). The true luminosities of planetary nebulae in M31’s Bulge: massive central stars from an old stellar population. ApJ 863, 189. doi:10.3847/1538-4357/aad3c4

CrossRef Full Text | Google Scholar

Dopita, M. A., and Meatheringham, S. J. (1991). Photoionization modeling of magellanic cloud planetary nebulae. II. ApJ 377, 480. doi:10.1086/170377

CrossRef Full Text | Google Scholar

Dopita, M. A., Jacoby, G. H., and Vassiliadis, E. (1992). A theoretical calibration of the planetary nebular cosmic distance scale. ApJ 389, 27. doi:10.1086/171186

CrossRef Full Text | Google Scholar

Jacoby, G. H. (1989). Planetary nebulae as standard candles. I. Evolutionary models. ApJ 339, 39. doi:10.1086/167274

CrossRef Full Text | Google Scholar

Jacoby, G. H., and Ciardullo, R. (2025). The critical role of dust on the [O III] planetary nebula luminosity function’s bright-end cutoff. ApJ 983, 129. doi:10.3847/1538-4357/adc0fb

CrossRef Full Text | Google Scholar

Jacoby, G. H., Hillwig, T. C., Jones, D., Martin, K., De Marco, O., Kronberger, M., et al. (2021). Binary central stars of planetary nebulae identified with Kepler/K2. MNRAS 506, 5223–5246. doi:10.1093/mnras/stab2045

CrossRef Full Text | Google Scholar

Jacoby, G. H., Ciardullo, R., Roth, M. M., Arnaboldi, M., and Weilbacher, P. M. (2024). Toward precision cosmology with improved planetary nebula luminosity function distances using VLT-MUSE. II. A test sample from archival data. ApJS 271, 40. doi:10.3847/1538-4365/ad2166

CrossRef Full Text | Google Scholar

Jones, D., and Boffin, H. M. J. (2017). Binary stars as the key to understanding planetary nebulae. Nature 1, 0117. doi:10.1038/s41550-017-0117

CrossRef Full Text | Google Scholar

Méndez, R. H., Thomas, D., Saglia, R. P., Maraston, C., Kudritzki, R. P., and Bender, R. (2005). Oxygen and neon abundances of planetary nebulae in the elliptical galaxy ngc 4697. ApJ 627, 767–781. doi:10.1086/430498

CrossRef Full Text | Google Scholar

Miller Bertolami, M. M. (2016). New models for the evolution of post-asymptotic giant branch stars and central stars of planetary nebulae. A&A 588, A25. doi:10.1051/0004-6361/201526577

CrossRef Full Text | Google Scholar

Reid, W. A., and Parker, Q. A. (2010). A new population of planetary nebulae discovered in the Large Magellanic Cloud - iii. The luminosity function. MNRAS 405, 1349–1374. doi:10.1111/j.1365-2966.2010.16635.x

CrossRef Full Text | Google Scholar

Rekola, R., Richer, M. G., McCall, M. L., Valtonen, M. J., Kotilainen, J. K., and Flynn, C. (2005). Distance to NGC 253 based on the planetary nebula luminosity function. MNRAS 361, 330–336. doi:10.1111/j.1365-2966.2005.09166.x

CrossRef Full Text | Google Scholar

Sarzi, M., Mamon, G. A., Cappellari, M., Emsellem, E., Bacon, R., Davies, R. L., et al. (2011). The planetary nebulae population in the central regions of M32: the SAURON view. MNRAS 415, 2832–2843. doi:10.1111/j.1365-2966.2011.18900.x

CrossRef Full Text | Google Scholar

Soker, N. (2006). Accreting white dwarfs among the planetary nebulae Most luminous in [O III] λ5007 emission. ApJ 640, 966–970. doi:10.1086/500291

CrossRef Full Text | Google Scholar

Souropanis, D., Chiotellis, A., Boumis, P., Jones, D., and Akras, S. (2023). Planetary nebulae hosting accreting white dwarfs: a possible solution for the mysterious cut-off of planetary nebula luminosity function? MNRAS 521, 1808–1816. doi:10.1093/mnras/stad564

CrossRef Full Text | Google Scholar

Valenzuela, L. M., Méndez, R. H., and Miller Bertolami, M. M. (2019). Revised simulations of the planetary nebulae luminosity function. ApJ 887, 65. doi:10.3847/1538-4357/ab4e96

CrossRef Full Text | Google Scholar

Valenzuela, L. M., Remus, R.-S., Miller Bertolami, M. M., and Méndez, R. H. (2025a). PICS: planetary nebulae in cosmological simulations – revelations of the planetary nebula luminosity function from realistic stellar populations. IAU Symposium 384, 69. doi:10.1017/S1743921323005458

CrossRef Full Text | Google Scholar

Valenzuela, L. M., Miller Bertolami, M. M., Remus, R.-S., and Méndez, R. H. (2025b). The PICS project. I. The impact of metallicity and helium abundance on the bright end of the planetary nebula luminosity function. A&A 699, A371. doi:10.1051/0004-6361/202553974

CrossRef Full Text | Google Scholar

Ventura, P., Dell’Agli, F., Schneider, R., Di Criscienzo, M., Rossi, C., La Franca, F., et al. (2014). Dust from asymptotic giant branch stars: relevant factors and modelling uncertainties. MNRAS 439, 977–989. doi:10.1093/mnras/stu028

CrossRef Full Text | Google Scholar

Walsh, J. R., Jacoby, G. H., Peletier, R. F., and Walton, N. A. (2012). The light element abundance distribution in ngc 5128 from planetary nebulae. A&A 544, A70. doi:10.1051/0004-6361/201118580

CrossRef Full Text | Google Scholar

Yao, P. Z., and Quataert, E. (2023). The origin of the consistent planetary nebula luminosity function bright-end cutoff. ApJ 957, 30. doi:10.3847/1538-4357/acfed9

CrossRef Full Text | Google Scholar