© The Authors 2025

1 Introduction

For more than 30 years, planetary nebulae (PNe) have been used as distance indicators in the cosmic distance ladder through the luminosity function of their [O III] λ5007 magnitude in extragalactic systems, referred to as the planetary nebula luminosity function (PNLF; Jacoby 1989). It was observed that the PNLF has a universal bright-end cutoff magnitude independent of the galaxy morphology, and it has been found to work well as a distance indicator (e.g., Ciardullo et al. 1989; Jacoby et al. 1990; Rekola et al. 2005). Despite wide usage, the underlying physics of the universality of the PNLF has remained a mystery, in particular due to the question how old elliptical galaxies would be able to host any PNe as bright as the ones observed at the magnitude cutoff (Ciardullo 2012). The argument in support of this ability has been that old stars at the main-sequence turnoff are too low in mass to produce such high emission in [O III]. The PNLF has been successfully used as a distance indicator for galaxies at distances of up to 10-15 Mpc, but its use has been limited to within 20 Mpc because of the limitations of narrowband filter imaging (Roth et al. 2021).

In the past 10 years, there have been two significant advances, with one on the modeling side and one on the observational side. With respect to modeling, Miller Bertolami (2016) computed updated post-asymptotic giant branch (AGB) tracks of potential PN central stars. These tracks were the first major update to the formerly computed stellar tracks from over 30 years ago (e.g., Schönberner 1983; Vassiliadis & Wood 1994; Blöcker 1995). The new tracks feature a much faster post-AGB evolution with timescales shorter by a factor of three to ten and with brighter post-AGB stars. This update led to several PN models being developed in recent years with a particular focus on improving the understanding of the universal bright end of the PNLF. Overall, the brighter post-AGB evolution has slightly diminished the large discrepancy between the observed bright PNe in old stellar populations and the models (e.g., Gesicki et al. 2018; Valenzuela et al. 2019). Further models have been developed with alternative attempts of explaining the bright end of the PNLF, including accreting white dwarfs (WDs; e.g., Soker 2006; Souropanis et al. 2023), blue straggler stars (e.g., Ciardullo et al. 2005), and WD mergers (e.g., Yao & Quataert 2023).

On the observational side, the integral field unit (IFU) measurements with the instrument MUSE have made it possible to obtain the full spectra for each individual pixel, thereby making the detection of PNe and the removal of contaminants possible without follow-up measurements in pencil beam observations. This has enabled the study of PNe at greater distances and complemented previously existing methods such as the PN Spectrograph (PN.S; Douglas et al. 2002; see also the extended PN.S early type galaxy survey, ePN.S; Pulsoni et al. 2018). In turn, the PN.S has the advantage of a significantly larger field of view (100 arcmin² versus 1 arcmin² for MUSE) and a higher spectral resolution, which has enabled detailed dynamical studies of galaxies using PNe (e.g., Aniyan et al. 2018, 2021, for face-on disks). Using MUSE, several works have analyzed extragalactic PN populations, obtaining better data for previously studied objects or investigating new systems (e.g., Kreckel et al. 2017; Spriggs et al. 2020; Roth et al. 2021; Scheuermann et al. 2022; Soemitro et al. 2023; Jacoby et al. 2024). Furthermore, the high spatial resolution of MUSE is expected to make it possible to reach distances of up to 40 Mpc, which is more than twice as far compared to previous limits (Roth et al. 2021).

To establish a theoretical foundation for the PNLF and the usefulness of its bright end as a distance indicator, numerous approaches have been taken throughout the years to explain the properties of the PNLF. Some authors have created models for the entire population of PNe based on mock stellar populations (e.g., Méndez et al. 1993; Méndez & Soffner 1997; Schönberner et al. 2007; Valenzuela et al. 2019), and others have employed analytic properties of the PNLF (e.g., Rodríguez-González et al. 2015). In contrast, the majority of studies have focused on the bright end specifically, which is the maximum brightnesses achievable by PNe of certain stellar populations (e.g., Dopita et al. 1992; Marigo et al. 2004; Gesicki et al. 2018; Souropanis et al. 2023; Yao & Quataert 2023). One of the findings from theory and observations has been that there should be a slight dependence of the PNLF bright end on metallicity, where the bright end becomes fainter for very metal-poor systems (e.g., Ciardullo & Jacoby 1992; Dopita et al. 1992; Ciardullo et al. 2002; Ciardullo 2010). However, the theoretical expectation that older stellar populations would produce PNLFs with a bright end that is significantly dimmer has not been observed (e.g., Marigo et al. 2004; Ciardullo 2012, stating that the PNLF bright end should already be more than 4 magnitudes dimmer for ages of 10 Gyr).

Not only has the bright end of the PNLF drawn attention, but the shape of the faint part of the PNLF has also been the subject of research. The general shape was introduced in its analytical form by Ciardullo et al. (1989) as an exponential cutoff function with a parameter describing the slope on the faint side based on theoretical derivations of Henize & Westerlund (1963) and Jacoby (1980). This slope has been predicted to encode clues about the star formation history (e.g., Ciardullo 2010; Longobardi et al. 2013; Hartke et al. 2017) and has been investigated in particular for the different substructures in the outskirts of M 31 (Bhattacharya et al. 2021). For this, Bhattacharya et al. (2021) also used the approach of fitting a superposition of two modes of the PNLF following Rodríguez-González et al. (2015). A further interesting feature that has been observed in the faint part of the PNLF is a dip for galaxies with recent star formation, that is a lower abundance of PNe at an intermediate magnitude (e.g., Jacoby & De Marco 2002; Reid & Parker 2010; Soemitro et al. 2023). Some models have even been able to find such a dip in the PNLF for younger stellar populations (e.g., Méndez et al. 2008; Valenzuela et al. 2019).

With the multitude of observations and modeling approaches, the continual lack of fundamental understanding concerning the PNLF suggests that some of the previously employed assumptions are insufficient to address the problem as a whole. Indeed there are two central assumptions that have been commonly made in PN models over the years that do not represent realistic stellar populations as they are found in actual galaxies. First, models have mostly only taken into account stellar evolution at metallicities around the solar metallicity, which has been done with respect to their main sequence (MS) and post-MS lifetimes, the initial-to-final mass relation (IFMR), and the post-AGB evolution (e.g., Gesicki et al. 2018; Valenzuela et al. 2019; Souropanis et al. 2023). Furthermore, even stellar evolution models employ many recipes based on measurements of the Sun and assume abundances based on the Milky Way (e.g., Miller Bertolami 2016). For example, the helium abundance has mostly been assumed to scale linearly with metallicity, though this has been put into question recently (e.g., Clontz et al. 2025). Second, the studies that have attempted to model the entire PN population have set up artificial mock stellar populations with the goal of mimicking selected real galaxies (e.g., Méndez & Soffner 1997; Valenzuela et al. 2019), but they again have relied on observations in the Milky Way, specifically the WD mass distribution. These models therefore fall short of accurately representing the large diversity of galaxies and star formation histories, which can strongly vary depending on intrinsic factors and the cosmological environment.

To overcome these limitations of current PN models, we introduce PICS (PNe In Cosmological Simulations), a modular framework for modeling PNe within any given stellar population that takes into account properties such as mass, age, metallicity, and the initial mass function (IMF). In the era of hydrodynamical cosmological simulations, a new approach is further possible. Realistic stellar populations with self-consistent properties in the cosmological context can be directly obtained for a large variety of galaxies from such simulations. With PICS, we present the very first method for bringing PNe into cosmological simulations and simultaneously improving on previous PN models by fully considering the metallicity and specific abundances beyond solar values. Finally, the modular nature of the framework allows for controlled testing of different parameters and individual models, which we also demonstrate in this paper.

In Sect. 2, we present the PICS model structure and the implementation of the fiducial model, including the metallicity dependence on the individual modules. By combining some of the modules, we then show how the metallicity and helium abundances affect the stellar evolution relevant for PNe in Sect. 3. In Sect. 4, we present grids of PNLFs resulting from different pairs of ages and metallicities and from varying the IFMR - these are the building blocks for any real PNLF. Finally, we summarize the results and conclude in Sect. 5.

2 Model

In the following, we describe the building blocks of the PICS model. Its target is to take a given single stellar population (SSP) parameterized by a set of properties and determine the contained PN population. The obtained PN properties in particular include the magnitude of the emitted [O III] line, as it is the primary intrinsic observable quantity of an extragalactic PN. This procedure will make it possible to apply PICS to both observed galaxy star formation histories discretized by SSPs and to the stellar particles extracted from cosmological simulations (Soemitro et al., in prep.; Valenzuela et al., in prep.). The end result is then a full population of PNe for an entire system like a galaxy or galaxy cluster.

The model is designed to be modular, meaning that the individual parts described in the following can be replaced with different underlying models, and the selection of models will also be continuously expanded in future versions. In this work, we present the fiducial model and some alternative relations for the IFMR and how they affect the PNLF. The source code and more information on using the model are made available online¹. A summary of the general model functionality is presented in the proceedings by Valenzuela et al. (2024). The correction for nebular metallicity and the inclusion of circumnebular extinction are presented here for the first time.

2.1 Single stellar population properties

An SSP represents the population of stars having formed from a gas cloud where all the member stars share the same underlying age and chemical composition. The distribution of stellar masses within the SSP is given by an IMF. In its simplest form, the input SSP properties need to include its total mass, an IMF, and its age. The age is needed to determine the initial mass range of stars now residing within the PN phase after having left the AGB track. The IMF is used to determine the fraction of stars that are within that initial mass range, and the total mass sets the absolute number of those stars. In addition to these properties, further quantities significantly influence the processes relevant for stellar evolution and PNe: in particular, we show in this work that the stellar metallicity is a vital parameter that affects the duration of stellar evolution on the MS and post-MS. In the future, individual abundances like helium or α-element abundances will also be added as quantities that are taken into account.

In this work, each SSP is parameterized with the following quantities:

• The IMF, for this work we only consider the Chabrier (2003) IMF.

• The total mass, M_tot, of the SSP.

• The age, t, measured as the time since the formation of the SSP stars.

• Metallicity, Z, which in this work we quantified as the initial metal mass fraction: Z = M_met,init/M_tot,init. As a point of reference, the solar metallicity lies between Z = 0.01 and 0.02 (e.g., Asplund et al. 2009; Caffau et al. 2011; von Steiger & Zurbuchen 2016; Asplund et al. 2021).

2.2 Stellar lifetime function

In the context of PN modeling, we define the stellar lifetime as the time from the birth of a star until it leaves the AGB phase and enters the post-AGB phase, during which a star can become a PN. It therefore roughly corresponds to the combined MS and post-MS lifetimes from other works (e.g., Renzini 1981; Renzini & Buzzoni 1986; Padovani & Matteucci 1993). The purpose of this module in PICS is to return the initial stellar mass of the stars that now are passing through the PN phase after the age t. Due to the short-lived nature of the PN phase (assumed as 30 000 yr here), the time spent after the departure from the AGB phase is negligible for the initial mass of the star. Then, disregarding binary interactions, we can assume that all PNe in a single stellar population are descendants of stars with the same initial mass.

For this work we took the lifetimes of the stellar evolution simulations from Miller Bertolami (2016), given in their Table 2. For each of their considered metallicities (Z = 0.0001, 0.001, 0.01, and 0.02) and initial masses (ranging from M_init = 0.8-4.0 M_⊙, depending on the metallicity), we added all the stated lifetimes: MS, red giant branch, core He-burning, early AGB, and thermally pulsing AGB phases as an M-type and carbon star. Taking these lifetimes and the initial masses from Miller Bertolami (2016), we fit a double power law for each of the four metallicities (Z = 0.0001, 0.001, 0.01, and 0.02) by fitting a power law to the lower masses up to and including 1.5 M_⊙ and a power law to the higher masses starting at 1.5 M_⊙ (i.e., with one value in common). Table 1 shows the fitting parameters for the free parameters of $$t/\mathrm{G}\mathrm{y}\mathrm{r}=\mathrm{a}\times ({\mathrm{M}}_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}/{\mathrm{M}}_{\odot }{)}^{\mathrm{k}}$$(1)

for the low and high masses at the four metallicities, which led to very good descriptions of the original data points, even in logarithmic space.

These fits only correspond to the four discrete metallicities from Miller Bertolami (2016). To obtain a continuous representation of the lifetime function across metallicities, we found the following relations to describe the change of the four double power-law parameters with metallicity²: $$a_{\mathrm{l}\mathrm{o}\mathrm{w}}=78.27\times Z^{0.649}+5.67,$$(2) $$a_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}=107.5\times Z^{0.880}+4.58,$$(3) $$k_{\mathrm{l}\mathrm{o}\mathrm{w}}=-0.943\times Z^{0.400}-3.38,$$(4) $$k_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}=-3.29-0.654\times (\log Z)-0.124\times (\log Z{)}^2.$$(5)

The lifetime is then determined as the maximum between t_low and t_high computed from the respective power law parameters obtained from the stated relations. The lifetime function is shown for the original and some selected further metallicities in Fig. 1, showing the continuous nature of the here derived expression. In Sect. 3, we present how the extrapolation of this relation to higher metallicities compares to detailed evolutionary models computed at such higher metallicities.

The obtained lifetime functions at the four original metallic-ities (Eqs. (2) to (5)) as well as at intermediate and one lower metallicity (Z = 0.00001) are shown in Fig. 1 as lines, colored by metallicity. The circular points indicate the values obtained from Table 2 of Miller Bertolami (2016) and lie directly on the respective thick solid functional lines. The thin lines, which are at interpolated and lower metallicities, lie between the thick lines, just as they should. The power law slope k_high becomes steeper with decreasing metallicity, such that the extrapolation to Z = 0.00001 leads to a very small change in slope toward higher masses. We therefore conclude that our determined analytical lifetime function is suitable for metallicities of at least Z = 0.00001-0.02.

While other PN studies have generally implicitly or explicitly assumed a fixed metallicity for their underlying mock stellar populations, oftentimes Z = 0.01 or Z = 0.02 (e.g., Gesicki et al. 2018; Souropanis et al. 2023; Valenzuela et al. 2019), Fig. 1 shows that the metallicity plays a very significant role for stars of any age. For a fixed stellar age, varying the metallicity from low values of Z = 0.0001 to Z = 0.02 results in a change of more than M_init = 0.2M. In contrast, large absolute changes of stellar age for t ≳ 5 Gyr only slightly affect the initial mass, which becomes relevant toward lower masses as the mass approaches the minimum viable mass for PNe to form, however. For old elliptical galaxies the metallicities tend to be high (e.g., Li et al. 2018), which means that especially for those types of galaxies it is important to take into account the actual underlying stellar metallicity. The implications of a higher metallicity are the following: at a fixed stellar age, more metal-rich stars now passing through the PN phase are more massive than metal-poor stars that are PNe because more metal-rich stars evolve more slowly at a fixed stellar mass. This leads to more massive central stars with higher temperatures and luminosities, thus giving rise to brighter PNe, which is a possible way of populating the bright side of the PNLF with PNe from older stellar populations.

Fig. 1 Lifetime function of stars until they reach the post-AGB phase as a function of their initial mass. Top: linear axes. Bottom: logarithmic axes. The logarithmic scale makes it easier to distinguish the individual lines and data points. The points are from the models by Miller Bertolami (2016), and the lines are the fit relations to those data points. The thick lines are at the same metallicities as the ones from the original models (Z = 0.0001, 0.001, 0.01, and 0.02), and the thin lines are at Z = 0.00001, 0.0003, 0.003, and 0.014, where the last three values are the logarithmic midpoints of the original metallicities.

2.3 Initial-to-final mass relation

Having obtained the initial mass of stars that are now PNe, the IFMR is next used to determine the final mass based on the initial mass and metallicity of the star, that is how massive the central star of the PN is. This relation has been mainly investigated through WD stars in star clusters in the solar neighborhood (e.g., Cummings et al. 2018), though it has also been estimated from composite stellar populations (e.g., El-Badry et al. 2018).

For the fiducial model, we create a metallicity-dependent IFMR based on the values obtained from the models of Miller Bertolami (2016), which includes the four metallicities Z = 0.0001, 0.001, 0.01, and 0.02, just as for the lifetime function (Sect. 2.2). To sample more data points between the simulated pairs of initial and final mass values at each metallicity, a forward-looking quadratic interpolation (Bhagavan et al. 2024) is used. Then, a 2D linear interpolation is applied between the four logarithmic metallicities to the more finely sampled values, which is the resulting metallicity-dependent IFMR. This IFMR can be applied in two extrapolation modes: for metal-licities outside of the range of Z = 0.0001-0.02, the values at the respective bounding metallicity can be used as a constant beyond that value, or the linear fit can be extrapolated beyond the original metallicity range. In Fig. 2, we show the resulting IFMR with extrapolated values at the same metallicities as shown for the lifetime function in Fig. 1, plus the higher metal-licity of Z = 0.05. The colored dots show the original data points from Table 2 of Miller Bertolami (2016), between which the interpolation took place.

The dependence on metallicity of the IFMR shows that more metal-rich stars at a fixed initial mass end up having lower final masses than metal-poor stars. As more massive stars are hotter and more luminous in the post-AGB phase and thus lead to brighter PNe, this counteracts the effect that metallicity has on the lifetime function, which leads to stars with higher initial masses entering the PN phase.

It is worth noting that the metallicity dependence of the IFMR of the models of Miller Bertolami (2016) is a consequence of the adopted winds. The metallicity dependence of cold winds is still poorly understood, in particular in the first red giant branch phase. The resulting theoretical IFMRs should therefore be taken as educated guesses. Most reliable IFMRs come from semi-empirical studies of star clusters in the solar neighborhood, which are only available for near-solar metallicities, however.

As the IFMR in general and its metallicity dependence are still poorly understood, it is important to better understand the dependence of the resulting PNe in a given SSP on the used IFMR. For this, we additionally consider four IFMRs derived from observed WDs in the Milky Way: Cummings et al. (2018) determined their IFMR from 80 WDs located in 13 star clusters, El-Badry et al. (2018) estimated the IFMR from Gaia measurements of over 1000 WDs within 100 pc fo the Sun, Cunningham et al. (2024) studied 40 WDs with spectroscopic measurements from Gaia data within 40 pc of the Sun, and Hollands et al. (2024) studied 90 WD binaries from Gaia data with distances of up to 250 pc. All four of these IFMRs correspond to nearsolar metallicities. We also show IFMRs from the four studies in Fig. 2, where we use the adopted MIST-based IFMR of Cummings et al. (2018) and linearly extrapolate the IFMRs down to lower initial masses, as Cunningham et al. (2024) and Hollands et al. (2024) only fit their IFMRs down to 1 M_⊙, and Cummings et al. (2018) to 0.83 M_⊙. The uncertainties of these IFMRs are on the order of σM_final ≈ 0.025 M_⊙, for which we also show the uncertainty region of Cunningham et al. (2024). Further, the systematic uncertainties are on the order of 0.1 M_⊙, as seen by using different underlying stellar evolution model is employed (Cummings et al. 2018). Only the low-mass IFMR from El-Badry et al. (2018) produces significantly lower final masses than the other relations. According to Cunningham et al. (2024), the reason for this is that El-Badry et al. (2018) use a different method for obtaining the minimum WD mass that can be created through single stellar evolution. While these relations are all monotonic, we note that it is still debated whether this is the case or if there is a dip around M_init ≈ 1.5-2M_⊙ (e.g., Marigo et al. 2020) due to the large uncertainties.

Overall, there is a general agreement on what the IFMR looks like within the measurement uncertainties and their large systematics, where the theoretical metallicity-dependent IFMR from Miller Bertolami (2016) predicts overall lower final masses at solar metallicity than the observationally based IFMRs, especially at final masses of M_init ≈ 1.5-3 M_⊙ around the potential dip. Still, the small differences in final mass can lead to significant changes in the resulting PN properties, especially moving toward final masses below 0.6M. In Sect. 4, we investigate how the different IFMRs impact the resulting PN population with respect to their PNLF.

Fig. 2 Initial-to-final mass relation of stars. The points are from the models by Miller Bertolami (2016), and the lines in the corresponding colors are the fit relations to those data points. The thick lines are at the same metallicities as the ones from the original models (Z = 0.0001, 0.001, 0.01, and 0.02) and the thin lines are at Z = 0.00001, 0.0003, 0.003, 0.014, 0.05, where the middle three values are the logarithmic midpoints of the original metallicities. Additionally, four observationally based IFMRs from nearby WD measurements at solar metallicity are shown from Cummings et al. (2018), El-Badry et al. (2018), Cunningham et al. (2024), and Hollands et al. (2024), as well as the 1σ uncertainty band of Cunningham et al. (2024).

Fig. 3 Post-AGB tracks in the Hertzsprung-Russell diagram of stars at two final masses Mfinal = 0.55 M⊙ and 0.58 M⊙ for a range of metallicities. The tracks are interpolated between the data tables from Miller Bertolami (2016). The metallicity colorbar corresponds to the same coloring as found in Figs. 1 and 2. The thick lines are at the same metallicities as the ones from the original models (Z = 0.0001, 0.001, 0.01, and 0.02) and the thin lines are at Z = 0.00001, 0.0003, 0.003, 0.014, 0.05, where the middle three values are the logarithmic midpoints of the original metallicities.

2.4 Post-AGB stellar evolution tracks

With the final mass of the star in question, the evolving properties of the central star can be obtained from post-AGB stellar evolution tracks. In particular, the luminosity L and effective temperature T_eff can be directly obtained with the tracks in the Hertzsprung-Russell diagram, given the final mass and an identifying time. For PICS, we assume the lifetime of a PN can be up to 30 000 yr (Valenzuela et al. 2019). For a given central star, the time since leaving the post-AGB phase is therefore randomly drawn from a uniform distribution between 0 and 30 000 yr, as the timescale is negligible compared to the timescale of the MS and post-MS.

In this work, we use the post-AGB tracks of the H-burning stars from Miller Bertolami (2016) and interpolate them to obtain the central star luminosity and effective temperature as a function of metallicity, final mass, and the post-AGB age (see Appendix A for details on the interpolation). Following the approach of Valenzuela et al. (2019), we take the post-AGB age to be the time since the temperature reached a value of T_eff = 25 000 K (log T_eff/K ≈ 4.40). The interpolated post-AGB tracks at two final masses of M_final = 0.55 and 0.58 M_⊙ (of which only the former mass was explicitly computed by Miller Bertolami 2016) for a range of different metallicities can be seen in Fig. 3. Clearly, the interpolations successfully produce intermedial values, in particular between the originally simulated metallicities. For metallicities outside of the range between Z = 0.0001 and 0.02, we simply use the respective tracks at Z = 0.0001 or Z = 0.02.

Additionally, it can here be seen that more metal-poor stars at a fixed final mass have a higher luminosity and can reach higher temperatures than metal-rich stars, as was also shown by Miller Bertolami (2016). Like the metallicity-dependent IFMR from Miller Bertolami (2016) presented in Sect. 2.3, this counteracts the effect the lifetime function has that one would expect more metal-rich stars to host brighter PNe, as well as the lower number of oxygen atoms in the nebula, which is also expected to limit the brightness of the [O III] line.. The dependence of the maximum luminosity and effective temperature of a star on its metallicity and age is therefore not trivial.

Fig. 4 Post-AGB tracks in the Hertzsprung-Russell diagram of stars at metallicities Z = 0.001, 0.01, and 0.02 (left, middle, and right, respectively) colored by the λ5007 magnitude. The tracks are shown for stars at different final masses Mfinal = 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.60, 0.65, 0.70, 0.80, 0.90, and 1.00M. For each final mass, stars at 600 evenly spaced post-AGB ages between 0 and 30 000 yr are plotted. The white lines indicate the tracks for the selected central star masses.

2.5 Planetary nebula model

The central star properties can finally be used to infer the PN emissions. Particularly, the [O III] intensity is of great importance for studies of extragalactic PNe. Additionally, properties such as the Ηβ intensity may be of further relevance and have been studied for extragalactic PN populations as well (e.g., Reid & Parker 2010).

In this work, we use the PN model from Valenzuela et al. (2019). To our knowledge it is the only current PN model that not only addresses the brightest PNe at the bright end of the PNLF, but can also be used to produce dim PNe that are consistent with observations, while being applicable to large numbers of central stars without exceedingly high computational cost. The model computes the maximum total intensity of Hβ emission assuming optical thickness, I(Hβ), based on the luminosity of the central star luminosity and its blackbody radiation given by the effective temperature. Next, the ratios between I(Hβ) and the intensity of the [O III] emission line I(5007) that have been observed for PNe in the Milky Way and the Large Magellanic Cloud (LMC) are used to empirically obtain I(5007) for a given PN. This ratio is drawn from a Gaussian distribution and is modified by some physically motivated prescriptions based on stellar mass and temperature. We give an overview of this prescription in Appendix B.1.

Finally, an absorbing factor μ is applied to the intensities to account for ionizing stellar photons that are not absorbed by the nebula. The absorbing factor is equivalent to the fraction of ionizing photons that are actually absorbed and is simply multiplied with the previously obtained intensity to get the observed value of I(5007). The absorbing factor for a given star is drawn randomly because of possible asymmetries in the nebula that can lead to a leakage of ionizing photons (for more details see Méndez & Soffner 1997; Valenzuela et al. 2019). Additionally, its maximally allowed value is decreased over the lifetime of the PN, which accounts for the nebula dissipating over time. The prescription for the absorbing factor is described in Appendix B.2. One aspect that is not included in the model of Valenzuela et al. (2019) is the internal dust extinction, which becomes more important for younger stars (e.g., Jacoby & Ciardullo 2025), which means in general for higher stellar masses. We add this effect as a further step in the PICS framework in Sect. 2.6.

In Fig. 4, we show the post-AGB tracks in the Hertzsprung-Russell diagram of stars with certain final central star masses at metallicities of Z = 0.001, 0.01, and 0.02. The data points show individual PNe at evenly spaced ages between 0 and 30 000 yr starting from the given track reaching T_eff = 25 000 K. The PNe are colored according to their [O III] λ5007 absolute magnitude, M(5007) = -2.5 log I(5007) - 13.74 with I(5007) being the intensity seen 10 pc away, as defined by Jacoby (1989). According to the model, the brightest PNe for a given central star mass are located on the rising temperature track, after which the intensity of the [O III] line decreases again due to the dissipation of the nebula over time and the eventual drop of the star’s luminosity. Here it also becomes apparent how quickly more massive stars (M_final ≳ 0.6 M_⊙) evolve along the tracks compared to less massive stars, for example at Z = 0.001: two data points are separated by 50 yr, which means that after only a little more than 100 yr, a central star of mass 0.8 M_⊙ reaches its peak temperature. As a result, even though more massive stars have the potential of developing brighter PNe, the likelihood of observing such a PN in its bright stage is much lower for more massive stars due to their fast evolution.

The PN model by Valenzuela et al. (2019) was originally developed for use with the post-AGB tracks from Miller Bertolami (2016) at a metallicity of Z = 0.01 and the [O III] to Hβ intensity ratio was empirically based on Milky Way and LMC PNe, which span a range of metallicities. As there is no obvious step in the model’s recipes where it would be trivial to add a direct dependence on metallicity, in this work we apply the PN model to central stars of all metallicities without modification. While the model therefore statistically works well enough for composite populations of PNe that roughly have the same metallicity distribution as the Milky Way and LMC PNe, the intrinsic dependence on metallicity is lost, where lower abundances of oxygen should generally lead to less [O III] emission.

For this reason, only applying the PN model by Valenzuela et al. (2019) to especially metal-poor PNe will result in overluminous PNe. Multiple studies have investigated the dependence of the maximum [O III] flux on the metallicity. This relation is not straightforward as [O III] is also one of the primary coolants of the nebula (e.g., Jacoby 1989). Using nebular models, Jacoby (1989) found that small changes in the metallicity at solar-like abundances resulted in variations of the [O III] flux by the square root. Later models by Dopita et al. (1992) and Schönberner et al. (2010) showed a quadratic form in logarithmic space around solar abundances for the peak luminosity of a given PN (see Fig. 27 from Schönberner et al. 2010). In this work we consider all three possible relations as direct corrections to the obtained [O III] fluxes, extrapolating them beyond solar abundances to both low and high metallicities. As the nebular model by Valenzuela et al. (2019) was not based on a particular metal-licity to which a simple nebular metallicity correction can be applied, we needed to estimate a base metallicity Z_base of the model from which it is possible to use the abundance differences to determine the correction required for the [O III] flux. Because Milky Way and LMC PNe were used for creating the empirical model by Valenzuela et al. (2019), the value has to be equal to or lower than the maximum metallicities of both galaxies. The overall more metal-poor LMC has an upper metal-licity of around Z = 0.008 (Harris & Zaritsky 2009), which is consistent with the upper oxygen abundance found for its PNe of 12 + log(O/H) = 8.5-8.6 (Dopita & Meatheringham 1991; Richer 1993). At the same time the base metallicity value may not be too low as the observed nebulae have to have a sufficient amount of oxygen available to become PNe in the first place. We therefore estimate this base metallicity as Z_base = 0.007, a value near the upper metallicity value of the LMC and also roughly half of the solar metallicity. We decided on this value because it is within the reasonable metallicity range of both the Milky Way and the LMC. However, it is important to note that this is not a physically derived value due to the nature of the empirical model being based on a distribution of metallicities.

For the square root relation of Jacoby (1989), we obtained a nebular metallicity correction factor of $$p_{5007,\mathrm{J}89}(Z)={\left(\frac{Z}{Z_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}}\right)}^{0.5}.$$(6)

For Dopita et al. (1992), we assumed a linear conversion between the metallicity and oxygen abundance [O/H], so with our choice of the base metallicity, we obtained a correction factor of $$p_{5007,\mathrm{D}92}={10}^{\log F_{5007,\mathrm{D}92}(\log Z/Z_{\odot })-\log F_{5007,\mathrm{D}92}(\log Z_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}/Z_{\odot })},$$(7) $$\log F_{5007,\mathrm{D}92}([Z])=-0.8791+0.1459[Z]-0.3013[Z{]}^2,$$(8)

where the latter equation is taken from their eq. (4.1) with [Z] = log Z/Z, and we took Z = 0.0134 from Asplund et al. (2009), which is also the solar metallicity that Schönberner et al. (2010) assumed in their models. Finally, for Schönberner et al. (2010) we fit two quadratic functions to their data points for their two final masses and used the mean of each of the fit parameters to obtain the following correction factor: $$p_{5007,\mathrm{S}10}={10}^{\log F_{5007,\mathrm{S}10}(\log Z/Z_{\odot })-\log F_{5007,\mathrm{S}10}(\log Z_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}/Z_{\odot })},$$(9) $$\log F_{5007,\mathrm{S}10}([Z])=-0.24([Z]-0.78{)}^2.$$(10)

We show the direct comparison of the three correction factors as a function of metallicity in Appendix C. As the correction factor only depends on metallicity, the values of M(5007) in Fig. 4 are simply offset by a fixed value per panel when applying one of the three nebular metallicity corrections.

While the nebular metallicity relations from the three works were originally constructed to describe how the fluxes of the brightest PNe vary with metallicity, we have decided to use these relations for correcting the [O III] fluxes of all PNe irrespective of their properties. We also note the caveat of applying the nebular metallicity correction to the resulting fluxes from the model by Valenzuela et al. (2019): as the model was tuned to PNe with a range of metallicities, having to pick a single base metallicity from which the correction is determined leads to an excessive brightness of PNe on the slightly metal-rich side and an excessively strong dimming on the slightly metal-poor side. As the disadvantages of not applying such a correction outweigh this caveat, we picked as the fiducial model the nebular metallic-ity correction by Dopita et al. (1992), which does not boost the [O III] flux as much as the other two relations (see Fig. C.1).

Finally, we comment on the fact that for our current fiducial nebular model, we only take into account the general metallic-ity mass fraction Z and not any specific abundances. On the one hand, this has the big advantage of keeping the model simple, especially in light of the large parameter space already spanned through the stellar evolution models, which has also been the focus for this paper. With the aim of applying PICS to cosmological simulations, most chemical evolution models are not yet able to find good agreement with observed abundances, whereas the total metallicity mass fractions are in general consistent with observations. On the other hand, the simple linear scaling assumptions of different abundances like helium, oxygen, and iron with Z are certainly not entirely correct and may lead to systematics and biases in the modeled PN emission lines. For instance, accounting for alpha elements independently of iron-peak elements may be relevant when treating galaxies with different star formation histories like typical elliptical versus disk galaxies. Similarly, increased helium abundances could affect the number of ionizing photons available to ionize hydrogen (e.g., Osterbrock & Ferland 2006). As the focus of this paper is on the advances of stellar evolution modeling in the context of PNe, treating the nebular physics with a more detailed model will be the subject of future work, further improving the PICS model.

2.6 Circumnebular extinction

During the AGB phase stars create and eject dust through their stellar winds, which may lead to extinction of the PN emission. As more massive stars evolve more quickly along the post-AGB tracks and thus reach high temperatures to ionize the nebula shortly after the AGB phase, these PNe would suffer a greater extinction than PNe of low-mass stars when they reach their brightest [O III] magnitude. Ciardullo & Jacoby (1999) suggested that this effect could explain how the bright end of the PNLF is even universal for PNe in younger stellar populations, where the central stars of the PNe would be more massive. For the bulge of M 31, Davis et al. (2018) found that a significant fraction of the brightest PNe are strongly extincted. Most recently, Jacoby & Ciardullo (2025) presented the observed relation between central star mass and circumnebular extinction measured for the brightest PNe in the LMC and M 31, which they fit through linear relations. Their work extends on Ciardullo & Jacoby (1999) and presents strong indications that this effect indeed could explain how younger galaxies have the same universal bright-end cutoff of the PNLF as the circumneb-ular extinction could bring the overly bright PNe down to the cutoff.

As PICS has the objective to produce PN properties for direct comparison with observations, accounting for circumneb-ular extinction is an important last step for obtaining the actually observed [O III] flux. While it is in part possible to determine the circumnebular extinction of PNe in the Local Group, this is impossible for most extragalactic systems. The relations for extinction as a function of final mass provided by Jacoby & Ciardullo (2025) are only valid for the PNe in the brightest magnitude. As the dust becomes more diffuse over the lifetime of a PN, the extinction is expected to decrease. Under the assumption that the evolution of dust extinction is the same for all PNe, we derived the extinction relation as a function of post-AGB age that recovers the relation by Jacoby & Ciardullo (2025) for the PNe at metallicity Z = 0.01 as they reach their peak brightness. We chose the metallicity of Z = 0.01 as a reasonable value between the expected metallicities of the brightest PNe in M 31 and the LMC. We note that using slightly lower or higher metallicities do not affect the resulting extinctions significantly.

We first determined the post-AGB age corresponding to the maximal brightness as a function of final mass at Z = 0.01 and fit the following broken linear relation to it (see Appendix D for the details of the derivation): $$M_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l},\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}=max\left\{\begin{array}{r}1.014-0.1704\log (t_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}-\mathrm{A}\mathrm{G}\mathrm{B}}/\mathrm{y}\mathrm{r})\\ 0.729-0.0479\log (t_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}-\mathrm{A}\mathrm{G}\mathrm{B}}/\mathrm{y}\mathrm{r})\end{array}\right\}.$$(11)

In this work, we use the combined orthogonal regression fit to the bright LMC and M 31 PNe of Jacoby & Ciardullo (2025, Combined OR in their Table 2), which takes the corresponding brightest final mass as its argument (independent of what the actual final mass is): $$c_{\mathrm{H}\beta }=7.05M_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l},\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}/M_{\odot }-3.97,$$(12)

where the extinction in Hβ, cHβ, is used to obtain the [O III] extinction, A₅₀₀₇, through the Cardelli et al. (1989) extinction law with A_V = R_VE(B - V), R_V = 3.1, and E(B - V) = c_Hβ/1.47 (Ciardullo & Jacoby 1999): $$A_{5007}=\frac{c_{\mathrm{H}\beta }}{1.47}(R_{V}a(x)+b(x)).$$(13)

Here, a(x) and b(x) are defined for the optical/near-infrared range in Eqs. (3a) and (3b) from Cardelli et al. (1989) with x = 1 μm/5007 A. In this way, the extinction only depends on the post-AGB age of a PN. While in this work we only apply the one circumnebular extinction relation from Jacoby & Ciardullo (2025), we will investigate the effects of variations of the relation in a forthcoming paper (Valenzuela et al., in prep.).

2.7 Number of planetary nebulae

The previous sections dealt with the modeling of a single PN based on the IMF, age and metallicity of the SSP. How many PN this SSP produces is additionally dependent on the total mass. For this the IMF is integrated between the initial mass obtained from the SSP age (Sect. 2.2) and the initial mass obtained from the SSP age minus 30 000 yr because those are the stars that are still in the PN phase. The integrated value is then multiplied with the total mass to obtain the expected total number of PNe in that SSP. The actual integer number is then determined by rounding down or up randomly based on the decimal places, for example an expected value of 10.2 would lead to 10 PNe being selected with an 80% chance and 11 PNe with a 20% chance.

The number of PNe is obtained after the initial stellar mass has been found from the lifetime function. The final mass and the PN properties are then determined individually for each of those PNe since values such as the post-AGB age or the absorbing factor can be different for PNe even within the same SSP.

3 The effect of metallicity and helium abundances

3.1 Metallicity

In the previous section, it can be seen that multiple aspects of the emergence of PNe depend on the metallicity: For a given initial mass, more metal-rich stars remain on the MS and post-MS for a longer period of time (Fig. 1), and they lose more mass assuming the IFMR from Miller Bertolami (2016), leading to smaller final masses after the AGB phase (Fig. 2). For the post-AGB tracks at a given final mass, more metal-rich stars have lower luminosities and reach lower maximum temperatures before entering the cooling track (Fig. 3). Finally, while not depicted in this work, stars of a given final mass move along the Hertzsprung-Russell diagram at roughly the same speed, that is they heat up and cool down on a similar timescale (Fig. 8 of Miller Bertolami 2016).

It is not directly clear from the individual relations themselves how the metallicity affects the central stars of the PNe when combining the lifetime function and the IFMR. To shed more light on this, we show how the final stellar mass depends on the age with varying metallicity in the left panel of Fig. 5. For a large age range, most of the metallicities lead to similar final masses: at ages of 2-7 Gyr all but Z = 0.001 lie closely together, and at ages of 5-11 Gyr, all but Z = 0.05 lie closely together. Based on this model, this implies that there is a nearly universal final mass relation with age independent of metallic-ity for these ages, with variations of around 0.01-0.02 M_⊙ at a given age. For lower and higher ages, however, the final masses begin to deviate more from one another. For stellar populations less than 2Gyr old, stars becoming PNe that are more metal rich tend to have smaller final masses. For ages above 8-10 Gyr and with Z > 0.001, the trend emerges under the current model assumptions that more metal-rich stellar populations produce PNe with larger final masses due to their longer lifetimes compared to their lower-mass counterparts. While these trends are not monotonic in the shown metallicity range, old stellar populations that are especially metal rich above solar metallicity produce significantly more massive central stars of PNe than at lower metallicities, still reaching final masses of 0.55 M_⊙ after 9 Gyr.

Based on these relations, we can now highlight the importance of taking the metallicity into account for especially metal rich systems: if one were to follow the approach of some previous studies of assuming a fixed roughly solar-like metallicity of Z = 0.01 (e.g., Gesicki et al. 2018; Valenzuela et al. 2019), then for ages above 5 Gyr the final masses at that fixed metallicity would be systematically lower than if taking the correct higher metallicity into account (see the left panel of Fig. 5). The information can also be presented in the IFMR plot from Sect. 2.3 by marking specific times along the IFMRs at each metallicity. The right panel of Fig. 5 shows the IFMR for a more limited range of initial and final masses than shown in Fig. 2 where the points along the IFMRs are marked at equally spaced lifetimes between 2 and 12 Gyr. Especially at higher metallicities it can be seen that very old stellar populations with t > 6 Gyr have stars of higher initial and final masses reach the PN phase, resulting in old and metal-rich elliptical galaxies having more massive central stars of PNe than the old stellar population ofa more metal-poor spiral galaxy, for example.

One caveat of the above conclusions based on the highest shown metallicity Z = 0.05 is that this metallicity lies outside the range originally modeled by Miller Bertolami (2016). Since those models only reached a metallicity of Z = 0.02, it cannot simply be assumed that the extrapolation of the lifetime function to higher metallicities is applicable using the fit analytic relation (Sect. 2.2). One difficulty in this matter is that stellar evolution and in particular the AGB phase are still poorly understood for super-solar metallicity stars. This is largely due to such metal-rich stars being mostly observed only in the Galactic bulge or in M 31 for massive stars (e.g., Pietrinferni et al. 2013; McDonald et al. 2022). Some studies have run models for super-solar metallicity stars, including the PARSEC (Bressan et al. 2012) and BaSTI (Pietrinferni et al. 2013) stellar evolution models, additionally Althaus et al. (2009) investigated potential high-metallicity progenitors of WDs, and Nanni et al. (2014) ran super-solar metallicity models to study the dust production in pulsing AGB stars. As stated by Pietrinferni et al. (2013), there are differences between their BaSTI and the PARSEC models in part because of different initial helium abundance assumed in the models.

Fig. 5 Final masses as a function of stellar lifetime (left) and of initial mass (right) for selected metallicities (Z = 0.0001, 0.001, 0.01, 0.02, and 0.05) based on the models of Miller Bertolami (2016). The final mass y-axis range is reduced to the primarily relevant values in the regime of older stellar populations of t > 2 Gyr. For the initial-to-final mass relation, the masses reaching the post-AGB phase after t = 2, 4, 6, 8, 10, and 12 Gyr are marked with circles.

3.2 Helium abundance

Not only is the initial helium abundance uncertain, its scaling with the overall metallicity has been primarily measured at lower metallicities in the Milky Way and in H II regions. Studies have mostly quantified a slope ΔY/ΔZ as the cosmic production rate of helium, which has been estimated to be between one and three (e.g., Jimenez et al. 2003; Izotov & Thuan 2004; Casagrande et al. 2007). However, not all authors assume a linear slope, especially at the lowest metallicities (e.g. Keszthelyi et al. 2024). While Miller Bertolami (2016) assumed a slope of two for the initial abundances of the stellar evolution models, it is unclear how exactly the helium-to-metallicity relation continues beyond the upper metallicity of Z = 0.02. Casagrande et al. (2007) even suggested that the slope becomes smaller toward higher metallicities, though the large uncertainties of the helium abundances and the metallicities make it difficult to draw definite conclusions on the true relation.

For these reasons, we ran two additional models following the procedure by Miller Bertolami (2016) with a super-solar metallicity of Z = 0.05. The first assumes the same linear relation between the helium abundance Y and metallicity Z of Y = 0.245 + 2Z (i.e., Y = 0.345, linear He), whereas the second assumes that a saturation of helium is reached at the linearly determined value at Z = 0.02 (i.e., Y = 0.285, cutoff He). The lifetime functions of the MS and post-MS for stars with these initial abundances as a function of the initial stellar mass are shown in Fig. 6 as the blue (linear He) and red (cutoff He) data points. Additionally, the values from the analytic relation obtained in Sect. 2.2 are shown as lines for the original four metallicities as well as for the extrapolated metallicity Z = 0.05. We find that the blue linear helium data points at Z = 0.05 lie very close to the relation at the lower metallicity Z = 0.02 (light-orange line), whereas the red cutoff helium data points at Z = 0.05 also correspond to the extrapolated analytical function at Z = 0.05 (yellow line).

The newly run models show that the lifetime not only depends on metallicity, but is also heavily influenced by the helium abundance at Z = 0.05. Another finding is that the overall trend of the lifetimes increasing with metallicity is not valid when extrapolating the linear relation between Y and Z beyond Z = 0.02. This can be qualitatively explained by the following: in general, an increase of Z leads to higher bound-free opacities in the star, which lowers the radiative energy transport and overall luminosity. More metal-rich stars thus spend more time on the MS and post-MS. This is the case at lower metallicities, where the hydrogen abundance only minimally increases with higher Z. However, at higher metallicities, a relative increase of Z results in a significantly greater change of the helium abundance following the linear relation of Y = 0.245 +2Z. As the hydrogen abundance is given by X = 1 - Y - Z, the higher value of Y means that the total amount of available hydrogen is decreased, resulting in less fuel for the MS and thus shorter MS lifetimes. Moreover, a higher helium abundance implies a higher mean molecular weight of the stellar material. Due to the high sensitivity of luminosity on the molecular weight (L ∝ μ⁴; Kippenhahn et al. 2013), this leads to a greater luminosity and faster fuel consumption. In summary, when extrapolating the linear relation between Y and Z , the stellar lifetime increases with metallicity until the rising helium abundance counteracts that effect, leading to similar lifetimes at Z = 0.05 and Z = 0.02. In contrast, letting the helium abundance saturate and stay constant for values of Z ≥ 0.02 means that the lifetimes continue to increase with metallicity.

For the saturated (cutoff) helium abundance, it is interesting that the simulated MS and post-MS lifetimes are overall very similar to the analytically computed values at Z = 0.05 for stars with initial masses M_init ≳ 1.5 M_⊙, but slightly higher at lower masses. This indicates that the extrapolation based on the lower metallicities (i.e., the analytic function that was fit to the original data points) corresponds to a slightly higher helium abundance than the saturated helium simulation. A thorough analysis of how the helium abundances affect the stellar lifetimes and the resulting PN populations is beyond the scope of this paper, however. It will be the subject ofa future study.

Due to the poor observational constraints with respect to the Y -Z relation, we conclude that the analytically derived lifetime function from Sect. 2.2 at super-solar metallicities is still consistent with the stellar evolution simulations within the uncertainties of the observed helium abundances. Moreover, there is recent observational evidence that the helium abundance saturates beyond a certain metallicity, as measured for Omega Centauri (Clontz et al. 2025). The first evidence supporting the usage of this function even at high metallicities has been presented by Valenzuela et al. (2024), where the fiducial model of PICS was applied to several galaxies from a cosmological simulation, spanning a wide range of star formation histories, ages, and metallicities. The thereby obtained PN populations and their PNLF are consistent with observations and give a robust view of PNe, as will be fully shown in detail in the next paper of this series (Valenzuela et al., in prep.). In a future study, we will further investigate how different assumptions regarding the helium abundances affect the PN populations of galaxies with different star formation histories in detail. Comparing the observed PN populations with modeled PNe based on simulated galaxies may even provide an independent way of probing the production rate of helium at high metallicities. In the following sections we apply the analytic function from Sect. 2.2 for all metallicities.

Fig. 6 Lifetime function of stars until they reach the post-AGB phase as a function of their initial mass. The lines are the fit relations to the models of Miller Bertolami (2016) with original metallicities of Z = 0.0001, 0.001, 0.01, and 0.02, as in Fig. 1. Here, the extrapolation of these fits to Z = 0.05 is also shown. The points are the results of newly run models following the procedure of Miller Bertolami (2016) for Z = 0.05, and the lines are the fit relations to those data points. The linear helium abundance is at Y = 0.345 and the cutoff helium abundance at Y = 0.285. The thick lines are at the same metallicities as the ones from the original models (Z = 0.0001, 0.001, 0.01, and 0.02) and the thin lines are at Z = 0.00001, 0.0003, 0.003, and 0.014, where the last three values are the logarithmic midpoints of the original metallicities.

4 Single stellar population PNLFs

In Sect. 3 we have highlighted the importance of the metallicity and helium abundance when considering some of the individual blocks required to model the PNe within a stellar population. We now showcase how the SSP age and metallicity affect the resulting PNe by exploring the parameter space for the fiducial PICS model. For this, we paired ten SSP ages between t_age = 0.25 and 13 Gyr with six metallicities between Z = 0.0001 and 0.08 to run PICS on and obtain corresponding PNLFs. To overcome the stochastic nature of the PN model itself with respect to the postAGB age and the absorbing factor that determines the opacity of the nebula, we ran PICS on 10⁶ stars for each pair of age and metallicity, with equally spaced post-AGB ages between 0 and 30 000 yr.

4.1 PNLFs with nebular metallicity correction

As a first step, Fig. 7 shows the PNLFs without accounting for circumnebular extinction for a grid of SSPs with different ages (constant age per row) and metallicities (constant metal-licity per column), that is the M(5007) histograms for the 10⁶ PNe. The colors represent the approaches to correcting for the nebular metallicity of Jacoby (1989), Dopita et al. (1992), and Schönberner et al. (2010) that were presented in Sect. 2.5, while the black PNLFs were created from PNe without applying such a correction and functions as a reference for what the corrections are doing. Statistically, the PNLFs in a given column are always shifted by the same amount relative to each other due to the way the metallicity corrections are implemented as a function solely of metallicity. For each panel, the final mass of the central stars is stated in the bottom right corner.

The effect of accounting for nebular metallicity is immediately obvious when comparing the left to the right columns: at low metallicities, the [O III] brightness drops considerably (colored lines), removing any PNe at Z = 0.0001 from bright magnitudes of M(5007) < -4 compared to the metallicity-independent nebular model (black line), while at high metallicities the correction leads to overall brighter PNe. The corrections are close to zero at Z = 0.01, which is due to the chosen base metallic-ity of Z_base = 0.007 at which the corrections are defined to not change the [O III] brightness. Among the three metallicity correction relations, the simple square root law by Jacoby (1989) always leads to the brightest PNe (yellow). At high metallicities, the relation by Dopita et al. (1992) enhances the PN brightness the least, as is expected from Appendix C, and again reaches a correction factor of 1 at around Z = 0.08, such that the purple and black lines are nearly the same in the last column of Fig. 7.

The first finding is that not all combinations of age and metallicity contribute to the bright end of the PNLF at around M(5007) = -4.5 (e.g., Ciardullo 2012), marked by the vertical gray dashed lines in Fig. 7. Old stellar populations with ages of t_age ≳ 8-10 Gyr generally produce fainter PNe. We note, however, that when nebular metallicity effects are not taken into account, the model predicts much brighter PNe at Z = 0.0001 and Z > 0.04. The final masses of the PNe in the stellar populations not reaching the bright end for the most part lie around M_final ~ 0.51-0.55 M_⊙. This range lies well below the expected maximum final masses of older observed stellar populations from previous studies needed to explain the brightest PNe (e.g., Méndez 2017; Valenzuela et al. 2019).

For intermediate-age stellar populations with t_age ~ 2-8 Gyr, the bright end of the PNLF is well populated with PNe. Brighter nebulae are found for lower ages. When accounting for the nebular metallicity (colored lines), this is only true for intermediate to high metallicities of Z ≳ 0.001. Even at Z = 0.001, the older intermediate-age stellar populations (t_age ≳ 6-8 Gyr) do not reach the bright end. Interestingly, there is a maximum of the non-metallicity-corrected PNLF between M(5007) = -4.5 and −4.0 for a number of the age-metallicity combinations in this intermediate-age regime, with a slow decline of the PNLF toward fainter PNe. This maximum is simply shifted to other magnitudes for the metallicity corrections. All of the SSP PNLFs here feature the recognizable typical shape of observed PNLFs, just with shifts toward brighter or fainter magnitudes. For the intermediate ages, the final masses lie between around 0.55 and 0.59 M_⊙. This approximately lies in the expected maximum central star final mass range for intermediate-to-old stellar populations as found by Valenzuela et al. (2019) based on modeling PNe in mock stellar populations.

The younger stellar populations at t_age < 0.5 Gyr only feature the typical PNLF shape at the higher metallicities when not accounting for nebular metallicity (black lines). In contrast, lower-metallicities lead to much brighter PNe existing even well beyond M(5007) = -5.0. At high metallicities, the nebular metallicity corrections of Jacoby (1989) and Schönberner et al. (2010) also lead to such bright PNe. The absolute numbers are overall much lower than for older stellar populations, which is a result of the much faster evolution along the post-AGB tracks for these more massive central stars (M_final ≳ 0.6). The stars reside on the heating track for only a short period of time (see the widely spaced data points in Fig. 4 for the more massive central stars) before cooling down and dropping in luminosity, thus making bright nebulae around those massive stars rather rare.

Fig. 7 Single stellar population PNLFs depending on the age and metallicity for the fiducial model without circumnebular extinction but varying the nebular metallicity correction. Rows correspond to the different ages, and columns to the metallicities. The final masses are computed from the IFMR of Miller Bertolami (2016) and are stated in the lower-right corner of each panel. In the cases where the initial mass values lie outside the modeled range, only a dash is shown (this is only the case for two panels in the top right and one in the bottom left). The PNLFs are obtained from 106 stars of the respective final mass with equally spaced post-AGB ages between 0 and 30 000 yr. For the black PNLFs no nebular metallicity correction factor was applied, whereas the colored lines include the corrections by Jacoby (1989), Dopita et al. (1992), and Schönberner et al. (2010). The bright-end cutoff is marked by the vertical gray dashed line at M(5007) = -4.5.

4.2 PNLFs with circumnebular extinction

The fact that the PNLFs discussed in the previous section are able to reach magnitudes much brighter than the observed cutoff value is very likely due to the PN model of Valenzuela et al. (2019) not incorporating circumnebular dust extinction. The circumnebular dust extinction has been observed to rise with increasing PN brightness, thus pushing the brightest PNe back down to the bright-end cutoff (e.g., Davis et al. 2018; Jacoby & Ciardullo 2025). In Fig. 8, we show the full grid of PNLFs with respect to their unextincted and dust-extincted [O III] magnitudes. The unextincted PNLFs (dotted) are the same as those shown in Fig. 7 for the non-metallicity dependent PNLFs (black) and those corrected for nebular metallicity (colored) through the relations of Dopita et al. (1992) and Schönberner et al. (2010). We do not show the PNLF by Jacoby (1989) to keep the figure from being too cluttered and since the simple square root relation is most likely the least accurate of the three. The solid PNLFs move to dimmer magnitudes by taking into account timedependent dust extinction according to the Combined OR linear relation for circumnebular extinction, which is stronger in the early post-AGB phase. As a result, the extinction is seen the strongest for the fast evolving PNe of younger ages (top rows), where for t_age = 0.25 Gyr none of the PNe reach the brightest two magnitudes of the universal PNLF anymore (M(5007) < −2.5). At older ages of t_age ≳ 6 Gyr, there is practically no more extinction that can be observed: the solid and dotted PNLFs are almost perfectly overlayed.

For lower-to-intermediate ages of 1-4 Gyr, the extinction affecting the brightest PNe is on the order of 0.25-0.75 mag and also varies slightly with metallicity. This is expected as the lifetime function and the IFMR both affect the resulting final masses and thus how fast they evolve. Figure 8 nicely shows that circumstellar extinction is able to lower the SSP cutoff magnitude to the actually observed cutoff value of M(5007) = -4.5 for the age range of 1-4 Gyr. At younger ages with generally higher final masses, the adopted prescription leads to much higher extinction for the brightest PNe and is on the order of 1.5-2.5 mag for all but the highest metallicity. Only at Z = 0.08 the lower final masses at young-to-intermediate ages lead to significantly lower circumnebular extinction of the bright end. The extent to which this is true naturally depends on the nebular metallicity correction applied, as for example the relation by Schönberner et al. (2010) leads to significantly brighter PNe at the metal-rich end of Z ≳ 0.02. Even the circumnebular extinction recipe does not bring down the brightest PNe to the bright end for intermediate-aged metal-rich populations (t_age ~ 0.56 Gyr, Z ≳ 0.02) for Schönberner et al. (2010). We conclude that the fiducial model involving the nebular model by Valenzuela et al. (2019) with the correction by Dopita et al. (1992) is the most realistic with respect to our expectations from observations as it properly accounts for metallicity and simultaneously does not over-enhance the brightness of PNe at high metallicities.

As a result, we find that correcting for nebular metallicity through the relation of Dopita et al. (1992) and additionally accounting for circumnebular extinction (Jacoby & Ciardullo 2025) leads to a grid of PNLFs that are in general consistent with observations with respect to the bright end. They offer the basis for interpreting the observed PNLFs of composite stellar populations. At the same time, the fact that the PNLFs vary significantly between different nebular metallicity corrections also indicates that improved modeling will be necessary to decrease the systematic uncertainties seen between methods. Developing a refined nebular model is beyond the scope of this work, however.

4.3 PNLFs with different IFMRs

In Sect. 2.3, we presented several IFMRs in addition to the one by Miller Bertolami (2016) used in the fiducial PICS model, namely those by Cummings et al. (2018), El-Badry et al. (2018), Cunningham et al. (2024), and Hollands et al. (2024). In Fig. 9 we show the consequence of adopting an IFMR that does not depend on metallicity by using the additional four IFMRs instead of that of the fiducial PICS model by Miller Bertolami (2016). The black lines are the same as the blackberry lines in Fig. 8 from the fiducial PICS model with the nebular metallicity correction by Dopita et al. (1992), while the colored lines correspond to the obtained PNLFs using the IFMRs of Cummings et al. (2018), El-Badry et al. (2018), Cunningham et al. (2024), and Hollands et al. (2024). The different IFMRs only show a similar PNLF for younger ages of around t_age ≲ 6 Gyr, and only the PICS fiducial IFMR shows a deviating behavior at higher metallicities in this age range with a greater number of PNe, but at the same time a fainter bright end of the PNLF. In the context of circumnebu-lar extinction, as the final masses now vary even within a single panel based on which IFMR was used, the effect of extinction is different for each PNLF because of the final mass influencing how fast the central star evolves. This can especially be seen for the young metal-rich panels, where the dotted and solid black lines lie almost on top of each other for the fiducial model using the IFMR by Miller Bertolami (2016), while the other IFMRs lead to much higher final masses, thus resulting in significantly stronger extinction at the bright end (the colored solid lines lie much further to the right than their dotted counterparts).

The values of age and metallicity where similar PNLFs are found for all IFMRs roughly correspond to the initial mass range of M_init ≳ 1.3 M_⊙ in which the alternative IFMRs do not lead to higher final masses than at even the lowest metallicities of the fiducial models (Fig. 2). In contrast, for lower initial masses and thus for older stellar populations, Fig. 2 shows that in particular the IFMR by Cummings et al. (2018) surpasses the final mass values of the fiducial model at all metallicities. At the lowest initial masses around M_init ≈ 0.9 M_⊙, both the models of Cummings et al. (2018) and Cunningham et al. (2024) show higher final masses than the fiducial models because their values result from linear fits in this mass range without a drop toward lower initial masses (in contrast to the fiducial model and that by Hollands et al. 2024). Since this leads to final masses of M_final ≳ 0.55 M_⊙, even such old low-mass stars are able to produce bright PNe, as seen by the bright PNe in turquoise and red in Fig. 9 for older systems of t_age ≳ 8 Gyr. Overall, the PNe produced with the IFMR by Cummings et al. (2018) can be the brightest at older ages. The PNLFs obtained with their IFMR are surprisingly similar in bright end and shape across most ages and metallicities with only minor shifts of the bright end. Among the other four IFMRs, differences in the bright-end cutoff of up to half a magnitude can be found between the different IFMRs for t_age = 6 Gyr, one magnitude at 8 Gyr, and more than 2-3 magnitudes in stellar populations older than 10 Gyr. Finally, for the IFMR by El-Badry et al. (2018) the especially low final masses obtained from initial masses of M_init ≲ 2 M_⊙ lead to PNe dropping out of the brightest magnitudes of the PNLF for t_age ≳ 4 Gyr for metallicities approximately up to the solar value (a lack of magenta PNLFs in the center to bottom left panels). At higher metallicities this is the case for t_age ≳ 6Gyr (Z = 0.04) and t_age ≳ 8 Gyr (Z = 0.08) as a result of the longer lifetimes of metal-rich stars leading to higher initial masses at a given age.

Despite the sharp drop in the IFMR of Hollands et al. (2024) toward lower masses in the range of M_init ≲ 1 M_⊙, which is even more extreme than the fiducial model shows, the light blue lines for Hollands et al. (2024) at high ages and metallicities show that bright PNe can still be produced (bottom right panels in Fig. 9). This is a consequence of the alternative IFMRs not depending on metallicity, while the metallicity still influences the lifetime function, where higher metallicities lead to more massive stars reaching the post-AGB phase for a given SSP age. Because of the drastic drop of the final mass with lower initial masses, old low-metallicity SSPs do not produce bright PNe when using the IFMR by Hollands et al. (2024), leading to a lack of the light blue PNLFs in the bottom left panels in Fig. 9.

For the younger ages at high metallicities (top right panels), it may seem curious that all alternative IFMRs lead to similar bright PNLFs before accounting for circumnebular extinction, whereas the fiducial IFMR by Miller Bertolami (2016) leads to a slightly dimmer bright end but a significantly higher number of PNe. The reason for this can be understood by considering Fig. 2. For young ages, the initial masses are above 2 M_⊙. In this range, the high metallicity fiducial IFMRs show much lower final masses than the alternative IFMRs. The final masses are still high enough to produce relatively bright PNe, but the PNe are not quite as bright as those obtained with higher final masses. However, the stars with lower final masses move more slowly across the Hertzsprung-Russell diagram before cooling down, thus allowing the nebula to survive for a longer period of time. The result is a greater number of PNe that can be observed.

The grid of PNLFs with different assumed IFMRs shows that the assumed IFMR strongly affects the PNLF that can be expected from a given stellar population, in particular for older systems. It is interesting that it is the regime of old stellar populations, where the abundance of PNe at the bright end has been a puzzle for a long time, that shows the highest sensitivity to the IFMR. This shows that understanding the IFMR dependence on metallicity (i.e., understanding the dependence of cold winds on metallicity) will play a significant role in the understanding of PNe and the PNLF in old stellar populations. We also note that here no potential intrinsic scatter of the IFMR is accounted for, which will also affect the resulting PN populations and the brightest PNe through the variations of the final mass (see also the discussion by Jacoby & Ciardullo 2025). Finally, even though the metallicity is not considered for the alternative IFMRs themselves, the metallicity dependence of the lifetime function and the post-AGB tracks leads to a shift of the bright end of the PNLFs, which can mostly be summarized as having a brighter PNLF cutoff at higher metallicities.

Fig. 8 Single stellar population PNLFs depending on the age and metallicity for the fiducial model with different nebular metallicity corrections applied and showing the effect on the PNLFs from circumnebular dust extinction. The structure of the grid and the amount of PNe sampled for each PNLF are the same as in Fig. 7. For the black PNLFs no nebular metallicity correction factor was applied, whereas the colored lines include the corrections by Dopita et al. (1992) and Schönberner et al. (2010). The dotted PNLFs are the same as those shown in Fig. 7 and the solid PNLFs additionally take into account the circumnebular extinction from Jacoby & Ciardullo (2025) as described in Sect. 2.6.

Fig. 9 Single stellar population PNLFs depending on the age and metallicity for the fiducial model using the nebular metallicity correction by Dopita et al. (1992) and with variations of the IFMR (colored). The IFMRs are from Miller Bertolami (2016), in black, and Cummings et al. (2018), El-Badry et al. (2018), Cunningham et al. (2024), and Hollands et al. (2024), in turquoise, magenta, red, and light blue, respectively. The structure of the grid and the amount of PNe sampled for each PNLF are the same as in Fig. 7. The dotted PNLFs are constructed from the intrinsic [O III] magnitudes and the solid PNLFs additionally take into account the circumnebular extinction from Jacoby & Ciardullo (2025) as described in Sect. 2.6.

4.4 Normalized PNLFs

While the PNLF grids shown until this point were each constructed for a fixed number of 10⁶ PNe and are useful for investigating the parameter space of the PICS model, their normalization still lacks the direct comparability to observations with respect to the underlying stellar populations. Since the PN central stars have different initial and final masses depending on the age and metallicity (the final masses are noted in the bottom left corner of each panel in Fig. 7), varying numbers of PNe are expected to be found if the total SSP mass or luminosity is kept constant, based on the IMF. In Fig. 10, we show the Dopita et al. (1992) nebular metallicity-corrected and circumnebular-extincted PNLFs for a fixed number of 10⁶ PNe as shown in the other PNLF grids (black), for a fixed total stellar mass of M_* = 10¹² M_⊙ assuming a Chabrier (2003) IMF (red), and for a fixed total bolometric luminosity of L_bol = 10¹² L_⊙ based on the GALAXEV CB07 SSP models by Bruzual & Charlot (2003) using a Chabrier (2003) IMF (orange).

For the luminosity-normalized PNLFs (orange), there are more PNe than for the fixed-number PNLFs (black) at older ages relative to the number of PNe at younger ages. The reason for this arises from the fact that younger stellar populations are intrinsically brighter than older ones, thus resulting in less total mass when keeping the total bolometric luminosity constant. In turn, this leads to a smaller number of PNe for young ages. The normalization to bolometric luminosity also gives us the possibility to determine the luminosity-specific PN numbers, typically denoted by α. As our model produces complete samples of PNe for a given SSP, it is possible to compute the actual value of α without having to assume an underlying PNLF shape based on which the number is estimated by extrapolation. Instead, the α value is obtained directly from counting the number of PNe in a given SSP brighter than a given magnitude. As the model can produce extremely dim magnitudes for objects that are essentially not PNe anymore, we limited the α determination to PNe with M(5007) ≤ +3.5, which corresponds to the 8 brightestmag-nitudes with an assumed bright-end cutoff of M(5007) = -4.5. We chose this magnitude range to align with the estimation of Buzzoni et al. (2006) on PNe spanning roughly 8 orders of magnitude. The values of α for the same grid of SSPs found in Fig. 10 are shown in Table 2, using the extincted PICS fiducial model and the same bolometric luminosity prescription by Bruzual & Charlot (2003) as described above. We note that this approach differs from that of Buzzoni et al. (2006) through our full modeling of PNe: we can simply count the PNe and do not have to make assumptions on the specific evolutionary flux, but rather can infer it from the model results themselves (see below).

These values for the most part lie between 5 × 10^-7 L_⊙^-1 for solar metallicities of Z = 0.01-0.02 and ages of t_age ≳ 2 Gyr, which drop toward lower metallicities and ages, and rise toward higher metallicities and ages. These results are generally on the same order of magnitude as what Buzzoni et al. (2006) found, except that they find α to correlate much more strongly with age and much less with metallicity. The trend with metallicity that we find is largely due to the metallicity-dependent stellar evolution models employed from Miller Bertolami (2016) with respect to the stellar lifetimes and the IFMR, which was not taken into account by Buzzoni et al. (2006). Our smaller trend with age, however, could arise from the assumption of the nebular model from Valenzuela et al. (2019) that PNe have a maximum lifetime of 30 000 yr. While this assumption is likely sufficient for most observed PNe, the dimmest PNe may very well also occur after this time frame. This could increase the total number of PNe for older stellar populations since their central stars have lower masses and thus evolve more slowly along the stellar tracks, leading to a slower loss of temperature and luminosity.

As Buzzoni et al. (2006) discuss, the α parameter can be linked to the PN visibility lifetime τ_PN. For the PICS model, it is possible to directly extract this lifetime for each age and metal-licity as the time that a PN has a magnitude of M(5007) ≤ +3.5. From the 10⁶ modeled PNe uniformly sampled across 30 000 yr for each SSP, we took the fraction of PNe in that magnitude range as the factor to be multiplied with 30 000 yr to obtain the respective PN visibility lifetime, which we show in Table 3 for the same SSP grid. Around solar metallicity of Z = 0.010.02, the visibility lifetimes roughly range between 25 000 and 28 000 yr for SSP ages t_age ≥ 1 Gyr. Lower metallicities and ages lead to shorter PN visibility lifetimes, similar to the trends seen with α as expected. We note that the values do not reach the full 30 000 yr because of the stochastic nature of the model by Valenzuela et al. (2019) through the absorption parameter, which accounts for possible optical transparency of the nebula. This can drop the [O III] flux below the 8 magnitudes for the faintest PNe. Further, these values of the PN visibility lifetime are affected by the model’s choice of 30 000 yr, a matter that will be addressed through improved nebular models in the future. From α and τ_PN the specific evolutionary flux can be determined as their fraction B = α/τ_PN (Buzzoni et al. 2006). For instance, at Z = 0.01 and t_age = 6.0 Gyr, the specific evolutionary flux is B = 1.95 × 10^-11 L_⊙^-1 yr^-1, which exactly lies in the upper range shown in fig. 1 from Buzzoni et al. (2006).

In the literature the use of α_2.5 has commonly been used due to the fact that observationally one is generally restricted to only the brightest magnitudes of the PNLF (e.g., Jacoby et al. 1990; Longobardi et al. 2013; Hartke et al. 2017). Hereby only the 2.5 brightest magnitudes are taken into account for the luminosityspecific PN number. We show the predicted values of α_2.5 and their corresponding visibility lifetimes τ_PN,2.5 in Appendix E.

For some nearby galaxies it has been possible to measure a statistically complete PNLF down to six magnitudes below the bright cutoff (e.g., Reid & Parker 2010; Bhattacharya et al. 2019), which has revealed significant variations of the PNLF shape at its faint end. We show the same grid of normalized PNLFs as in Fig. 10 up to a magnitude of M(5007) = +2 in Appendix F, which is around 6.5 magnitudes below the bright end of the PNLF. These will be relevant in upcoming studies to compare the models with such nearby deep observations and provide a means for interpreting the variations of the PNLF shape. We highlight two insights from the extended magnitude PNLF grid: first, we find the dip or “camel” shape previously found in some galaxies with recent star formation like the Magellanic Clouds or NGC 300 (Jacoby & De Marco 2002; Reid & Parker 2010; Peña et al. 2012; Soemitro et al. 2023) and simulations of mock stellar populations with more massive central stars (Méndez et al. 2008; Valenzuela et al. 2019). Here, we find that this dip already arises from the PNLF of a sole SSP and is not necessarily the result of multiple stellar populations overlapping with each other, as has for example been proposed by Rodríguez-González et al. (2015). As proposed by Ciardullo (2010), the physical reason for the dip is that more massive post-AGB stars can host brighter PNe on their heating track than lower-mass stars, but their extremely fast evolution leads to them cooling down very quickly, resulting in a rapid decline of the [O III] emission before reaching a moderate brightness, where their evolution slows down and PNe accumulate, which is where the PNLF rises again. Second, the differences in shape and slope found at the faint end of the PNLF in observations indicate that these individual SSP PNLFs will help in better understanding and interpreting the different behaviors found at the faint end of the PNLF, for example for PN populations in M 31 and its outskirts (Bhattacharya et al. 2021).

As a final comment, since any given stellar population is made up of many SSPs, any observed or simulated PNLF is in fact a superposition of these SSP PNLFs. In an upcoming study, we will use this fact to present a test using resolved star formation histories of observed galaxies to compare the PICS-predicted PNLF based on the observationally inferred stellar populations to the actually observed PNLF. With those tests we will demonstrate how successfully PICS is capable of modeling accurate PN populations self-consistently.

Fig. 10 Single stellar population PNLFs depending on the age and metallicity for the fiducial model with the nebular metallicity correction by Dopita et al. (1992) and circumnebular extinction by Jacoby & Ciardullo (2025), normalized to different quantities. The structure of the grid and the amount of PNe sampled for each PNLF are the same as in Fig. 7. The black lines are the PNLFs obtained from 106 stars of the respective final mass with equally spaced post-AGB ages between 0 and 30 000 yr. The red PNLFs are normalized to a total mass of 1012 M⊙ assuming a Chabrier IMF (Chabrier 2003). The orange PNLFs are normalized to a total bolometric luminosity of 1012 L⊙ based on the GALAXEV CB07 SSP models by Bruzual & Charlot (2003) using a Chabrier IMF.

5 Summary and conclusion

This first paper of a series has introduced PICS, a modular framework that models the PN population contained within an SSP. The SSP is described by its total mass, age, and IMF and optionally by its metallicity and/or further abundances. For our fiducial model, we considered the SSP metallicity mass ratio and used the metallicity-dependent lifetime function, IFMR, and post-AGB tracks from Miller Bertolami (2016) to model the central stars of the PNe. The nebula and its [O III] emission were modeled using the method from Valenzuela et al. (2019), and we accounted for the nebular metallicity according to the results by Dopita et al. (1992). Finally, time-dependent circumnebular dust extinction was applied to the [O III] emission derived from the observed relation by Jacoby & Ciardullo (2025). PICS thus produces the main observable quantity of extragalactic PNe and can be employed to compare theoretical models and hydrodynamical cosmological simulations with observations, which can in turn be used to interpret the properties of the PNe and underlying stellar populations.

We have shown how the initial metallicity of the star affects the stellar evolution relevant to PNe. Namely, higher metallic-ities lead to more massive stars reaching the post-AGB phase later than their lower-metallicity counterparts, while the final masses after the AGB phase are lower at a constant initial mass. Higher metallicities also lead to lower luminosities and temperatures of the central stars for a fixed final mass. Considering both the lifetime function and the IFMR, the effect of the metallic-ity is strongest at ages below 3-4 Gyr, where lower metallicities lead to more massive PN central stars, and at ages above 6 Gyr, where higher metallicities lead to more massive central stars (Fig. 5). The effect on the final mass becomes strongest at very young and at very old ages. The final mass is important for the maximal brightness of M(5007), which can be reached throughout the post-AGB phase, and reaching higher final masses will increase the probability of finding such a PN at the bright end of the PNLF. Furthermore, these relations provide the first indications that the metallicity is relevant for older stellar populations, as these populations are oftentimes more metal rich in massive elliptical galaxies than in disk galaxies such as the Milky Way.

One important finding of this work is that for high-metallicity stars, the metallicity is not the only property affecting the lifetime of stars. While higher metallicities generally do increase the opacity of the star, leading to longer lifetimes, the accompanying increase of the helium abundance counteracts this trend due to higher MS luminosities and less available fuel. As the relation of helium abundance with metallicity is still poorly constrained from observations, particularly at high metallicities, we explored the parameter space ranging from a linearly increasing helium abundance with metallicity to a constant helium abundance of Y = 0.285 above a metallicity of Z = 0.02. Such a constant helium abundance above a certain metallicity is in line with recent observations of the metallicity of the stellar population of the globular cluster Omega Centauri (Clontz et al. 2025).

Assuming the linear relation, the lifetime function does not lead to longer timescales at higher metallicities above Z = 0.02 anymore, whereas assuming a saturation of the helium abundance results in longer lifetimes at Z > 0.02. The helium abundance relation is therefore especially relevant for metalrich stellar populations, as it governs how the lifetime function behaves for the individual stars. A higher helium abundance causes shorter lifetimes, leading to less massive stars becoming PNe after a given time at a fixed metallicity. This in turn means that the brightest possible magnitudes M(5007) are dimmer. If the helium abundance is capped, however, the stellar lifetimes increase, leading to more massive central stars at older ages for super-solar metallicity stellar populations. Based on the metallicities originally used by Miller Bertolami (2016) of Z = 0.0001-0.02, we derived an analytic expression for the lifetime function, whose extrapolation lies between the simulated values for linear and constant helium abundance relations but is closer to the latter. Thus, for the fiducial model, we picked the constant helium abundance relation, as it lies closer to the analytic function. Since longer lifetimes of metal-rich stars increase the occurrence of bright-end PNe in old stellar populations, our findings indicate that a saturation of the helium abundance at high metallicities could be part of the solution for the universal bright-end cutoff of the PNLF. The implications would be that high-metallicity stars have longer lifetimes than expected from models assuming non-saturated helium abundances. In future studies, we will extend our stellar grids to include the effect of the helium abundance consistently with the prediction from cosmological models.

To obtain a complete overview of how the SSP metallic-ity and age affect the resulting PN population, we produced a grid of PNLFs, each corresponding to a fixed pair of metallic-ity and age. The PNLF shapes are generally consistent with each other and have a varying bright-end cutoff. In particular, older SSPs generally contain less luminous PNe. However, at older ages above 6 Gyr, higher metallicities again lead to brighter PNe forming. In addition, at almost all ages, the lowest metallicities of Z = 0.0001 and Z = 0.001 lead to a significant dimming of the PNe through any of the three nebular metallicity corrections we implemented in this work (Jacoby 1989; Dopita et al. 1992; Schönberner et al. 2010), especially compared to not accounting for the nebular metallicity. Notably, the brighter old and metalrich PNe could in part explain how old elliptical galaxies still retain a similar bright end of the PNLF compared to disk galaxies, as they tend to be especially metal rich. The grids of SSP

PNLFs presented in this work set the theoretical foundation for interpreting the relations of the observed PNLF bright end with properties of the underlying stellar populations with respect to their age and metallicity.

We have also presented a grid of extended PNLFs down to six magnitudes below the bright-end cutoff, which are needed for comparison and interpretation of PN studies of Local Universe galaxies, where it has been possible to obtain statistically complete samples even of very faint PNe, for example in the LMC or M31 (Reid & Parker 2010; Bhattacharya et al. 2021). The extended PNLFs revealed that the “dip” or “camel shape” of the PNLF with a minimum at around M(5007) ≈ -2-0, which has previously been seen in observations (e.g., Jacoby & De Marco 2002; Reid & Parker 2010; Peña et al. 2012; Soemitro et al. 2023) and models (e.g., Méndez et al. 2008; Valenzuela et al. 2019), is in fact an intrinsic property of PNe with more massive central stars and arises even in SSPs. It is thereby not merely a result of two or more modes of PNLFs from multiple stellar populations, as has been suggested by Rodríguez-González et al. (2015). Of course, any real PNLF is a superposition of SSP PNLFs, which can lead to an even larger variety of observed faint PNLF shapes. This carries the potential of being another independent diagnostic tool for investigating the underlying star formation history and metallicity properties of galaxies and their formation histories, which we will also be presenting with a self-consistent test on observations in an upcoming study.

Through the PNLF grids, we investigated the effects of applying different nebular metallicity corrections, which are critical to avoid overly luminous, extremely metal-poor PNe. The square root law by Jacoby (1989) is the relation first found for solar-like metallicities. It enhances the brightness of metal-rich PNe very strongly when extrapolated and leads to the smallest dimming of metal-poor PNe. In contrast, the relations by Dopita et al. (1992) and Schönberner et al. (2010) enhance it less significantly at the metal-rich end, and the relation by Dopita et al. 1992 offers the least enhancement. Applying circumnebular extinction further decreases the brightness of the brightest PNe of very young stellar populations such that the bright end of the PNLF can only be found at solar-like and higher metallicities for stellar ages of ≳1 Gyr.

Finally, we studied how different IFMRs affect the resulting PN populations by considering metallicity-independent IFMRs determined from the solar neighborhood in the Milky Way in addition to the relation used for the fiducial model from Miller Bertolami (2016). Most of the various IFMRs (except that by El-Badry et al. 2018) do not change the SSP PNLFs much for younger stellar populations of ≲6 Gyr. However, at higher ages the metallicity-independent IFMRs mostly lead to brighter PNe being produced, and some even reach the bright-end cutoff for the IFMR by Cummings et al. (2018). The reason for this is that some of these IFMRs do not experience a final stellar mass drop as steep toward low-initial masses as the relations by Miller Bertolami (2016). Only the IFMR by El-Badry et al. (2018) predicts much lower final masses in the low-mass regime, that is for old stellar ages. As a result, no PNe are produced in the four brightest magnitudes for ages of ≳6 Gyr for this IFMR. Furthermore, the lack of metallicity dependence means that at high metallicities, the PN central stars are more massive than when employing the metallicity-dependent IFMR, in turn leading to an on average brighter PNe. The implications for the existence of bright PNe in old stellar populations are very important. The IFMR, in particular at initial masses of M_init ≲ 1.2M_⊙, plays a significant role in how massive the central stars are and therefore how bright the PNe can get. Obtaining better constraints and developing a more robust theoretical understanding of the mass loss in the AGB phase is therefore essential for explaining the observed bright end of the PNLF in old stellar populations.

With PICS, we have laid the groundwork for populating cosmological simulations with PNe, thereby producing PN populations based on self-consistent star formation histories within a cosmological context. On the theoretical side, this allows for predictions of the PNLF behavior to be made based on statistically significant samples of galaxies with diverse morphologies and star formation histories, which will aid in further understanding the universal bright end of the PNLF. Furthermore, the modular structure of PICS makes parameter studies where controlled changes to one specific block in the pipeline can be used to understand its effect on the resulting PN population possible, just as we did for the IFMR in this work. For the observational side, it will be possible to compare the modeled PN populations to observations by identifying appropriate simulated galaxy analogs and to offer an interpretation for the observed PNe and their properties.

To conclude, we have shown in this study that it is essential to account for both the metallicity and helium abundances in PNe models as well as to improve the observations and models of the IFMR. These factors all affect the brightest PNe in particular, bringing the community one step further in our understanding of the universal bright end of the PNLF, which we will explore in more detail in future studies.

Acknowledgements

We thank Souradeep Bhattacharya, Tadziu Hoffmann, George Jacoby, Rolf-Peter Kudritzki, and Joachim Puls for helpful comments and discussions, and Tim Cunningham for providing us with his data tables of the IFMR. We thank the referee, Robin Ciardullo, for his constructive feedback that helped us significantly improve the paper. LMV acknowledges support by the German Academic Scholarship Foundation (Studienstiftung des deutschen Volkes) and the Marianne-Plehn-Program of the Elite Network of Bavaria. M3B is partially funded by CONICET and Agencia I+D+i through grants PIP-2971 and PICT 2020-03316, and by CONICET-DAAD 2022 bilateral cooperation grant number 80726. This research was supported by the Excellence Cluster ORIGINS, funded by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy - EXC-2094-390783311. The following software was used for this work: Julia (Bezanson et al. 2017), DataInterpolations.jl (Bhagavan et al. 2024), Interpolations.jl (Holy et al. 2022), matplotlib (Hunter 2007), and numpy (Harris et al. 2020).

Appendix A Post-AGB track interpolation

Here we present how the post-AGB tracks provided by Miller Bertolami (2016) are interpolated (see Sect. 2.4 for the corresponding section in the main text). For the post-AGB tracks, we take the output tables from Miller Bertolami (2016), which each consist of 101 points of equivalent position along the tracks in the Hertzsprung-Russel diagram, with precomputed luminosity, effective temperature, and post-AGB time since a temperature of log T_eff/K = 3.85 was reached. The tables are available for between 12 and 20 final masses for each of the four metallicities Z = 0.0001 (12 between M_final = 0.525 and 0.800 M_⊙), Z = 0.001, Z = 0.01, and Z = 0.02 (each of them have 20 tables for final masses between M_final = 0.525 and 1.000 M_⊙). As we follow the approach of Valenzuela et al. (2019) in this work to define the post-AGB age as the time since a temperature of T_eff = 25 000 K was reached, the precomputed post-AGB ages are shifted accordingly in the individual tables. To interpolate the tracks for arbitrary metallicities, final masses, and post-AGB ages, we perform linear interpolations on the logarithmic metallicity and the final mass to find the eight relevant equivalent points on the tracks corresponding to the two adjacent metallicities, final masses, and post-AGB ages. The luminosity and effective temperatures are then determined from these eight points by trilinear interpolation, the result of which can be seen in Fig. 3.

Appendix B Nebular model recipes

The nebular model by Valenzuela et al. (2019) as described in Sect. 2.5 employs empirical relations for the intensity ratio between I(5007) and I(Hβ) and for the absorption factor that accounts for the possibility of optically thin nebulae and their dissipation over time. In this appendix we describe these recipes in more detail to present the final model from Valenzuela et al. (2019), which was built on previous versions by Méndez et al. (1993) and Méndez & Soffner (1997).

B.1 Intensity ratio recipe

After having obtained I(Hβ) assuming a blackbody spectrum based on the luminosity and effective temperature of the central star, the relative intensity of I(5007) is determined empirically. For this we consider three different cases based on T_eff:

• T_eff ≤ 30152 K: I(5007) is set to zero.

• 30152K < T_eff ≤ 60000K: The intensity is set to 100I(5007)/I(Hβ) = 0.0396T_eff/K - 1194 based on the linear relation derived by Méndez et al. (1993) in their eq. 5.

• T_eff > 60 000 K: First, a random value I_rel is drawn from a normal distribution with a mean of 9.50 and a standard deviation of 3.75, truncated to 0 if a negative value is drawn. As Valenzuela et al. (2019) describe, the steep drop-off at the high end of the distribution has to be taken into account to avoid overly high values of I(5007): First, for I_rel > 12 the value is decreased by ((I_rel - 12)/8)^1.5. If it is still larger than 15, it is then further decreased by (I_rel - 15)/5. Finally, if the value is still larger than 18, I_rel is simply set to a uniformly random value between 18 and 19 and then truncated to a maximum value of 18.5. In a last step, for low central star masses of M_final < 0.55 M_⊙ that are hot and luminous with log(L/L) > 3 and T_eff > 75 000 K, the value of I_rel is multiplied with 0.4 as described by Valenzuela et al. (2019). This value of I_rel is finally used to obtain I(5007) by I(5007) = I_relI(Hβ).

B.2 Absorption factor recipe

The absorption factor μ used by Valenzuela et al. (2019) is determined through a linear combination based on central star mass, effective temperature, and post-AGB age relative to reaching the highest combination of effective temperature and luminosity at the knee of the post-AGB tracks. The post-AGB time to the knee is approximated through the analytical relation $t_{\mathrm{k}\mathrm{n}\mathrm{e}\mathrm{e}}=0.0239M_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}^{-21.82}$. The absorption factor is a value between 0 and 1 that is multiplied with the originally obtained value for I(5007) to account for a lack of optical thickness and nebular dissipation. First, a generator value μ_gen is randomly drawn between 0 and 1 from a uniform distribution. For the cooling tracks with post-AGB ages larger than 1.1 times the time to reach the knee, the absorption factor is determined as $\mu =(0.1+0.9{\mu }_{\mathrm{g}\mathrm{e}\mathrm{n}}^{0.9})\times (1-t_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}-\mathrm{A}\mathrm{G}\mathrm{B}}/30000\mathrm{y}\mathrm{r})$ to account for the nebulae becoming transparent over time. For the heating tracks until the knee, three temperature ranges are considered, where linear interpolation is applied for temperatures between the respective limits:

• log T_eff/K < 4.51: For the beginning of the heating track, μ is set to min(0.4 + μ_gen, 1) to account for generally more opaque PNe.

• 4.68 < log T_eff/K < 4.85: For low-mass central stars below 0.535 M_⊙, μ = 0.03 + 0.17μ_gen is used to obtain overall low values of μ because of their slow post-AGB evolution, such that the nebulae already dissipate during this time. In contrast, for higher-mass central stars above 0.557 M_⊙, $\mu =0.2+0.8(1-{\mu }_{\mathrm{g}\mathrm{e}\mathrm{n}}^2)$ is applied, which allows for high values of μ. Linear interpolation is applied for the mass range in between these mass limits.

• log T_eff/K > 4.95: For low-mass central stars below 0.535 M_⊙, μ = 0.03 + 0.17μ_gen is again used to account for their slow evolution along the post-AGB tracks. For highermass central stars above 0.557 M_⊙, μ = 0.1 + 0.9μ_gen is applied, which now leads to lower values of μ with a higher probability compared to at lower effective temperatures due to the time that has passed. Linear interpolation is again applied for the mass range in between these mass limits.

For the time between 1 and 1.1 times the time to reach the knee of the post-AGB track, and for the effective temperatures between the ranges listed above, linear interpolation is applied between the respective obtained values for μ. As a final step, μ is truncated to the range between 10^-8 and 1.

Appendix C Nebular metallicity

In Fig. C.1 we show the nebular metallicity correction factor p₅₀₀₇, which is multiplied by the obtained [O III] fluxes from the model by Valenzuela et al. (2019), as a function of metallicity. The relations are taken from Jacoby (1989), Dopita et al. (1992), and Schönberner et al. (2010) and are described in detail in Sect. 2.5. If no metallicity dependence is taken into account, the flux remains the same with a factor of 1 (black line). As we selected the base metallicity Z_base = 0.007 based on which the corrections are computed, all lines run through p₅₀₀₇ at that metallicity. Because the values of p₅₀₀₇ do not lead to as heavily boosted or dimmed [O III] fluxes, we selected the relation by Dopita et al. (1992) for the fiducial PICS model.

Appendix D Circumnebular extinction

The relations for circumnebular extinction by Jacoby & Ciardullo (2025) are only derived from the brightest PNe of M 31 and the LMC, and they are only dependent on the final stellar mass. In Sect. 2.6, we describe our choice of using PNe of metallicity Z = 0.01 to obtain the post-AGB ages at which the maximal brightness can be reached. Here we present the detailed derivation of this relation. We first modeled the intrinsic [O III] flux evolution throughout the PN lifetime of 30 000 yr for final masses M_final = 0.53-1.00 M_⊙ at Z = 0.01, taking at each timestep the maximal brightness based on 10 000 realizations of the nebular model. For each final mass, we then determined the post-AGB age t_post-AGB at which the highest flux is reached. We selected 50 equally spaced values in the final mass range, for which the brightest post-AGB ages are shown in Fig. D.1, using the stellar evolution relations by Miller Bertolami (2016) and the nebular model by Valenzuela et al. (2019).

Fig. C.1 Nebular metallicity correction factor p5007 as a function of metallicity for the relations given by Jacoby (1989), Dopita et al. (1992), and Schönberner et al. (2010), as well as for having no metallicity dependence. The relations are set up to reach a value of p5007 = 1 at the base metallicity Zbase = 0.007.

Fig. D.1 Post-AGB age tpost-AGB at which PNe at metallicity Z = 0.01 reach the brightest possible [O III] flux as a function of final mass Mfinal. The black data points mark the determined values from the model for 50 equally spaced final masses between 0.53 and 1.00M, and the red and orange lines are the linear fits to the values below and above Mfinal = 0.60 M⊙, respectively.

As the relation for the logarithmic brightest post-AGB ages appears to be well-described by a broken-line fit, we fit two linear relations to the final masses smaller and larger than 0.60 M_⊙ (red and orange lines, respectively). Since in our recipe we need to invert this relation to obtain the corresponding final mass for a given post-AGB age, we perform the fits as a function of postAGB age. The analytical form is given in Eq. 11. This approach ensures that a PN at its maximum brightness experiences the extinction from the derived relation from Jacoby & Ciardullo (2025), while simultaneously accounting for dust dissipating over time.

Appendix E Luminosity-specific bright PN numbers

In Sect. 4.4 we presented the luminosity-specific PN numbers (α) and the visibility lifetimes (τ_PN) for all PNe with magnitudes brighter than M(5007) = +3.5, as obtained for the individual SSPs. This limiting magnitude is based on the assumption that PNe span roughly 8 magnitudes of brightness (e.g., Buzzoni et al. 2006). However, in practice it is impossible to measure the entire population of PNe in a galaxy, where for most galaxies outside the Local Group one is limited to the brightest couple of magnitudes. For this reason, α_2.5 has become a relevant measure for the luminosity-specific PN number for those PNe with magnitudes within the brightest two and a half magnitudes of the PNLF (e.g., Jacoby et al. 1990; Buzzoni et al. 2006; Longobardi et al. 2013; Hartke et al. 2017). For the values of each of the SSP PNLFs, we assume a fixed PNLF bright-end cutoff of M(5007) = -4.5 and count all PNe brighter with M(5007) ≤ -2.0. In Table E.1 we show the values of α_2.5 for the SSP PNLFs, analogously to Table 2 with the bolometric luminosities obtained through the GALAXEV CB07 SSP models of Bruzual & Charlot (2003) using a Chabrier (2003) IMF. Compared to α, these values are necessarily lower as the magnitude range is smaller. The values are generally compatible with results from observational work, where typically observed values are on the order of α_2.5 ~ 10^-8 L_⊙^-1 (e.g., Longobardi et al. 2013), which corresponds to values for ages of t_age ~ 10-12 Gyr, depending on metallicity (see Table E.1). The corresponding PN visibility lifetimes τ_PN,2.5 are shown in Table E.2, which indicate the time a PN spends being brighter than M(5007) = -2.0. These are also smaller than the values of τ_PN (Table 3) by necessity.

Appendix F Extended PNLFs

In Fig. F.1 we show the same set of SSP PNLFs as in Fig. 10 for a grid of different ages and metallicities, just increasing the plotted magnitude range to fainter PNe at M(5007) = +2. Here it is much more obvious that most of the individual PNLFs drop toward fainter PNe and many in the young to intermediate age range (t_age ≲ 6 Gyr) then again have a rising number of PNe toward the faintest magnitudes. This “camel” shape has previously been found in both observations and simulations of PNe, in particular in younger systems (e.g., Jacoby & De Marco 2002; Reid & Parker 2010; Valenzuela et al. 2019), where the dip in the SMC occurs at a slightly dimmer magnitude of M(5007) ≈ -0.5 (Jacoby & De Marco 2002). This is consistent with the grid of PNLFs shown here, where the SSP models in Fig. F.1 show that this shape can also arise within a pure SSP and it is not dependent on the correct superposition of different SSP PNLFs. In contrast, it is possible to wash out the camel shape by combining multiple SSP PNLFs as their peaks are located at different magnitudes and their shapes also differ in declining and increasing slope. The extreme rise in PN number at dim magnitudes for the youngest stellar populations of t_age ≲ 0.5 Gyr arises from the circumnebular extinction decreasing with post-AGB age. This means that the extinction affects the intrinsically brightest PNe in the early post-AGB phase the strongest. In contrast, the older, but much dimmer PNe are hardly affected by the extinction. The PNLF is thereby effectively squashed toward dim magnitudes.

Fig. F.1 Single stellar population PNLFs depending on the age and metallicity for the fiducial model with the nebular metallicity correction by Dopita et al. (1992) and circumnebular extinction by Jacoby & Ciardullo (2025), normalized to different quantities. This grid of PNLFs shows the same plots as seen in Fig. 10, except that the λ5007 magnitude range is extended to dimmer PNe by four magnitudes. The black lines are the PNLFs obtained from 106 stars of the respective final mass with equally spaced post-AGB ages between 0 and 30 000 yr. The red PNLFs are normalized to a total mass of 1012 M⊙ assuming a Chabrier IMF (Chabrier 2003). The orange PNLFs are normalized to a total bolometric luminosity of 1012 L⊙ based on the GALAXEV CB07 SSP models by Bruzual & Charlot (2003) using a Chabrier IMF.

References

All Tables

All Figures

Fig. 1 Lifetime function of stars until they reach the post-AGB phase as a function of their initial mass. Top: linear axes. Bottom: logarithmic axes. The logarithmic scale makes it easier to distinguish the individual lines and data points. The points are from the models by Miller Bertolami (2016), and the lines are the fit relations to those data points. The thick lines are at the same metallicities as the ones from the original models (Z = 0.0001, 0.001, 0.01, and 0.02), and the thin lines are at Z = 0.00001, 0.0003, 0.003, and 0.014, where the last three values are the logarithmic midpoints of the original metallicities.

In the text

Fig. 2 Initial-to-final mass relation of stars. The points are from the models by Miller Bertolami (2016), and the lines in the corresponding colors are the fit relations to those data points. The thick lines are at the same metallicities as the ones from the original models (Z = 0.0001, 0.001, 0.01, and 0.02) and the thin lines are at Z = 0.00001, 0.0003, 0.003, 0.014, 0.05, where the middle three values are the logarithmic midpoints of the original metallicities. Additionally, four observationally based IFMRs from nearby WD measurements at solar metallicity are shown from Cummings et al. (2018), El-Badry et al. (2018), Cunningham et al. (2024), and Hollands et al. (2024), as well as the 1σ uncertainty band of Cunningham et al. (2024).

In the text

Fig. 3 Post-AGB tracks in the Hertzsprung-Russell diagram of stars at two final masses Mfinal = 0.55 M⊙ and 0.58 M⊙ for a range of metallicities. The tracks are interpolated between the data tables from Miller Bertolami (2016). The metallicity colorbar corresponds to the same coloring as found in Figs. 1 and 2. The thick lines are at the same metallicities as the ones from the original models (Z = 0.0001, 0.001, 0.01, and 0.02) and the thin lines are at Z = 0.00001, 0.0003, 0.003, 0.014, 0.05, where the middle three values are the logarithmic midpoints of the original metallicities.

In the text

Fig. 4 Post-AGB tracks in the Hertzsprung-Russell diagram of stars at metallicities Z = 0.001, 0.01, and 0.02 (left, middle, and right, respectively) colored by the λ5007 magnitude. The tracks are shown for stars at different final masses Mfinal = 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.60, 0.65, 0.70, 0.80, 0.90, and 1.00M. For each final mass, stars at 600 evenly spaced post-AGB ages between 0 and 30 000 yr are plotted. The white lines indicate the tracks for the selected central star masses.

In the text

Fig. 5 Final masses as a function of stellar lifetime (left) and of initial mass (right) for selected metallicities (Z = 0.0001, 0.001, 0.01, 0.02, and 0.05) based on the models of Miller Bertolami (2016). The final mass y-axis range is reduced to the primarily relevant values in the regime of older stellar populations of t > 2 Gyr. For the initial-to-final mass relation, the masses reaching the post-AGB phase after t = 2, 4, 6, 8, 10, and 12 Gyr are marked with circles.

In the text

Fig. 6 Lifetime function of stars until they reach the post-AGB phase as a function of their initial mass. The lines are the fit relations to the models of Miller Bertolami (2016) with original metallicities of Z = 0.0001, 0.001, 0.01, and 0.02, as in Fig. 1. Here, the extrapolation of these fits to Z = 0.05 is also shown. The points are the results of newly run models following the procedure of Miller Bertolami (2016) for Z = 0.05, and the lines are the fit relations to those data points. The linear helium abundance is at Y = 0.345 and the cutoff helium abundance at Y = 0.285. The thick lines are at the same metallicities as the ones from the original models (Z = 0.0001, 0.001, 0.01, and 0.02) and the thin lines are at Z = 0.00001, 0.0003, 0.003, and 0.014, where the last three values are the logarithmic midpoints of the original metallicities.

In the text

Fig. 7 Single stellar population PNLFs depending on the age and metallicity for the fiducial model without circumnebular extinction but varying the nebular metallicity correction. Rows correspond to the different ages, and columns to the metallicities. The final masses are computed from the IFMR of Miller Bertolami (2016) and are stated in the lower-right corner of each panel. In the cases where the initial mass values lie outside the modeled range, only a dash is shown (this is only the case for two panels in the top right and one in the bottom left). The PNLFs are obtained from 106 stars of the respective final mass with equally spaced post-AGB ages between 0 and 30 000 yr. For the black PNLFs no nebular metallicity correction factor was applied, whereas the colored lines include the corrections by Jacoby (1989), Dopita et al. (1992), and Schönberner et al. (2010). The bright-end cutoff is marked by the vertical gray dashed line at M(5007) = -4.5.

In the text

Fig. 8 Single stellar population PNLFs depending on the age and metallicity for the fiducial model with different nebular metallicity corrections applied and showing the effect on the PNLFs from circumnebular dust extinction. The structure of the grid and the amount of PNe sampled for each PNLF are the same as in Fig. 7. For the black PNLFs no nebular metallicity correction factor was applied, whereas the colored lines include the corrections by Dopita et al. (1992) and Schönberner et al. (2010). The dotted PNLFs are the same as those shown in Fig. 7 and the solid PNLFs additionally take into account the circumnebular extinction from Jacoby & Ciardullo (2025) as described in Sect. 2.6.

In the text

Fig. 9 Single stellar population PNLFs depending on the age and metallicity for the fiducial model using the nebular metallicity correction by Dopita et al. (1992) and with variations of the IFMR (colored). The IFMRs are from Miller Bertolami (2016), in black, and Cummings et al. (2018), El-Badry et al. (2018), Cunningham et al. (2024), and Hollands et al. (2024), in turquoise, magenta, red, and light blue, respectively. The structure of the grid and the amount of PNe sampled for each PNLF are the same as in Fig. 7. The dotted PNLFs are constructed from the intrinsic [O III] magnitudes and the solid PNLFs additionally take into account the circumnebular extinction from Jacoby & Ciardullo (2025) as described in Sect. 2.6.

In the text

Fig. 10 Single stellar population PNLFs depending on the age and metallicity for the fiducial model with the nebular metallicity correction by Dopita et al. (1992) and circumnebular extinction by Jacoby & Ciardullo (2025), normalized to different quantities. The structure of the grid and the amount of PNe sampled for each PNLF are the same as in Fig. 7. The black lines are the PNLFs obtained from 106 stars of the respective final mass with equally spaced post-AGB ages between 0 and 30 000 yr. The red PNLFs are normalized to a total mass of 1012 M⊙ assuming a Chabrier IMF (Chabrier 2003). The orange PNLFs are normalized to a total bolometric luminosity of 1012 L⊙ based on the GALAXEV CB07 SSP models by Bruzual & Charlot (2003) using a Chabrier IMF.

In the text

	Fig. C.1 Nebular metallicity correction factor p5007 as a function of metallicity for the relations given by Jacoby (1989), Dopita et al. (1992), and Schönberner et al. (2010), as well as for having no metallicity dependence. The relations are set up to reach a value of p5007 = 1 at the base metallicity Zbase = 0.007.
In the text

	Fig. D.1 Post-AGB age tpost-AGB at which PNe at metallicity Z = 0.01 reach the brightest possible [O III] flux as a function of final mass Mfinal. The black data points mark the determined values from the model for 50 equally spaced final masses between 0.53 and 1.00M, and the red and orange lines are the linear fits to the values below and above Mfinal = 0.60 M⊙, respectively.
In the text

Fig. F.1 Single stellar population PNLFs depending on the age and metallicity for the fiducial model with the nebular metallicity correction by Dopita et al. (1992) and circumnebular extinction by Jacoby & Ciardullo (2025), normalized to different quantities. This grid of PNLFs shows the same plots as seen in Fig. 10, except that the λ5007 magnitude range is extended to dimmer PNe by four magnitudes. The black lines are the PNLFs obtained from 106 stars of the respective final mass with equally spaced post-AGB ages between 0 and 30 000 yr. The red PNLFs are normalized to a total mass of 1012 M⊙ assuming a Chabrier IMF (Chabrier 2003). The orange PNLFs are normalized to a total bolometric luminosity of 1012 L⊙ based on the GALAXEV CB07 SSP models by Bruzual & Charlot (2003) using a Chabrier IMF.

In the text