Tracing the Galactic Disk with Planetary Nebulae using Gaia DR3 Distance Catalog, Velocities, Populations, and Radial Metallicity Gradients of Galactic Planetary Nebulae

Aims. We want to study the population of Galactic planetary nebulae (PNe) and their central stars (CSs) through the analysis of their distances and Galactic distribution. The PN distances are obtained by means of a revised statistical distance scale, based on an astrometrically-defined sample of CSs from Gaia DR3 as calibrators. The new statistical distances, together with the proper motion of the CSs – also from DR3, and with published PN abundances as well as radial velocities, are used to characterize the PN populations in the Galaxy, and to derive the radial metallicity gradient. Methods. The statistical scale is applied to infer distances of a significant number ( ∼ 850) of Galactic PNe, for which we deliver a new catalog of PN distances. By adopting a circular velocity curve of the Galaxy, we also obtain peculiar 3D velocities for a large sample of PNe ( ∼ 300). Elemental abundances of the PNe are culled from the literature for an updated catalog, to be used in our analysis and other external applications. Results. The radial chemical gradient of the Galactic Disk is traced by PNe with available chemical abundances and distances, and kinematic data of the CSs are employed to identify the Halo PN population. We date PN progenitors based both on abundances and on kinematic properties, finding a confirmation of the first method with the second. For all PNe with at least one oxygen determination in the literature, we find a slope of the radial oxygen gradient equal to ∆ log(O / H) / ∆ R G = − 0 . 0144 ± 0 . 00385 [dex kpc − 1 ]. Furthermore, we estimate radial oxygen gradients for the PNe with old ( > 7.5 Gyr) and young ( < 1 Gyr) progenitors to be respectively ∆ log(O / H) / ∆ R G = − 0 . 0121 ± 0 . 00465 and -0.022 ± 0.00758 [dex kpc − 1 ], thus disclosing a mild steepening of the gradient since Galaxy formation, with a slope change of 0.01 dex. The time evolution is slightly higher ( ∼ 0 . 015 dex) if we select the best available abundances in the literature. This result is in broad agreement with previous PN results, but now based on DR3 Gaia analysis, and also in agreement with what traced by most other Galactic probes. We also find a moderate oxygen enrichment when comparing the PNe with young and old progenitors.


Introduction
The third data release of the Gaia ESA mission, DR3 (Gaia Collaboration 2016, 2023b) provides unprecedented astrometric parameters, complemented by photometric and spectroscopic data, for hundreds of central stars of planetary nebulae (CSPNe) that, further to witnessing a short-lived phase of their life, are suitable probes of chemical abundances in the Galaxy (e.g., Peimbert 1978;Henry et al. 2010), as well as in external galaxies (see Stanghellini 2019, and references therein).Metal contents of planetary nebulae (PNe) with young and old progenitors can be examined to constrain different mechanisms of gas accretion and star formation history (e.g., Gibson et al. 2013;Curti et al. 2020).
In particular, metallicity gradients along galactic disks are crucial observables in galaxy evolution studies.They are the outcome of complex interacting physical processes involving metal Full Tables 1, 3-5 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https:// cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/680/A104production, consumption, transport, and loss (Matteucci 2021;Prantzos et al. 2023).Their radial profiles reconstructed from different metal tracers currently show some diversities that need to be investigated further (Stanghellini et al. 2019).The basic observational step in this direction is to build a metallicity gradient as function of the galactocentric radius (R G ).The profile and slope of the gradient can be inspected for a coeval stellar population, and its time evolution can be tested with different age populations (e.g., Magrini et al. 2016;Stanghellini & Haywood 2018).To this end, being able to gauge the precise PN location as function of R G is clearly highly important.
In our previous paper (Stanghellini et al. 2020, hereafter Paper I), we performed a recalibration of the PN statistical distance scale from a sample of ≈100 CSPNe with parallaxes from Gaia Data Release 2 (DR2) with a relative uncertainty better than 20%, demonstrating the potential of Gaia astrometry in improving the extant knowledge on PN distances.Following this path, we present a revision of the PN distance scale using the best parallaxes available from Gaia DR3 (Lindegren et al. 2021b).In addition, the highly accurate Gaia DR3 proper motions are valuable parameters in the context of this study: Space velocities can be inferred by the CSPN proper motions using their distances.Because the velocity dispersion of a disk population tends to increase with stellar age (e.g., Strömberg 1925Strömberg , 1946;;Wielen 1977;Mackereth et al. 2019), the kinematic properties of CSPNe independently provide information about the progenitor age, which can complement and support the elemental abundance dating method.Therefore, in this work we have exploited the full astrometric data provided by Gaia DR3 for a bona fide set of Galactic PNe, in combination with their best available chemical abundances and radial velocities, to investigate the gas-phase metallicity gradient of the Milky Way and its evolution in time.
The structure of the paper is as follows.In Sect. 2 we characterize the Galactic PN sample and assess the new distance scale calibration based on DR3 parallaxes of their central stars and other nebular parameters.Section 3 describes the data and methods for the derivation of radial metallicity gradients.Finally, Sect. 4 is devoted to the discussion of the results and the conclusions.

Revised distance scale
We explored the statistical distance scale for Galactic PNe, which hinges on the commensurate relation between their surface brightness and physical radius advocated by Shklovsky (1956), which has been widely used and refined since then (see Smith 2015 for a thorough review).The similarity between the expansion radius of a fully ionized nebula and the decreasing density of its (nearly) constant ionized mass provides an estimation of the PN distance when its angular radius and observed flux have been measured.
Despite the wide range of PN mass progenitors, the socalled Shklovsky mass, that is, the apparent mass that would be derived from the measured flux according to simulated PN models (Buckley & Schneider 1995), has a fairly constant value of ≈0.1−0.2M .This means that the underlying assumption that all nebulae have the same ionized mass is viable.Nonetheless, the Shklovsky method tends to overestimate the distances of less evolved PNe that are in an ionization-bound state, that is, whose ionization front is still within the nebular material, far from the edge of the gaseous cloud (Schneider & Buckley 1996).
The wealth and unprecedented accuracy of Gaia DR3 parallaxes allow us to calibrate the physical radius-surface brightness distance scale on well-constrained empirical evidence while allowing us to test some of the crucial aspects of the underlying physical models.In logarithmic units, the scale is represented by the equation where the surface brightness S Hβ , a distance-independent parameter, is given in terms of the extinction-corrected absolute Hβ flux from the nebula I Hβ (in erg cm −2 s −1 ) and its angular radius θ (in arcsec) as S Hβ = I Hβ /πθ 2 , whereas the physical radius of the PN, expressed as R PN = θ /(206 265 p ) (in pc), depends on the parallax p of the CSPN, that is, its distance.

Identifying central stars of planetary nebulae in Gaia DR3
Starting from a list of 2556 true PNe extracted from the HASH database (Parker et al. 2016) 1 , we used the HASH coordinates of the nebular centers to searched for Gaia DR3 sources inside a radius of 3 arcsec and up to 5 arcsec if the smaller radius returned no matches.We resolved multiple matches by applying a ranking system based on the astrometric quality indicators provided with the Gaia solution.In this way, we identified 1674 CSPN candidates lying near the geometric center of the nebula with a measured DR3 parallax.Furthermore, to verify the reliability of our purely astrometric criterion, we made comparisons with the CSPN catalogs that were recently published by González-Santamaría et al. (2021) and Chornay & Walton (2021), whose matching algorithms relied on the Gaia BP-RP color to distinguish between hot and cool star candidates and weight them accordingly.Because our main aim in this context was sample purity and not completeness, we retained only objects for which our matched Gaia source agreed with the González-Santamaría et al. ( 2021) and Chornay & Walton (2021) catalogs.We obtained a sample of 942 CSPNe, and this constitutes our calibrator sample.

Selection of the distance-scale calibrators
To determine the parameters a and b of Eq. ( 1), we matched the calibrator sample with the physical parameters that are needed to build the scale.Table 1 contains the calibrator sample PNe for which angular radius, Hβ flux, extinction correction, and an estimation of the ionized mass are available.For each of these targets, we list the PN G name, the Gaia DR3 ID name, which is different from the DR2 ID, the common name used for the PN, the corrected (see below) DR3 parallax c , the angular radius θ, the logarithmic Hβ flux, the logarithmic extinction constant, and the ionized mass logM i , the latter calculated following the prescription described in Paper I. When no uncertainties on the physical parameters were available, we assumed a typical error for each parameter.In particular, we adopted a 20% uncertainty on the angular radius, and 0.1 as the uncertainty on the logarithmic extinction.
All parameters, excluding those derived from Gaia DR3, were selected from the literature (references in Stanghellini & Haywood 2018).The parallaxes in Table 1 are those used to fit the scale and were corrected by applying a global offset of 0.017 mas to the DR3 parallaxes, plus an additional bias term that depends on the stellar G magnitude, its color G BP − G RP (or pseudo-color 2 , for a six-parameter (6p) astrometric solution), and ecliptic latitude (see details in Lindegren et al. 2021a).Moreover, the parallax errors were multiplied by 1.05/1.22for a 5p/6p parameter solution to compensate for an underestimation (bias) of the error as described in Fabricius et al. (2021).In the following, the symbols c and σ c indicate that the corresponding DR3 parallax and standard error σ were corrected for these systematic effects.
When we matched our calibrator sample with the PNe whose radius, Hβ flux, extinction, and ionized mass were also available, we obtained the 401 calibrators given in Table 1.For σ c / c < 20% (10%), we still have 137 (74) calibrators, which doubles the corresponding samples in the DR2 calibration of Paper I.

Building the distance scale for planetary nebulae
The semi-empirical relation expressed by Eq. (1) sets the stage for the inference problem relating the sought-for posterior density distribution of the unknown parameters a and b to the likelihood of the observed physical quantities characterizing our PN calibrators.We briefly describe the method here and refer to the appendix of Paper I for a complete derivation.
To formulate our likelihood function, we first introduce the measured radius θ ∼ N(φ, σ θ ), flux I ∼ N(J, σ I ), and parallax ∼ N(p, σ ) as Gaussian-distributed uncorrelated variables, then marginalize the likelihood of the ith set of measurements with respect to the parameters radius (φ), flux (J), and parallax (p) and treat the parallax as a dependent variable through Eq. (1), obtaining the expression of the likelihood for the ith calibrator as where we used uniform priors for φ and J inside some reasonable intervals.Then, after forming the log-likelihood of the complete set of calibrators, by virtue of Bayes theorem, we derive the bivariate posterior probability density p pdf (a, b) assigning uniform priors to the a and b parameters' space as We approach this problem by numerical integration over the range of the φ and J parameters and by sampling the posterior distribution of a and b inside a bidimensional grid with step sizes of 0.0005 and 0.005 for the slope and intercept, respectively.The estimated â and b correspond to the mode of the discrete posterior probability density function.The 68% confidence intervals for the two parameters are calculated separately by projecting onto the a and b axes the isoprobability contour of the posterior distribution for which |p(a, b) max − p(a, b)| ≤ 0.5 (recalling that the log-likelihood ratio statistics is asymptotically distributed as χ 2 /2).Finally, we estimate the correlation term from the complete posterior sample as obtaining in all cases a highly positive correlation between the estimated slope and intercept, as illustrated in Fig. 1 by the density levels of the bivariate posterior distribution p pdf (a, b) for the σ c / c < 20% sample-cut case.
By adopting different thresholds on the relative DR3 parallax error σ c / c , we inferred different scale parameters, giving consistent estimates of the scale slope a and intercept b, as reported in Table 2. Figure 2 shows the log R−log S Hβ plane with the subsets of calibrators for which σ / ≤ 10% (top panel) and σ c / c ≤ 20% (bottom panel) along with their error bars, and the corresponding estimated linear relation (see caption).
An inspection of these plots shows that both data sets present a somewhat large dispersion around the scale.Moreover, from a comparison of the 20% sample of Fig. 2 with that of Fig. 2 in Paper I (where the calibration used DR2 parallaxes), we note that the numerous stragglers populating the lower right end of the plane, which correspond to PN with very low ionized mass (log M i < −2), are not so prominent in the new plots.
In the complete set of potential calibrators, 37 objects had log M i < −2, but the parallaxes of none of them were accurate enough at the 10% level, and only a few (indicated in red in the lower panel of Fig. 2) met the 20% requirement.When we placed these 37 objects on the log R PN −log S Hβ plane, they were scattered mainly in the lower part of the graph, and their distance did not conform with the statistical scale.This corroborates the findings of our previous paper based on DR2 data.
To further investigate the ionized-mass effect, we inspected the PNe with σ c / c < 20% by displaying their log M i on a color-density scale.The result, shown in Fig. 3, is quite striking and can be fully appreciated through the quality of the calibrators: The dispersion around the linear scale is clearly linked to the (estimated) log M i of the PN shell, manifesting the deviation of its actual evolutionary state.In the transition from ionizationbound, or optically thick, to density-bound, or optically thin, the ionized mass of PN increases, as discussed in Pottasch (1980).Therefore, the model assumption of fully ionized, constant mass, on which Eq. (1) stands, does not strictly hold.
A104, page 3 of 10 Remarkably, when we tried to tighten the constraint on the goodness of the parallax relative error even more, we obtained a similar dispersion of the calibrators as in Fig. 3. Finally, by subtracting the effect due to the nominal uncertainties of the observed angular radius and flux, we tentatively estimated an average intrinsic (cosmic) scatter of ∼0.1 dex around the distance scale, and we conclude that this is due to an actual observed nebular evolution that limits the potential accuracy of this statistical distance scale, unless the ionized mass of the PN is known a priori.
The results of the calibration are reported in Table 2, where for each of the three relative uncertainty thresholds on σ c / c , we give the number of calibrators N cal and the estimated slope and intercept along with their correlation ρ ab .We also list the averaged parameter K = statistical_scale_distance × c , which ideally, is equal to 1, and its dispersion σ K , which is an indicator of the goodness of the scale (see Smith 2015).
We used the third scale listed in Table 2, which is calibrated with central stars whose DR3 parallaxes were measured with an uncertainty better than 20%, which has a K ∼ 0.964.Through the analysis presented in this paper, we have confirmed that using another of the three scales would not affect the science results.Our adopted scale is log R PN = (−0.242± 0.0042) × log S Hβ − (4.2 ± 0.057). (4) In Table 3 we list the new catalog of distances of 843 Galactic PNe, based on Eq. ( 4).For each PN in the HASH catalog whose Hβ intensity and apparent radius is available, we calculated its heliocentric distance, D, based on our scale.
In the catalog, we list the PN G name, the angular radius θ, the logarithmic Hβ flux and its uncertainty, the logarithmic extinction constant and its uncertainty, the heliocentric distances from our scale of Eq. (4), D [kpc], and their left and right formal uncertainties.The latter were computed by estimating σ log D via first-order error propagation of the linear relations among the variables log R, log S Hβ θ, and D, accounting for uncertainties in both observed and fitted parameters; then, asymmetric uncertainties were defined as σ D − = 10 (log D−σ log D ) and σ D + = 10 (log D+σ log D ) .Finally, we included the PN peculiar velocity V pec and its uncertainty calculated from Eq. ( 12), plus a disk, bulge, and halo population classification, as discussed in the following section.

Planetary nebulae as tracers of the Galactic disk
An important scientific motivation for the Galactic PN distance scale is the study of the radial (and vertical) metallicity gradients in the Galaxy.Planetary nebulae are the ejecta of evolved 1-8 M stars, thus they represent a continuum of stellar populations from very old through relatively young, with ages ∼0 to ∼10 Gyr, which makes them very versatile probes for Galactic evolution.Gas-phase abundances of α-elements observed through emission lines and measured in their shell should be the same as those of their progenitors at the time of formation because these elements do not change much during the evolution of these stars.Thus, tracing the abundances of these elements is equivalent to probing the Galaxy abundance with a look-back time from ∼0 to ∼10 Gyr, depending on the initial progenitor mass of the PNe.

Elemental abundances
We used PN elemental abundances both for progenitordating purposes and to determine the metallicity gradients.We selected the elemental abundances from all references included in Stanghellini & Haywood (2018), to which we added the abundances from the data sets by Tsamis et al. (2003), Miller et al. (2019), McNabb et al. (2016), and Delgado-Inglada et al. (2015), which were not included therein.In the references where the abundances are given both for the whole nebula and for different parts of the PN, we used only the measures that were based on spectra of the whole PN, so that the abundances from different references are readily comparable with one another.The final abundances were curated, that is, we recalculated the ionization correction factor (ICF) uniformly, as described in Stanghellini & Haywood (2018).
In Table 4 we list the average abundances used in this paper.Column (1) gives the PN G name, and in Cols.(2) through (4) we give the C, N, and O abundances in the usual scale log(X/H) + 12, averaged from all the available references.The uncertainties, given in dex, are the dispersion of that elemental abundance across the literature, estimated with σ(log(O/H)) = 0.434 × [σ(O/H)/ O/H ].We also produced a curated catalog of selected PN abundances, shown in Table 5.For this selection, we prioritized abundances from space-based observations unless they were deemed uncertain in the original papers; otherwise, we used the published ground-based abundances, prioritizing the most recent, or recently revised, abundances.Curated abundances are given with their uncertainty when available in the original reference.

Velocities of the central stars of PNe
For a sizable sample of CSPNe with a statistical distance from our catalog, Gaia DR3 provides equatorial proper motions µ α * ± σ µ α *3 , µ δ ±σ µ δ .We used this kinematic information as insight for their populations and ages.After converting DR3 proper motions and their standard deviations into galactic coordinates µ l ± σ µ l , µ b ± σ µ b via the transformation matrices A G and C Gal defined by Eqs.(4.62) and (4.79) of the Gaia EDR3 online documentation 4 , we computed the corresponding observed spatial velocities plus errors, given in [km s −1 ], as where l and b are the Galactic longitude and latitude.We also adopted σ D ≡ (σ D − + σ D + )/2.The third component of the velocity is the radial component, V r , which has been observed independently for 498 PNe of our sample (Durand et al. 1998) 5 .
Having secured distances and kinematic data for these stars, we can estimate their peculiar velocity V pec as the residual stellar motion that does not conform to the general Galactic rotation.To this end, we assumed an empirical model for the Galactic rotation curve that follows Eilers et al. (2019) for galactocentric distances larger than 5 kpc, while for inner distances, we made use of the results of Bhattacharjee et al. (2014), who resorted to H I and H II regions, as well as CO emission lines, as tracers of the inner Galactic disk.We adapted the results of Bhattacharjee et al. (2014) to the fundamental parameters found by Eilers et al. (2019), that is, R = 8.122 kpc, the galactocentric distance of the Sun, and V LSR = 229 [km s −1 ], the circular velocity of the local standard of rest (LSR).We combined these data into a grid of Galactic radii and corresponding rotation velocities, plus uncertainties, as given in Table 6.Moreover, in the following calculations, the velocity of the Sun with respect to the LSR was set to (U , V , W ) = (11.1,12.4, 7.25) [km s −1 ], according to Schönrich et al. (2010).
We first transformed the observed spatial velocities into Cartesian components along the local triad at the star as To derive galactocentric velocities, we must account for the galactocentric velocity of the Sun, that is, u = (U , V + V LSR , W ), so that For the following analysis, it is useful to express spatial velocities in cylindrical coordinates, for which we need to compute the galactocentric azimuth φ as (e.g., López-Corredoira & Sylos Labini 2019)  where D is the heliocentric distance of the star, and φ is taken as positive in the direction of galactic rotation; the radial, azimuthal, and vertical components of the spatial velocity are then given by Finally, we are able to estimate the peculiar velocity of each star by subtracting the circular motion due to differential galactic rotation from the azimuthal component V φ .To do this, we computed for each PN the galactocentric radius R G projected onto the Galactic plane as and we interpolated linearly between the nearest two values of R G of Table 6 to calculate the circular stellar velocity V c (R G ).
The absolute space motion of the stars (or the PNe), which we call V pec , can be obtained by quadrature of the three velocity components as where σ V pec is a lower-limit uncertainty, which assumes uncorrelated variables and a perfect model.The peculiar velocities for the PNe are reported in Table 3.
In Fig. 4 we placed the inferred peculiar velocities in the socalled Toomre diagram, which can be used to identify regions in which the probability of finding disk versus halo stars is high.In particular, following Bonaca et al. (2017), we used V pec > 220 [km s −1 ] as the threshold value for rejecting halo stars.

Galactic PN populations, and dating PN progenitors
We derived heliocentric distances D and their confidence limits, as described in Sect.2.3, for 843 Galactic PNe.The range of distances spanned by these PNe is ∼150 to ∼27 000 [pc], and most of them can be included in a general disk population, thus they are adequate for our science scope.We populate the halo PN sample selecting them by peculiar velocity, where halo PNe have V pec > 220 [km s −1 ].This selection yields 13 halo PNe.We also tentatively defined the bulge population similarly to Stanghellini & Haywood (2018) by selecting the PNe in the central 3 × 3 [kpc], that is, with |z| < 3 and 0 < R G < 3 [kpc].To constrain the bulge sample even more, we selected in this sample only the PNe whose apparent radius was θ < 5 .Our selection gives 102 bulge PNe.Whether a PN belongs to the halo or bulge populations according to our selection is noted in Table 3.
In Fig. 5 we plot the distribution of the peculiar velocity versus vertical height z, where the color indicates the distance from the Galactic center.In the area between the two gray lines (±0.175 kpc) PNe with high peculiar velocity are away from the Galactic center.It is worth noting that all bulge PNe selected as in Sect.3.3 are outside the gray lines and have an intermediate V pec .
A104, page 6 of 10 Dating PN progenitors is not trivial, and it has been attempted by Peimbert and Torres-Peimbert by dividing Galactic PNe into populations (i.e., PN types I through V).This approach uses both the chemistry of the PNe and their kinematics.Stanghellini & Haywood (2018) used the final yields of the evolutionary elements compared to nebular abundances of carbon and nitrogen to select two PN progenitor populations with extreme ages, called OPPNe (old progenitors PNe, with progenitor ages > 7.5 Gyr), and YPPNe (young progenitors PNe, with progenitor ages < 1 Gyr).This selection is based on the observation that PNe with old and young progenitors occupy markedly different loci on the N/H versus C/H, and N/H versus O/H, planes.4 and 5, Col. 5, we list the PN progenitor age group as derived above, when available.
More insight into dating PN progenitors is derived from the comparison of Gaia proper motions described in the previous subsection, with the ages derived purely from the PN chemistry above.In Fig. 6 we show the resulting distribution of V pec for the general PN population (which includes OPPNe, YPPNe, and the intermediate-age progenitor PNe, and the PNe whose age from chemistry was unavailable), then by plotting separately OPPNe and YPPNe.We find that the OPPNe distribution peaks at V pec ∼ 60 [km s −1 ], while the YPPNe distribution peaks at a lower velocity (∼25−30 [km s −1 ]).This plot shows the correla- Fig. 5. Galactic distribution of V pec from Eq. ( 12), with Galactic altitude, where the color is coded for the distance from the Galactic center R G .The vertical thin lines encompass the |z| < 0.175 zone.tion between completely independent dating systems from the chemistry and the kinematics of the PNe, giving us confidence about using the OPPNe and YPPNe distributions.

Galactic metallicity gradients and other metallicity variations
The radial metallicity gradients in this paper were calculated by orthogonal distance regression, which makes use of the uncertainties in the measured quantities R G and log(O/H) + 12.
The uncertainties in the elemental abundances are described in Sect.3.1, where we indicate the literature fonts from which we chose the abundances and how the errors were estimated both in the average and chosen abundances (see Table 4).The uncertainties in R G were computed by first-order propagation of the symmetrized error of the distance scale, The regression yields a slope and an intercept, along with their formal errors.
In Table 7 we present the estimates of the gradient slopes obtained using different subsets of PNe, and with both averaged A104, page 7 of 10  7.
and selected chemical abundances.We used all PNe for which at least one oxygen measurement was available in the literature, excluding only the extreme outliers (log(O/H) + 12 < 7.5 and log(O/H) + 12 > 9).All radial oxygen gradients that we derive in this paper have negative slopes.The general PN population (Fig. 7) has a gradient slope of −0.0144 ± 0.00385 [dex kpc −1 ] when we use the average abundances.When the halo population (3 PNe) is excluded, the gradient does not change significantly.For all comparable samples, the gradients steepen with time since Galaxy formation.For the OPPN and YPPN populations, we derived gradient slopes ∆log(O/H)/∆R G = −0.0121± 0.00465 and −0.022 ± 0.00758 [dex kpc −1 ], respectively, when we used average abundances, and ∆log(O/H)/∆R G = −0.0132± 0.00464 and −0.0278 ± 0.00789 [dex kpc −1 ], respectively, when we used the selected abundances (plots in Fig. 8).This result indicates a mild evolution, that is significant when compared to the uncertainties, of the metallicity gradient with time.Our results point to a steeper gradient slope for the younger stellar population than for the older stellar population, with a slope difference of ∼0.15 dex.This implies a slow time evolution of the gradient.The time evolution of the gradient slope does not vary,  7. Notes.Gradients are estimated using chemical abundances from Tables 4 (average abundances) and 5 (selected abundances), and for the listed PN samples, of numerosity N.
regardless of whether we include halo PNe in the OPPN sample (see Table 7).
It is beyond the scope of this paper to discuss the causes for the steepening gradient.Chemical evolutionary models of star-forming galaxies (Gibson et al. 2013) indicate that the radial metallicity gradients steepen with time when the chemical evolution proceeds with feedback of fresh, metal-poor gas from the environment.Models like this cannot predict the gradient ab initio.The metallicity gradients for old populations such as those from OPPNe therefore represent a valuable constraint.
Finally, we tested the oxygen evolution (or chemical enrichment) by comparing average oxygen abundances in the OPPNe and YPPNe samples to derive an average Galactic enrichment of ∼0.06 dex between the age > 7.5 Gyr and the age < 1 Gyr populations.

Results and conclusions
This paper yields essentially three novel results for PNe in the Galactic environment.First, we produced a new distance scale based on DR3 parallaxes, which is given in Eq. ( 4).The distances derived from the new scale are published as a catalog in this paper (Table 3).By comparing our new distances with the CSPN DR3 parallaxes, as in the analysis by Smith (2015), A104, page 8 of 10 we find K = 0.964, where K = D× c , and σ K = 0.154 when we exclude the PNe with logM i < −2 from the comparison, and K = 0.931 when we include all calibrators.Comparing these results with those of Paper I, where the best scale obtained with DR2 parallaxes gave K = 0.948 and σ K = 0.25, we conclude that the use of DR3 parallaxes as calibrators successfully improved upon the Galactic PN distance scale.
We have calculated K based on DR3 parallaxes for the most commonly used PN distance scales, such as Frew et al. (2016, hereafter FPB) and Stanghellini et al. (2008, hereafter SSV).We found K FPB = D FPB × c = 1.272 and K SSV = D SSV × c = 1.203 with the exclusion of calibrators with logM i < −2, and K FPB = 1.36 and < K SSV = 1.29 when all calibrators are included.We conclude that the scale presented in Eq. ( 4) gives the most reliable distances to date, and the distances in Table 3 should be used, unless an accurate parallax measurement of the central star of PN under study is available, from DR3 or otherwise.Moreover, our analysis revealed that the logR PN −logS Hβ distance scale of Galactic PNe has an intrinsic mean dispersion of ≈ 0.1 dex around the log R PN − log S Hβ linear relation, which can be accounted for by the variation in the PN ionized mass (excluding objects with log M i < −2), bringing us to the conclusion that statistical distance scales such as this can hardly be improved by increasing the accuracy of calibrators.
Second, we used the Gaia DR3 proper motion for thorough description of the peculiar velocities of the CSPNe, and to characterize the halo population.The peculiar velocity derivation has the additional advantage to be a completely independent assessment of the progenitor ages of the PNe, thus validating even further the OPPN and YPPN classes based on elemental abundance ratios C/O and N/O.The proper motion and peculiar velocity analysis of Sect.3.2 provided a tool for excluding halo PNe from the gradient analysis.Even more importantly, it gave us additional confidence in the progenitor dating from chemical analysis, as described in Sect.3.3.To our knowledge, this is the first time that dating PNe is consistent with two completely independent approaches, while in the past, the dating scheme has used either approach, or both approaches for different age classes, such as in defining Peimbert's PN types (Peimbert 1978).Dating PN progenitors is essential to study the time evolution of the radial metallicity gradients in the Milky Way.
Third, we measured the radial oxygen gradients based on DR3 distance scale, with halo PNe selected through their kinematic properties measured by Gaia DR3, and progenitor ages from chemical evolutionary assessment, confirmed via the peculiar velocities determined by Gaia.The negative shallow gradients that we determined and the mild evolution of the gradients, which steepened since the formation of the Galaxy, are a confirmation of past results both in the Galaxy (Stanghellini & Haywood 2018) and in the nearby spiral galaxies examined (Peña & Flores-Durán 2019;Magrini et al. 2016), where all emission-line probes in star-forming galaxies (including the Milky Way) indicate that the oxygen radial gradients are steeper in the younger populations.For the gradients, we used curated catalogs of published abundances, which we give in Tables 4 and 5.
It is worth noting that while our results confirm most of the metallicity gradients measurements based on gas-phase data published to date, recent cluster and stellar results show a more controversial picture (Anders et al. 2017;Magrini et al. 2023;Gaia Collaboration 2023a).Stellar metallicity are mainly inferred from iron rather than oxygen.Because iron atoms in PNe are mostly found in condensed state (e.g., Delgado-Inglada et al. 2015), a direct comparison between the stellar and nebular samples is typically problematic.
The work of Magrini et al. (2023) also included α-element gradients, indicating very little oxygen gradient evolution, compatible with flattening of gradients with time.Their oxygen gradient slopes of the age < 1 Gyr and age > 3 Gyr bins in their Table A.10. are the same within the uncertainties.We underline that these results are compatible with ours if we take into account the different age ranges used: our PN progenitors are in the age < 1 Gyr and age > 7.5 Gyr bins, whereas their oldest clusters are within 7 Gyr.Furthermore, migration and evolution due to inflow of clusters and stars may be different in ways that are unpredictable when measuring radial gradients.
Interestingly, Bhattacharya et al. (2022), who studied the radial oxygen gradients in M 31, found that young thin-disk PNe describe a steeper radial gradient than the older thick-disk population.This agrees qualitatively with our results, where radial metallicity gradients steepen with galaxy evolution.Quantitatively, they attempted to unify the gradients slopes across galaxies in their Fig. 13 by taking disk scale lengths into account, and they showed remarkable similarities between oxygen gradients in the Galaxy and M 31.
More inspection of Galactic gradients from data and modeling comes from the recent analysis in Lian et al. (2023), who used integrated radial metallicities to find (see their Fig. 1) that radial Fe/H gradients from older populations are flatter than for the young population.This agrees with our results.They also compared the Galactic gradient from present-day stellar population with that obtained from gas-phase oxygen (H II region), which is consistent with the PNe abundance tracers.This reinforces the confidence in our findings.
The distance scale that we presented here and the distances in our catalog will be used in the future to study the behavior of dust and other elements in Galactic PNe in the framework of the chemical evolution of the Galaxy.

Fig. 1 .
Fig. 1.Normalized 2D posterior probability of the distance-scale parameters obtained from the sample of PN calibrators with σ c / c < 20%.

Fig. 2 .
Fig.2.Log(R)[pc]  vs. log(S Hβ ) [erg cm −2 s −1 ] plot for calibrators with σ c / c < 0.1 (top panel) and σ c / c < 0.2 (bottom panel).In both plots, the blue points are the calibrators.Their logarithmic (asymmetric) error bars reflect the 1σ uncertainties of the observed parameters.In the bottom plot, the red points are the calibrators with log M i < −2.The red line shows the two distance scales from the calibrators, and in cyan, we plot the 2σ confidence band of the scales, representing the uncertainty of the logarithmic radius R induced by the correlated errors on â and b.

Fig. 3 .
Fig. 3. Log(R) [pc] vs. log(S Hβ ) [erg cm −2 s −1 ] plot for the PN sample with σ c / c < 0.2, color-coded according to their ionized mass.The error bars have been omitted for clarity.The cross in the bottom left corner is representative of the typical error on the respective axis.
7) where e T l = (− sin l, cos l, 0), e T b = (− sin b cos l, − sin b sin l, cos b), e T r = (cos b cos l, cos b sin l, sin b), and A is the so-called direction cosine matrix representing the rotation to the local frame.

Fig. 4 .
Fig. 4. Halo PNe selected through the velocity analysis of the CSPNe by means of the Toomre diagram.The color bar is coded for the peculiar velocities as in Eq. (12).The dots outside the region encircled by the red curve correspond to halo PNe.

Fig. 6 .Fig. 7 .
Fig.6.Histogram of the spatial peculiar velocity V pec [km s −1 ] for the general disk population and for the OPPNe and YPPNe separately.The velocity bin is 10 km s −1 .The y-axis is in logarithmic scale.

Fig. 8 .
Fig.8.Radial oxygen gradients for the OPPN (top) and the YPPN (bottom) populations.All PNe have been included in these plots, and halo and bulge populations are indicated with different colors.Oxygen abundances are those selected here.The red lines are the linear fits of the data.The slopes and intercepts for the different populations are listed in Table7.

Table 1 .
PN calibrator parameters.Notes.See the main text for a detailed column description.The full table is available electronically at the CDS.

Table 2 .
Summary of the distance-scale analysis.Notes.N cal = number of calibrators; N comp = number of PNe used to compute K and σ K : in group (1), we used the calibrator sample with the additional requirement of log M i > −2, and in group (2), we included all PNe with a distance and a parallax determination, and again log M i > −2.

Table 3 .
Catalog of statistical distances and peculiar velocities of Galactic PNe.Distances are obtained from scale 3 of Table2; see the main text for a detailed column description.The full table is available at the CDS.

Table 5 .
Selected elemental abundances of Galactic PNe.
PN G log(C/H) + 12 log(N/H) + 12 log(O/H) + 12 Prog.Notes.In the computation of abundances, space-based observations and most recent ground-based abundances have been prioritized (see main text) The full table is available at the CDS.