© The Authors 2026

1 Introduction

Planetary nebulae (PNe) are short-lived ionized shells ejected by low- and intermediate-mass stars as they evolve from asymptotic giant branch (AGB) to white dwarfs (Shklovsky 1957; Jacob et al. 2013). As key tracers of stellar evolution, chemical enrichment, and the interstellar medium (Kwitter & Henry 2022), accurate distances to PNe are crucial for determining their intrinsic properties (Gathier et al. 1986; Pottasch 1996; Parker 2022). However, obtaining reliable PNe distances remains a long-standing challenge, as existing methods struggle to achieve both wide applicability and high measurement reliability.

Distances to PNe are commonly derived through individual or statistical approaches. The individual approach measures the distance to a single PN via trigonometric or expansion parallaxes, photoionization modelling, kinematics, cluster membership, or extinction. Although generally regarded as the primary method owing to its higher reliability, it applies to only very few PNe under specific observational or physical conditions (Phillips & Márquez-Lugo 2011; Kwitter & Henry 2022; Parker 2022). The statistical approach infers distances indirectly for large PNe samples from assumed physical properties and empirical relations. While well-suited to population studies, the statistical method introduces systematic uncertainties from calibration datasets, methodologies, and assumptions about PNe evolution and is commonly regarded as a secondary method in the literature (Smith 2015; Kwitter & Henry 2022).

The high-precision Gaia mission has greatly advanced PNe studies. About 70% of the ~3800 known Galactic PNe (Frew et al. 2016, hereafter F16) now have identified central stars of planetary nebulae (CSPNe; Chornay & Walton 2021, hereafter CW21; González-Santamaría et al. 2021, hereafter GS21). Gaia parallaxes for these CSPNe provide the most direct and reliable distance determinations (Parker 2022; Kwitter & Henry 2022). Nevertheless, more than 75% of CSPNe are too faint for precise parallaxes, and misidentifications remain non-negligible, limiting the overall utility of CSPNe-based distances. To expand distance coverage, several recent studies have combined CSPN parallaxes with statistical scales, producing samples of 96 PNe (Ali et al. 2022, hereafter A22), 843 PNe (Bucciarelli & Stanghellini 2023, hereafter BS23) and 2211 PNe (Hernández-Juárez et al. 2024, hereafter HJ24). The largest of these, HJ24, compiled distances from CSPN parallaxes and multiple statistical scales, adopting a rule-based scheme to select the most reliable value for each PN. The associated uncertainties were assigned empirically from comparisons among different distance estimates, ranging from a minimum of 38% to typical values around 70% in the published catalogue. Large discrepancies between distances derived from different methods further highlight the difficulty of establishing a reliable distance scale. For instance, the reported distance to Abell 38 ranges from 474 pc to 6201 pc (F16; CW21; GS21; BS23; HJ24), a spread of nearly a factor of 13. Even the most reliable CSPNe-based distances can disagree substantially. For example, the distance to PN HaTr 14 is reported as 12124 pc and 4954 pc (CW21; GS21), differing by roughly a factor of 3. This illustrates that substantial discrepancies can persist even among results from the most trusted method, most likely due to CSPN misidentification. Thus, even with many identified CSPNe, reliable distances still remain challenging, calling for improved methods.

Among individual methods, the extinction–distance approach is physically well-motivated and independent of assumptions about nebular evolution. Traditionally, PN extinction is inferred from emission-line or radio flux ratios, assuming a fixed total-to-selective extinction optical ratio (R_V ≈ 2.5–3.2; e.g. Acker 1978; Saurer 1995), with typical values of A_V ~ 0.1–8.5 mag (Dharmawardena et al. 2021, hereafter D21). Distances are then estimated either by comparison with field stars along the same sightline (Gathier et al. 1986) or by matching to 3D extinction maps (D21). Despite decades of effort, this method has produced distances for only a few dozen PNe (6–70 objects; Lutz 1973; Acker 1978; Martin 1994; Saurer 1995; Giammanco et al. 2011, hereafter G11; D21). The main limitations of this method arise from incomplete or uncertain PN extinction estimates, which result from restricted Balmer coverage, missing or low-quality radio fluxes, uncorrected intrinsic colours, and the use of a constant R_V that in reality varies significantly across sightlines (D21, Wang & Jiang 2014; Wang & Chen 2019, 2023; Green et al. 2025). Its applicability is further constrained by the small number of comparison stars and by the coarse angular resolution of current extinction maps.

Recent advances offer a more robust foundation for extinction-based distances. The blue-edge method, used to estimate stellar intrinsic colours, has been continuously refined (Wang & Jiang 2014; Wang & Chen 2019; Zhao et al. 2024; Deng et al. 2025), yielding accurate stellar extinction without assuming a fixed R_V. The extinction-jump model enables distance estimates via extinction jumps and has been successfully applied to supernova remnants and molecular clouds (Wang et al. 2020; Cao & Jiang 2024; Chen et al. 2025). When combined with Gaia’s all-sky coverage and high-precision parallaxes, which supply large and reliable stellar samples, these developments enable substantially improved extinction–distance determinations for PNe.

In this work, we combined the refined blue-edge method with the extinction-jump model to develop a Gaia-based extinction–distance method for PNe with identified CSPNe. The resulting distances are constrained jointly by multiple indicators, including CSPN-based estimates, parallaxes, and statistical scales, and are further validated using 3D extinction maps and the empirical R_PN–distance relation. This method yields robust distances for over one thousand PNe, addresses the difficulties posed by missing or unreliable CSPNe parallaxes, and helps identify more reliable distances when published estimates differ widely. This paper is organized as follows. Section 2 describes the sample and data. Section 3 presents the methodology. Section 4 discusses validation. Section 5 provides the results and discussion, and Section 6 summarises our conclusions.

2 Sample and data

2.1 PN sample

To address the shortage of reliable distances for PNe with identified CSPNe, we constructed our working sample from the two largest available CSPNe catalogues, CW21 and GS21. Both catalogues are based on the HASH (Hong Kong/AAO/Strasbourg/Hα) PN database (Parker et al. 2016), which includes nearly all known Galactic PNe from previous compilations and provides CSPNe identifications using Gaia EDR3. CW21 compiled a catalogue of 2117 CSPNe and derived distances for 733 PNe by combining Gaia parallaxes with the widely used statistical distance catalogue of F16. GS21 reported 2035 CSPNe, among which 405 distances from Bailer-Jones et al. (2021, hereafter BJ21) were classified as reliable and recommended for use. Combining both catalogues yields 2570 Galactic PNe including all CSPNe, with 897 PNe having reliable distances.

For these 2570 PNe, we retained all reported distances and extracted their central coordinates and angular radii (R_PN) from the HASH database. Each PN was assigned an individually defined circular field for extracting stars to construct a reliable extinction–distance profile (see Section 2.3). Planetary nebulae lacking a defined R_PN or containing fewer than ten stars within the field were excluded, leading to the final sample of 2225 PNe.

2.2 Data

For each PN, we selected stellar samples along the line of sight to trace the cumulative extinction towards it. Gaia DR3 provides high-precision photometry in the G_BP, G_RP, and G bands, together with spectroscopic and astrometric data for over 1.8 billion sources (Gaia Collaboration 2023). From its lowresolution XP spectra, Andrae et al. (2023) derived robust stellar atmospheric parameters for ~175 million stars, including effective temperature (T_eff), surface gravity (log g), and metallicity ([M/H]), with typical uncertainties of ~50 K, 0.08 dex, and 0.1 dex, respectively. For this study, we extracted these stellar parameters from Andrae et al. (2023), as they are required for intrinsic-colour estimation using the blue-edge method. Stellar distances were taken from BJ21, which provides geometric distances from Gaia parallaxes and photogeometric distances incorporating photometry and astrophysical priors, yielding improved accuracy for stars with low-quality parallaxes. Distances for CSPNe were adopted from either the geometric or photogeometric solution depending on parallax quality, whereas non-CSPNe used only geometric distances. We excluded stars that satisfied any of the following criteria:

• negative or missing parallaxes,

• fractional parallax error >20%,

• G > 17.65 mag,

• T_eff outside the range 4000–8000 K,

• uncertainty in log g or [M/H] > 0.3,

• T_eff uncertainty >0.1 (BJ21, Andrae et al. 2023; Zhao et al. 2024).

Combined with the depth and precision of Gaia DR3, these selection criteria ensured a clean and reliable stellar sample, providing a solid foundation for constructing extinction–distance profiles.

Fig. 1 Left: stars overlaid on the Digitized Sky Survey image with RPN (dashed yellow circle), the CSPN (yellow star), as well as core, nebular, and peripheral sources (red, dark blue, and blue dots). Right: extinction vs distance.

2.3 PN stellar field

The extinction–distance method requires a stellar sample for each PN to build an extinction–distance profile along its line of sight. Typically, a stellar field is defined around each PN to select nearby stars. The field size must balance statistical robustness with sensitivity to extinction jumps. Previous studies typically adopted constant stellar fields on angular scales ranging (R_field) from 1.5° to 10′ (Lutz 1973; G11) to ensure sufficiently large stellar samples. Although simple, this strategy ignores PN size variations (typically ~1′−2.5′) and introduces background contaminants, which reduce spatial sensitivity and limit the applicability of the method.

Using Gaia DR3, we defined R_field for each PN with a PN-size–dependent approach guided by spatial-consistency checks to construct reliable stellar samples and minimise contamination from nearby regions. The planetary nebula angular radii (R_PN) were taken from the HASH database (Parker et al. 2016). Stars within each R_field were then divided into three concentric zones based on spatial distribution. The core zone includes stars within 0.5 R_PN, the nebular zone covers 0.5–1 R_PN, and the peripheral zone extends from R_PN to R_field. This division enabled us to check the spatial consistency of extinction–distance trends and validating the chosen R_field.

As shown for Sh2-123 (Figure 1), the extinction–distance relations across these zones are consistent, confirming the reliability of our R_field selection. Systematic tests across the full sample yield the optimised selection criteria summarised in Table 1. For compact PNe, R_field was extended to several times R_PN. These criteria were designed to maximise completeness while minimising contamination.

3 Method

We adopted a three-step approach to estimate PN distances based on extinction–distance profiles. Extinction–distance profiles were first constructed for each PN sightline. Significant extinction jumps were then identified in these profiles through model fitting or by directly tracing features near the CSPNe. Finally, distances inferred from extinction jumps were adopted as PN distances only after satisfying multiple independent consistency checks.

3.1 Extinction–distance profile construction

The first and most critical step is to determine the extinction for each star within R_field, expressed as $$E_{{\mathrm{G}}_{\mathrm{B}\mathrm{P}},{\mathrm{G}}_{\mathrm{R}\mathrm{P}}}=(G_{\mathrm{B}\mathrm{P}}-G_{\mathrm{R}\mathrm{P}})-{(G_{\mathrm{B}\mathrm{P}}-G_{\mathrm{R}\mathrm{P}})}_0,$$(1)

where the observed colour (G_BP − G_RP) is taken from Gaia DR3 photometry, and the intrinsic colour (G_BP − G_RP)₀ is obtained via the widely used blue-edge method (Wang & Jiang 2014; Wang & Chen 2019, 2024; Deng et al. 2025), which estimates intrinsic colours from stellar atmospheric parameters (T_eff, log g, and [M/H]). We adopted the machine-learning-enhanced blue-edge method of Zhao et al. (2024), with a typical intrinsic-colour uncertainty of 0.014 mag. The derived extinction values do not rely on a fixed R_V and offer improved accuracy over previous methods. Stellar distances were taken from BJ21. These extinction and distance estimates form the basis of the extinction–distance profiles used in the subsequent analysis.

3.2 Extinction jump identification

We identified significant extinction jumps in the extinction–distance profiles using an extinction-jump model. These jumps appear as sharp increases in extinction over narrow distance ranges and define extinction-jump distances. Planetary nebulae can produce such detectable extinction features through dust within the nebula itself and in its surrounding extended AGB dust, whose spatial extent far exceeds that of the ionised nebula.

Our model, following Chen et al. (2017), Wang et al. (2020), Zhao et al. (2020), Sun et al. (2021), and Chen et al. (2025), assumes that cumulative extinction along the line of sight arises from two components: diffuse interstellar dust and PN-associated dust. The total extinction at distance d is then expressed as $$A(d)=A_0(d)+A_1(d),$$(2)

where A₀(d) represents the foreground diffuse extinction and A₁(d) accounts for the PN-associated extinction. For PNe, we approximated the foreground extinction as a constant, A₀(d) = A₀, following the treatment adopted for nearby molecular clouds by Cao et al. (2023). This approximation is motivated by the weak foreground extinction towards nearby PNe and the limited number of foreground stars, which makes a distance-dependent foreground difficult to constrain. The parameter A₀ was fitted as a free parameter and represents the foreground extinction level prior to the jump. The PN-associated extinction is modelled as $$A_1(d)=\mathrm{\Delta }A\frac{1+\mathrm{erf}\left[\frac{d-d_0}{\delta d}\right]}{2},$$(3)

where ΔA is the amplitude of the extinction jump and d₀ the distance to the PN centre. Here, δd is set to the PN diameter as a smoothing scale, but in practice it mainly reflects the propagation of stellar distance uncertainties (see Section 5.1). The fitting parameters are A₀, ΔA, and d₀. We employed Markov chain Monte Carlo (MCMC) to fit the extinction–distance profiles, deriving best-fit parameters (A₀, ΔA, d₀). Parallax uncertainties were propagated into the distance errors and incorporated into the MCMC likelihood. The uncertainty associated with each candidate distance solution was quantified by the 1σ statistical uncertainty of the fitted jump location d₀, as inferred from the MCMC posterior distribution. An extinction jump was considered significant when ΔA > 0.07 mag (five times the intrinsic colour uncertainty) and the distance parameter d₀ was well constrained with a single dominant peak in the posterior distribution. For cases with multiple peaks in the MCMC posterior indicating several extinction jumps, multiple components were fitted simultaneously.

3.3 PN distance estimation

In the previous step, each significant extinction jump identified in the extinction–distance profile yields an extinction–jump distance. Not all extinction jumps detected along a given line of sight are associated with the PN, as similar extinction features may also arise from interstellar clouds. Thus, these distances were evaluated using available comparisons with 3D dust extinction maps, analysis of the stellar spatial distribution relative to the PN centre, and the use of published distance constraints with different priority levels. The final adopted PN distance was selected based on the combined results of these checks.

All profiles were compared with two representative 3D interstellar dust extinction maps. The first is the recent map from Zhang & Green (2025, hereafter ZG25) which provides extinction curves, individual stellar extinction values, and distances for 130 million stars. The second is the widely used map from Green et al. (2019, hereafter G19), covering three-quarters of the sky and providing extinction–distance curves at 3.4′–13.7′ resolution. Within each PN’s R_field, ZG25 single-star data were compared point-by-point with our measurements, while G19 curves were used to verify overall trends and jump locations at larger scales.

For profiles exhibiting a single extinction jump, the distance was considered reliable when the jump was further supported either by agreement with the CSPN distance or by a strong extinction increase traced by stars nearest the PN centre. As illustrated in Figure 2, when a CSPN distance is available (top left panel), our extinction–distance measurements (coloured dots and black line) are consistent with the discrete ZG25 data (grey dots) and the G19 curve (green line) and are further supported by the CSPN distance. When no CSPN distance is available (top right panel), the extinction jump is instead supported by a clear extinction increase traced by stars closest to the PN centre (red dots, ΔE_GBP,GRP = 1.0 mag), providing an independent indication of the PN location.

When multiple extinction jumps produce more than one extinction-jump distance, these distances were further evaluated using distance constraints with different priority levels. These constraints are based on multiple independent distance estimates, together with their associated uncertainties. The highest-priority indicators are CSPN-based distances (CW21, GS21). When available, no other constraints were considered. The next priority is parallax-based CSPN distances (BJ21). The remaining indicators include canonical statistical distances (F16, with mean values adopted following CW21), CSPN-calibrated statistical distances (A22, BS23), and the latest distances (HJ24), all considered of equal weight. When indicators carry equal weight, distances supported by more independent constraints are considered more reliable. The bottom panels of Figure 2 present two examples. The lower-left panel shows a PN exhibiting two extinction jumps along the line of sight, at ~1 kpc (ΔE_GBP,GRP = 0.4 mag) and ~2.5 kpc (ΔE_GBP,GRP = 0.2 mag) respectively. The first jump is identified as the PN location based on the highest-priority indicators (CW21, GS21). The lower-right panel shows the high-confidence region where the HJ24 and F16 constraints overlap, with the core-source jump clearly detected (ΔE_GBP,GRP = 1.0 mag), confirming the PN location.

Throughout this paper, our results are shown in black and literature results in colour, with stellar distances as dots and cloud distances as diamonds. The colour coding assigns yellow to HJ24, cyan to BS23, purple to A22, orange to CW21, light green to GS21, blue to F16, and grey to BJ21. In rare cases where the distance constraints are insufficient to estimate the PN distance, the R_PN–distance relation offers auxiliary validation (Section 5.3).

4 Method validation

To evaluate the reliability of our method, we selected a tight sample of 137 PNe with both CW21 and GS21 distances consistent within 3σ. Distances derived for this sample using our method were then compared with the reference values, yielding residuals with a mean of 0.00 and dispersions of approximately 0.08 kpc for GS21 and 0.11 kpc for CW21, confirming the robustness of our distance estimates (Figure 3).

In addition to tight-sample validation using CSPNe distances, we also include other representative individual PN distance estimates, such as kinematic, outburst-brightness, pulsation-based, and convergent parallax distances. Kinematic distances, derived from sky position and radial velocity under an assumed Galactic rotation model, are particularly suited for young, disk PNe (typically Type I). As in F16, we adopted the best available kinematic distances (e.g. Jones et al. 2010; Sabin et al. 2012). Several historically observed final-flash CSPNe were measured via the outburst-brightness method (F16). Additional representative distances include FG Sge from pulsation theory (Mayor & Acker 1980), Sh 2-188 from PN–ISM interaction modelling combined with CSPN proper motion (Wareing et al. 2006), and NGC 7293 from convergent parallax assuming Hyades membership (Eggen 1984). We also included three PNe hosting rare eclipsing-binary CSPNe, whose distances were derived from the eclipsing-binary method (Bell et al. 1994; Exter et al. 2005, 2010). Table 2 summarises these distances with their methods and references. The strong agreement among these diverse techniques in Figure 4 further supports the reliability and robustness of our method.

Fig. 2 Extinction–distance profiles for four PNe with single jumps in the top row and multiple jumps in the bottom row. PN names are labelled at the lower right of each panel. Distance (in kiloparsec) vs extinction EGBP,GRP (in magnitude). The source colours are the same as those in Figure 1, and the solid black lines show the fits obtained with our extinction-jump model. PN positions are shown with vertical dashed black lines. The grey points represent ZG25 stars, and the dashed green lines indicate G19 curves. PN distance constraints are shown as diamonds with error bars from GS21 (light green), CW21 (orange), HJ24 (yellow), and F16 (blue). Conventions are consistent throughout.

Fig. 3 Comparison of distances for the tight sample. Distances derived in this work are compared with those from CW21 (orange circles) and GS21 (blue circles). All data points are shown with their associated uncertainties. The solid black line shows y = x, while the coloured lines represent the best-fitting relations. The inset displays the distribution of residuals.

5 Results and discussion

5.1 PNe distances

We applied our method to 2225 PNe (Section 2.1), obtaining extinction distances for 1066 objects, 87% of which have relative distance uncertainties below 20%, with a median value of 13%, as listed in Table 3. To quantify the degree of independent external support, we assigned a quality flag Q (1–3) to each adopted distance: Q = 3 denotes cross-validation by CSPN-based distance estimates together with at least one additional independent distance indicator; Q = 2 denotes support from CSPN-based distance estimates alone or, in the absence of CSPN information, by at least two mutually consistent independent distance indicators; Q = 1 denotes support from only one independent indicator. The Q flag characterizes the strength of cross-validation rather than the intrinsic robustness of the extinction-based distances themselves. Among these, 765 PNe have CSPN parallaxes with relative uncertainties ≥20%, 128 lack CSPN parallaxes, and 124 distances were determined for the first time. Figure 5 compares our distance distribution with CSPN-based distances (CW21 + GS21, 897 PNe), as well as with the largest-to-date and the most recent extinction–distance samples (G11, 70 PNe and D21, 17 PNe). Although our observational depth remains within the limits of previous work, our results substantially increase the number of distances to PNe derived using individual methods.

Extinction jumps detected by our method span 0.15 mag < ΔE_GBP,GRP < 2.8 mag, consistent with the largest recent PN extinction sample from D21 (ΔE_GBP,GRP ≈ 0.04–3.6 mag using the A_V–ΔE_GBP,GRP relation of Wang & Chen 2019). Planetary nebulae with higher extinction (ΔE_GBP,GRP > 2.8 mag) were not detected by our method, likely due to the paucity of visible stars in such dense regions, marking the method’s detection limits.

Planetary nebulae typically have physical radii below 1.5 pc, reaching 2.5–3 pc at late evolutionary stages (Frew & Parker 2010; Parker 2022), and our results are consistent with this, with a median radius of 0.45 pc and 90% of the sample spanning 0.12–1.51 pc. The apparent width of extinction jumps in distance space, however, is generally larger than the physical scale of PNe, with a median of about 150 pc, indicating that it reflects the effective distance resolution of the extinction–distance method rather than any intrinsic physical thickness. We performed a simple simulation by introducing a 10% relative uncertainty in tracer-star extinctions and distance uncertainties of 0%, 10%, and 20% for extinction jumps d at 1–3 kpc. The inferred jump width, Δd, increases systematically with distance uncertainty, growing from ~6 pc at d = 1 kpc with no distance uncertainty to ~110–150 pc for 10–20% uncertainties and reaching ~0.5 kpc at d = 3 kpc for a 20% uncertainty. This trend confirms that the jump width primarily reflects the propagation of distance uncertainties and is consistent with our earlier conclusions.

Our distances are compared with literature values in Figure 6. Orthogonal fits were performed with a single 3σ rejection and refitting, and the discarded sources are shown in grey. They agree best with CSPN-based results from GS21 and CW21. A22 and BS23 tend to underestimate distances, whereas HJ24 exhibits a pronounced systematic error growth with distance. F16 matches well. The grey outliers in GS21 and CW21 mainly reflect inconsistent CSPN identifications (Section 5.3.2), while other outliers trace method-related uncertainties (Section 5.3.1). In general, the final fits are close to one-to-one (μ = 0), with rejected outliers accounting for only 1.5 of our sample. Our method yields smaller uncertainties while remaining consistent with the literature, thereby providing more robust distance estimates.

To further assess the reliability of our distance estimates for distant PNe, particularly along low-latitude sightlines where extinction is strong, we additionally compared our results with the 3D extinction map of Marshall et al. (2006). Constructed from 2MASS near-infrared photometry combined with a Galactic stellar population model, this map remains effective in high-extinction regions and therefore provides a targeted complementary test for distant PNe. For a meaningful comparison, a PN must lie within both the sky coverage of the map, and its sampled distance must range along the relevant line of sight. Within this common range, 37 distant PNe with d > 4.5 kpc can be directly compared. Among these, ten coincide with or lie close to clear extinction jumps in the Marshall profiles, while the others fall in flat segments of the curves without a distinct jump. Within this subset, the jump-associated objects tend to have larger angular radii (mostly >28″), whereas the remaining ones are smaller. Given the angular resolution of the map (15′ × 15′) and its discrete distance sampling, small-scale extinction structures are expected to be smoothed by spatial averaging. No systematic offset relative to our distance scale is found for this subset.

Fig. 4 Comparison of our distances with selected individual PN measurements from Table 2. The yellow triangles indicate kinematic distances. The blue circles denote other independent methods.

Fig. 5 Step histograms of PN distance distributions for this work (black), CSPN-based samples (blue), and extinction–distance samples (red).

Fig. 6 Comparison of our distances (x-axis, kiloparsec) with literature values (y-axis, kiloparsec). Each panel shows a reference scale (CW21, GS21, F16, A22, BS23, and HJ24). The solid black line shows y = x, and the coloured lines represent orthogonal linear fits obtained with iterative 3σ rejection for each dataset. The grey points indicate the objects excluded (>3σ) from the final fit. The insets show the residual distributions with the mean (μ) and scatter (σ).

Fig. 7 Spatial distribution of the PN sample in Galactic coordinates (l, b). The colour indicates distance, and the point size denotes RPN (10″, 100″, 500″).

5.2 PNe spatial distribution

Contrary to the long-held view that the extinction–distance method could apply only within a few degrees of the Galactic plane (Phillips 2006; Frew et al. 2016; Kwitter & Henry 2022), we find that its limitation stems from data coverage. Previous PN extinction studies relied mainly on Hα surveys (e.g. G11; D21), which cover only the Galactic plane (e.g. IPHAS, |b|<5°; Drew et al. 2005). As shown in Figure 7, Gaia’s all-sky coverage and improved method allow us to probe higher latitudes and even halo PNe. About 85.5 of PNe lie within |b| < 10°, concentrated towards the plane. At higher latitudes, detections drop sharply, leaving larger, nearby PNe mostly within ≲1 kpc.

5.3 PN radius–distance relation

A clear trend between R_PN and distance is observed in our sample, with a power-law exponent of γ = −1.62. Applying the same procedure to the full datasets of HJ24, BS23, A22, CW21, GS21, and F16 produces exponents from −2.3 to −1.3 (mean residual 0.17), with our value at the median. Figure 8 displays all R_PN–distance relations. Our fit runs nearly parallel to the others, with all results lying within or near our 1σ range, confirming consistency with the established relation. In this work, the R_PN–distance relation was not used as a distance scale, but only as a supplementary check in a very few exceptional cases.

Fig. 8 Comparison of fitted RPN–distance relations. The grey points indicate our PNe. The black line and shaded area show our best fit and the 1σ range. The dashed coloured lines represent literature fits obtained from their original complete datasets, using the same colour coding as in Figure 2.

5.3.1 Applications: PN distance diagnostics

As noted in Section 3.3, when distance constraints are insufficient to estimate the PN distance, or when derived distances differ markedly from literature values (e.g. the grey points in Figure 6), the R_PN–distance relation serves as a supplementary diagnostic. The applications discussed below are illustrated by three representative cases (panels a–c of Figure 9), with the corresponding R_PN–distance residuals shown in panel d.

Figure 9a shows two extinction jumps within the constraint range of HJ24, illustrating a complex sightline in which conventional methods fail to uniquely locate the PN. The R_PN–distance relation favours the first jump as the plausible solution (residual 1.7 vs 6.1). Figure 9b presents a PN exhibiting a large methodological discrepancy, where our estimate differs from HJ24 by a factor of 5.3. Observational limits prevent detection at the HJ24 position. The R_PN check (panel d) favours our distance (residual 2.67) over HJ24 (17.27; ≫ 3σ from the relation, with a 68% distance uncertainty), demonstrating robustness for large-distance deviations.

Figure 9c shows PN Fr 2-36 whose distance differs significantly from GS21 and HJ24 (by factors of 10 and 13, respectively), likely due to a misidentified CSPN. First, the PN lies at high Galactic latitudes (~23°) along a clean, low-extinction sightline free from background contamination. The single, sharp extinction jump at ~0.49 kpc (maximum ΔE_GBP,GRP ≈ 1 mag), consistent with the PN core and the trend in ZG25, thus reliably traces the PN. Second, only one of the four required priors (GS21) applies to the HJ24 distance estimate for this PN, yielding a 74% uncertainty. The resulting overlap with GS21 therefore reflects a shared prior rather than independent evidence for distance reliability. Finally, our distance minimises the R_PN–distance residual (0.23 kpc), while the literature values deviate by ≫ 3σ (6.57 and 4.32 kpc). With R_PN = 1500″ adopted from HASH, as also used in CW21), a ~5 kpc distance would imply an unphysically large PN. All indications strongly suggest that this CSPN has been misidentified.

5.3.2 Applications: Discordant CSPN identifications

When CW21 and GS21 assigned different CSPNe to the same PN, we selected the most plausible one based on (1) stellar parameters (T_eff, log g) consistent with white dwarfs, (2) consistency among parallax distance, extinction jumps, and published constraints, and (3) proximity to the R_PN–distance relation. A CSPN was regarded as well-supported if it met at least two of these criteria. Among 33 such PNe, we identified the more reliable CSPN in 15 cases (Table 4).

As an illustration, in the top panel of Figure 10, the two CSPNe of NGC 2899 can be distinguished using the latter two criteria. An analysis of extinction jumps shows three jumps along the line of sight. One CSPN lies at 2.0 kpc with a 5.9 uncertainty, precisely matching the second jump. This jump encloses the core star, rises sharply by 0.4 mag, and is consistently covered by three independent distance measurements, showing good agreement. The other CSPN is at 4.4 kpc with a 69% uncertainty, with its error range covering only the third jump. The 2.0 kpc CSPN also has a smaller residual with respect to the R_PN–distance relation, indicating that it is the more reliable CSPN. The bottom panel shows the R_PN–distance relations for all PNe listed in Table 4, shown as triangles, with the PN discussed here indicated by circles. Each PN is shown in the same colour for its associated CSPNe. In all cases, the adopted CSPN (solid dot) lies closer to the fitted relation than the other (open circle), confirming the internal consistency of our identifications and the robustness of the R_PN–distance relation.

Collectively, these tests confirm the reliability and internal consistency of the method. They demonstrate that the extinction–distance approach provides a robust basis for constructing a large and homogeneous PN distance catalogue.

6 Summary and conclusions

We presented an optimized extinction-distance method based on Gaia data, providing a CSPNe parallax-independent and physically motivated framework for determining distances to PNe with confirmed CSPNe. Our main results are as follows:

1. We provided extinction-based distance estimates for 1066 PNe, including 124 objects measured for the first time. The relative distance uncertainty has a median value of 13% and is below 20% for 87% of the sample. Among them, 765 CSPNe have relative parallax errors exceeding 20%, and 128 lack parallax measurements, demonstrating that our method provides a robust complement when parallaxes are uncertain or unavailable;

2. We significantly extended the reliability and applicability of traditional extinction–distance measurements, enabling distance determinations for high-Galactic-latitude PNe. The sample of PNe with individually measured distances increases from a few dozen to over one thousand, forming the largest catalogue of extinction-based PN distances to date;

3. We examined 33 PNe with multiple and inconsistent CSPNe identifications, recommend more reliable CSPNe for 15 objects, and find that the uniquely reported CSPN of Fr2-36 is likely misidentified.

Future work will further refine extinction estimates, expand the sample, and incorporate forthcoming Gaia DR4 data to investigate the extinction properties and local environments of PNe in greater detail.

Fig. 9 Panels (a–c): extinction–distance profiles for three PNe, with names labelled at the lower right. Sources and lines are colour-coded, as in Figure 2. Panel (d): residuals ΔD = D − DRfit, where DRfit is the distance predicted by the RPN–distance fit for each PN. PN names are listed along the x-axis, with the adopted literature value of RPN indicated below each PN. The data points adopt the colours from panels (a) to (c), with our determinations highlighted as black diamonds. All residuals are labelled above their respective points.

Fig. 10 CSPNe validation. Top: extinction–distance profile of NGC 2899. The sources, reference distances, and lines follow the colour coding in Figure 2. The two published CSPNe are shown in grey with error bars. Bottom: discordant CSPNe. The black line represents the fitted RPN–distance relation with the shaded area marking the ±1σ range. All PNe from Table 4 are included. The filled symbols indicate the adopted CSPN, and the open symbols indicate other published CSPNe. NGC 2899 is shown in grey circles, whereas CSPNe in other PNe are shown as coloured triangles.

Data availability

Full Table 3 is available at the CDS via the VizieR online service and on Zenodo at https://doi.org/10.5281/zenodo.17010405.

Acknowledgements

We sincerely thank the anonymous referee for the careful reading of the manuscript and for the very constructive and helpful comments, which we greatly appreciate. We thank Prof. Xiaodian Chen for very helpful discussions and suggestions. This work is supported by the National Natural Science Foundation of China (NSFC) projects 12133002 and 12373028. S.W. acknowledges the support from the Youth Innovation Promotion Association of the CAS (grant No. 2023065). This work has made use of data from Gaia survey.

References

All Tables

All Figures

	Fig. 1 Left: stars overlaid on the Digitized Sky Survey image with RPN (dashed yellow circle), the CSPN (yellow star), as well as core, nebular, and peripheral sources (red, dark blue, and blue dots). Right: extinction vs distance.
In the text

Fig. 2 Extinction–distance profiles for four PNe with single jumps in the top row and multiple jumps in the bottom row. PN names are labelled at the lower right of each panel. Distance (in kiloparsec) vs extinction EGBP,GRP (in magnitude). The source colours are the same as those in Figure 1, and the solid black lines show the fits obtained with our extinction-jump model. PN positions are shown with vertical dashed black lines. The grey points represent ZG25 stars, and the dashed green lines indicate G19 curves. PN distance constraints are shown as diamonds with error bars from GS21 (light green), CW21 (orange), HJ24 (yellow), and F16 (blue). Conventions are consistent throughout.

In the text

	Fig. 3 Comparison of distances for the tight sample. Distances derived in this work are compared with those from CW21 (orange circles) and GS21 (blue circles). All data points are shown with their associated uncertainties. The solid black line shows y = x, while the coloured lines represent the best-fitting relations. The inset displays the distribution of residuals.
In the text

	Fig. 4 Comparison of our distances with selected individual PN measurements from Table 2. The yellow triangles indicate kinematic distances. The blue circles denote other independent methods.
In the text

	Fig. 5 Step histograms of PN distance distributions for this work (black), CSPN-based samples (blue), and extinction–distance samples (red).
In the text

Fig. 6 Comparison of our distances (x-axis, kiloparsec) with literature values (y-axis, kiloparsec). Each panel shows a reference scale (CW21, GS21, F16, A22, BS23, and HJ24). The solid black line shows y = x, and the coloured lines represent orthogonal linear fits obtained with iterative 3σ rejection for each dataset. The grey points indicate the objects excluded (>3σ) from the final fit. The insets show the residual distributions with the mean (μ) and scatter (σ).

In the text

	Fig. 7 Spatial distribution of the PN sample in Galactic coordinates (l, b). The colour indicates distance, and the point size denotes RPN (10″, 100″, 500″).
In the text

	Fig. 8 Comparison of fitted RPN–distance relations. The grey points indicate our PNe. The black line and shaded area show our best fit and the 1σ range. The dashed coloured lines represent literature fits obtained from their original complete datasets, using the same colour coding as in Figure 2.
In the text

Fig. 9 Panels (a–c): extinction–distance profiles for three PNe, with names labelled at the lower right. Sources and lines are colour-coded, as in Figure 2. Panel (d): residuals ΔD = D − DRfit, where DRfit is the distance predicted by the RPN–distance fit for each PN. PN names are listed along the x-axis, with the adopted literature value of RPN indicated below each PN. The data points adopt the colours from panels (a) to (c), with our determinations highlighted as black diamonds. All residuals are labelled above their respective points.

In the text

Fig. 10 CSPNe validation. Top: extinction–distance profile of NGC 2899. The sources, reference distances, and lines follow the colour coding in Figure 2. The two published CSPNe are shown in grey with error bars. Bottom: discordant CSPNe. The black line represents the fitted RPN–distance relation with the shaded area marking the ±1σ range. All PNe from Table 4 are included. The filled symbols indicate the adopted CSPN, and the open symbols indicate other published CSPNe. NGC 2899 is shown in grey circles, whereas CSPNe in other PNe are shown as coloured triangles.

In the text