Resolving diffusion signatures in distant pulsar halos with current and future experiments

In 2007, the water Cherenkov experiment Milagro discovered γ-ray emission extending $ \approx3^\circ $ around the bright γ-ray pulsar Geminga [1]. This large extended halo was difficult to attribute to the pulsar itself or its associated bow-shock pulsar wind nebula [2, 3]. A decade later, the High-Altitude Water Cherenkov Observatory (HAWC) observed this halo with improved angular resolution and sensitivity, revealing that an electron diffusion model could well explain its spatial distribution [4]. This finding pointed to the most likely origin of the Geminga halo: electron-positron pairs ($ e^\pm $) produced by the pulsar escape into the interstellar medium (ISM) and, as they diffuse through the turbulent magnetic field, produce extended γ-ray emission through inverse Compton scattering with background photons. Consequently, a new class of astrophysical sources emerged, known as pulsar halos or TeV halos. The morphology of a pulsar halo provides a clear projection of the spatial distribution of its parent $ e^\pm $, making it a unique probe for studying cosmic-ray propagation in a specific ISM region [5−8].

Based on surveys of TeV γ-ray sources conducted by the wide-field experiments HAWC and the Large High Altitude Air Shower Observatory (LHAASO), the number of confirmed or candidate pulsar halos has reached approximately ten [9, 10]. However, apart from the Geminga halo and the Monogem halo [4], the pulsars associated with other potential halos are all located at distances exceeding $ 1\ \text{kpc} $, resulting in significantly smaller angular extensions of these halos. In the case of LHAASO J0621$ + $3755, morphological measurements limited by instrumental angular resolution cannot distinguish between a diffusion model and mathematical templates such as disk or Gaussian profiles [11]. The lack of clear diffusion signatures in the observed morphology complicates the confirmation of these sources as pulsar halos.

Enhanced instrumental angular resolution plays a crucial role in reliably expanding the sample size of pulsar halos. The Cherenkov Telescope Array (CTA, [12]), currently under construction, will achieve unprecedented angular resolution ($ <0.05^\circ $ above $ 10\ \text{TeV} $ [13]) in γ-ray astronomy, along with significantly improved sensitivity compared to current atmospheric Cherenkov experiments. This capability will serve as the primary means of identifying diffusion signatures in distant pulsar halos. Meanwhile, LHAASO, with its wide field of view and long exposure time, remains indispensable for precise observations of relatively nearby pulsar halos. In the future, LHAASO plans to integrate deep learning techniques to improve directional reconstruction, thereby achieving the enhanced angular resolution necessary to resolve the morphologies of more distant pulsar halos. In light of these prospects, this study examines the potential of experiments to decipher diffusion signatures in candidate pulsar halos and evaluates how enhancements in angular resolution or exposure time can improve their identification.

The paper is organized as follows. Section II presents a simplified template of pulsar halos that effectively captures their morphological characteristics. Section III introduces the γ-ray experiments and the methods for simulating pulsar halo signals based on experimental parameters. Section IV presents the statistical approach used for morphological discrimination. Section 5=V demonstrates the capability of these experiments to conduct morphological tests on pulsar halos. Section VI provides the conclusion.

Electron-positron pairs produced by pulsars can be accelerated to energies as high as $ 100\ \text{TeV} $ by relativistic shocks in PWNe, from which they continuously escape into the ISM. When a PWN evolves into the bow-shock stage, its size remains on the order of $ 1\ \text{pc} $ for an extended period, confined by the ISM ram pressure [14]; this is much smaller than the characteristic diffusion scale of the escaped $ e^\pm $. Consequently, the propagation of the parent $ e^\pm $ in pulsar halos can be treated as a diffusion–radiation-loss problem under the point-source assumption. In this scenario, the γ-ray surface brightness of a pulsar halo can be approximately described by the following parameterization [4]:

where $ f_{0} $ is the normalization at $ 30\; \text{TeV} $, α is the spectral index, $ E_c $ is the cutoff energy, θ is the angular distance from the pulsar, and $ \theta_d $ is the characteristic angular extension of the halo, which can depend on $ E_\gamma $. The γ-ray spectrum is assumed to be a power-law with an exponential cutoff (ECPL) form, reflecting the fact that the parent $ e^\pm $ originate from a shock-acceleration process.

The spatial distribution term in Eq. (1) exhibits a pronounced centrally peaked profile, which differs significantly from the uniform disk template commonly used to describe the profiles of supernova remnants. Additionally, a Gaussian template is generally used in γ-ray astronomy to characterize the morphology of extended sources with a centrally peaked profile. However, it also differs remarkably from Eq. (1): within the range $ 0.06\theta_d \lt \theta \lt \theta_d $, Eq. (1) decreases sharply with increasing θ, exhibiting a steeper decline than the Gaussian distribution. Nevertheless, all these distinctions may become blurred due to insufficient instrumental angular resolution or limited photon statistics.

In this work, we aim to investigate the capability of γ-ray experiments to identify pulsar halos through morphological discrimination. We adopt the diffusion-based parameterization (hereafter referred to as the diffusion model) presented in Eq. (1) to describe the morphology of pulsar halos, and consider the disk and Gaussian models as representative alternatives for non-halo sources.

The variables that primarily reflect differences among pulsar halos are the pulsar distance d and its spin-down luminosity $ \dot{E} $. Generally, morphological discrimination becomes more difficult for more distant sources or those with lower spin-down luminosity. At present, the γ-ray energy spectrum and morphology of the Geminga halo are the most precisely measured. For simplicity, we assume that all halos share the same spectral shape as the Geminga halo. The γ-ray energy spectrum of the Geminga halo is obtained by fitting the HAWC data [15], yielding $ \alpha=1.49 $ and $ E_c=22.7 $ TeV, which serves as the benchmark model to generate mock data for all halo candidates.

Additionally, we assume that $ f_0 $ and $ \theta_d $ for a pulsar halo scale with the ratio of the pulsar's distance to that of Geminga: $ f_0=f_{0,G}(d/d_G)^{-2} $ and $ \theta_d=\theta_{d,G}(d/d_G)^{-1} $, where $ \theta_{d,G} $ is set to $ 5.5^\circ $ as measured by HAWC [4]. This scaling implies the underlying assumption that all pulsar halo candidates share the same pair-conversion efficiency from spin-down power as the Geminga pulsar halo. For simplicity, $ \theta_d $ is assumed to be energy independent.

We acknowledge that adopting a unified spectral shape and an energy-independent extension for all pulsar halos is a simplification. However, known pulsar halos exhibit similar spectral shapes, and existing energy-dependent morphological measurements indicate that the angular extension of pulsar halos in the TeV regime does not vary significantly [15, 16]. This suggests that our assumption is consistent with current observations.

Ground-based observations of high-energy γ rays are conducted indirectly by measuring secondary particles generated in extensive air showers initiated by primary γ rays. Two main experimental techniques are employed to detect this emission. The first is the Imaging Atmospheric Cherenkov Telescope (IACT) technique, which observes Cherenkov light flashes generated by relativistic shower particles propagating through the atmosphere during clear, moonless nights. The second is the extensive air shower (EAS) technique, which detects secondary shower particles directly at ground level. In this work, we focus on two representative instruments employing these approaches: CTA and the Kilometer Square Array of LHAASO (LHAASO-KM2A).

CTA will be the flagship next-generation instrument in γ-ray astronomy, covering photon energies from 20 GeV to 300 TeV. It is designed to achieve an energy resolution better than $ 10\% $ and a substantially improved angular resolution compared to existing γ-ray telescopes. CTA will operate as a global observatory with two sites: a Northern Hemisphere array (CTA-North) located on La Palma, Spain, and a Southern Hemisphere array (CTA-South) located at Paranal, Chile. In the Alpha Configuration, CTA-North will consist of four Large-Sized Telescopes (LSTs) and nine Medium-Sized Telescopes (MSTs), while CTA-South will comprise 14 MSTs and 37 Small-Sized Telescopes (SSTs). The primary mirror diameters of the LSTs, MSTs, and SSTs are 23 m, 11.5 m, and 4.3 m, respectively. Their corresponding fields of view are greater than $ 4.5^{\circ} $, $ 7^{\circ} $, and $ 8^{\circ} $ [17].

Instrument Response Functions (IRFs)—including angular and energy resolutions, effective area, and expected background rates—are commonly used to derive differential sensitivities and to enable comparisons between different γ-ray instruments. To facilitate a direct comparison with LHAASO-KM2A under similar sky coverage, we adopt the CTA-North IRFs for on-axis observations. These IRFs are publicly available through the Zenodo repository [18]. Following standard event selection procedures, including γ–hadron separation cuts, we obtain the photon effective area and residual background rates, consistent with the methodology described in Ref. [19].

LHAASO is a new-generation γ-ray and cosmic-ray observatory located at an altitude of 4410 m in Daocheng, Sichuan Province, China ($ 29^\circ 21' 31'' $ N, $ 100^\circ 08' 15'' $ E) [20]. It consists of three sub-arrays: the LHAASO-KM2A, the Water Cherenkov Detector Array (LHAASO-WCDA), and the Wide Field-of-View Cherenkov Telescope Array (LHAASO-WFCTA). Partial operation of LHAASO began in 2019, and full construction was completed in 2021. LHAASO-KM2A is composed of 5216 electromagnetic detectors (EDs) and 1188 muon detectors (MDs), covering an area of 1.3 $ \rm{km}^2 $. The EDs (MDs) are deployed with a spacing of 15 m (30 m). The EDs are plastic scintillation detectors that measure the electromagnetic component of EASs and are designed primarily for event reconstruction. The MDs detect the muonic component of EASs and are mainly used for cosmic-ray background discrimination. High-energy γ-ray events above 10 TeV are predominantly recorded by KM2A. The array has a wide field of view of approximately 2 sr and operates with a nearly full duty cycle, providing unprecedented sensitivity for surveying the γ-ray sky at energies above 20 TeV.

We adopt the LHAASO-KM2A instrument response functions (IRFs) calibrated using observations of the Crab Nebula [21]. Specifically, we use the IRFs corresponding to a zenith angle of $ \theta = 20^\circ $, ensuring consistency with those released by the other experiments considered in this work. The expected background rate is derived by convolving the all-particle cosmic-ray flux [22] with the effective area. The cosmic-ray (CR) effective area is derived from the reference γ-ray effective area by accounting for their respective selection efficiencies. Assuming the initial effective areas for CRs and γ-rays are identical prior to selection, the CR effective area is calculated as $ A_{\text{eff}}^{\text{CR}} = (A_{\text{eff}}^{\gamma} / \epsilon_{\gamma}) \cdot \epsilon_{\text{CR}} $, where $ \epsilon_{\gamma} $ and $ \epsilon_{\text{CR}} $ denote the selection efficiencies for γ-rays and CRs after cuts, respectively [21]. This procedure ensures that the angular resolution, effective area, and expected background rate remain consistent with the benchmarks established in Ref. [19].

A direct comparison between the performances of LHAASO-KM2A and CTA is now possible, as the same energy binning—namely 0.2 per decade on a logarithmic scale—is adopted for both instruments. Based on their respective performance characteristics and the energy ranges relevant for halo spectral measurements, we consider energies from 10 TeV to 100 TeV for LHAASO-KM2A and from 1 TeV to 100 TeV for CTA.

Owing to its pointing observation mode and imaging camera, CTA achieves excellent angular resolution. At 30 TeV, the $ 1\sigma $ width of the point spread function (PSF) is $ \sigma_{\rm PSF} \sim 0.03^\circ $ for CTA, compared to $ \sigma_{\rm PSF} \sim 0.3^\circ $ for LHAASO-KM2A. However, because CTA detects Cherenkov light produced by air showers, it can operate only under clear, dark-sky conditions. In contrast, LHAASO-KM2A has a nearly full duty cycle and can operate continuously, both day and night.

For a given extended halo, the differential flux at any spatial position and energy is determined by the assumed spectral and morphological models. The expected number of signal events, $ N_s $, per energy bin and per spatial pixel within the region of interest (ROI) is calculated by folding the model flux with the instrument response functions, including the PSF, and multiplying by the effective area, observation time, and solid angle. In this work, the model flux is defined by the diffusion morphology and the ECPL spectrum introduced in Section 2.1. The number of background events, $ N_b $, is computed from the background rate, which is assumed to be constant over time. Using a toy Monte Carlo (toyMC) approach, the mock observed event counts, $ N_{\rm on} $, are generated by drawing Poisson realizations of $ N_s + N_b $.

For a given d and $ \dot{E} $, mock data can be generated based on the diffusion model using the toy MC technique outlined in Section III C. Next, the simulated data are fitted with the model to be tested using the maximum-likelihood method to determine the best-fit parameters. Subsequently, we evaluate the best-fit model using a specific statistical test to reach a decision for this single trial.

We repeat this procedure multiple times and compile the decisions from each trial. The proportion of trials in which the diffusion model is correctly chosen (or the incorrect model is rejected) relative to the total number of trials is defined as the power of the test. The power reflects the capability of the considered experiment to distinguish the diffusion model.

In the ROI, a joint maximum-likelihood fit is performed over energy and spatial bins. Assuming Poisson statistics, the likelihood function is written as

where P denotes the Poisson probability distribution, i labels the i-th energy bin, and j labels the j-th spatial pixel. Here, $ N_{ij}^{\rm on} $ is the number of observed events, while $ N_{ij}^{\rm b} $ and $ N_{ij}^{s} $ are the expected background and signal counts, respectively.

The signal expectation $ N_{ij}^{s} $ is determined by the spatial model under test (Gaussian, disk, or diffusion) and the ECPL spectral model outlined in Section 2, folded with the instrument response functions. To simulate a realistic data analysis pipeline, we treat $ f_0 $, α, $ E_c $, and the angular width of the spatial model as unconstrained free parameters during the fitting process, identical to the procedure applied to experimental data. The likelihood is maximized using the $\mathrm{MIGRAD}$ algorithm implemented in the Minuit Minimization package ¹. Higher $ \ln \mathcal{L} $ values indicate a stronger preference for the presence of an extended source with the assumed morphology.

We adopt two distinct statistical metrics at the single-trial level to discriminate among different spatial models: the Akaike Information Criterion (AIC) and the $ \chi^2 $ goodness-of-fit test. When data are sufficient, discriminating in favor of the diffusion model over an incorrect model is generally easier than outright rejecting that incorrect model. Therefore, the AIC typically yields higher statistical power than the $ \chi^2 $ test. However, compared to the AIC, the decision criterion of the $ \chi^2 $ test can be explicitly linked to a probability, providing a valuable reference for our study.

1. AIC: Since the diffusion model and the alternative models are non-nested, the likelihood ratio test and Wilks' theorem cannot be directly applied to assess the statistical significance of the differences between them. Therefore, we employ the AIC as a model selection criterion, which is widely used in data analysis for comparing non-nested models. The AIC is defined as

where k is the number of free parameters in the model and $ \mathcal{L}_{\max} $ is the maximum likelihood. To compare different spatial hypotheses, we define

where $ \mathrm{AIC}_{\rm diffusion} $ corresponds to the diffusion model. Smaller values of $ \Delta \mathrm{AIC} $ indicate support for the alternative model. We adopt $ \Delta \mathrm{AIC} \gt 6 $ as the criterion for model discrimination, as this threshold represents strong evidence [23, 24] in favor of the diffusion model.

2. $ \chi^2 $ test: The $ \chi^2 $ statistic provides a straightforward and intuitive way to assess the compatibility between different spatial hypotheses and the mock data. Specifically, the $ \chi^2 $ test is performed on one-dimensional (1D) radial brightness profiles, with bin widths chosen to ensure at least 10 counts per bin. This threshold allows Poisson fluctuations to be reliably approximated by a Gaussian distribution. For a given spatial template, the $ \chi^2 $ value is calculated as:

where $ O_i $ represents the observed excess in the i-th bin, $ T_i $ is the corresponding model prediction, and $ \sigma_i $ denotes the statistical uncertainty. When alternative spatial models, such as Gaussian or disk profiles, are fitted to the 1D radial brightness profiles generated by the diffusion model, they generally yield higher $ \chi^2 $ values than the diffusion model itself, as shown in Fig. 1. As the statistical quality of the data improves, these discrepancies become more pronounced, allowing incorrect models to be excluded based on their reduced goodness of fit. We adopt a $ 2\sigma $ exclusion criterion, corresponding to a p-value of approximately $ 0.05 $ for the $ \chi^2 $ distribution with $ n-k $ degrees of freedom, where n is the number of spatial bins and k is the number of free parameters in the fitted model. Models resulting in a $ \chi^2 $ value exceeding the critical threshold at this significance level are considered statistically distinguishable from the diffusion scenario.

Figure 1.  (color online) Radial γ-ray brightness profile of a pulsar halo using the Geminga parameters. The black points represent mock profile data for CTA-North ($ 50 $ h) generated using the diffusion model. The red, blue, and green curves show the best-fit diffusion, Gaussian, and disk models, respectively.

To mitigate the impact of statistical fluctuations arising from a single Monte Carlo realization, we perform 1,000 independent toyMC simulations for each parameter configuration and evaluate the corresponding statistical power.

We establish a statistical power threshold of $ 0.8 $, which requires that at least 800 out of 1,000 realizations satisfy the discrimination criterion (either the $ \Delta \mathrm{AIC} \gt 6 $ condition or the $ 2\sigma $ $ \chi^2 $ exclusion). When a specific combination of d and $ \dot{E} $ yields a statistical power of $ 0.8 $, this indicates that the diffusion signature of the pulsar halo for these parameters can be resolved by the experiment under consideration. Consequently, we define a critical curve in the $ d-\dot{E} $ plane, above which lies the parameter space where the diffusion signature can be distinguished.

The results of the morphological discrimination analysis are summarized in Fig. 2, showing the critical curves in the $ d-\dot{E} $ plane introduced in Section IV C. The four subfigures correspond to different combinations of statistical criteria and competing spatial templates: the AIC and the $ \chi^{2} $ test, each applied to compare the diffusion-based halo model with Gaussian or disk spatial templates. Each subfigure displays four discrimination curves representing different experimental configurations: LHAASO-KM2A (blue solid line), LHAASO-KM2A with a PSF improved by 40% (blue dashed line), CTA-North with a 50 h exposure (red solid line), and CTA-North with a 200 h exposure (red dashed line). Notably, for LHAASO-KM2A, the 1-year exposure time adopted in this study refers to the cumulative effective exposure time rather than the calendar duration. Given that a source at $ \delta = 30^\circ $ has an effective daily observation window of $ 7.8\,\mathrm{h} $ ($ \theta_{\mathrm{zen}} \lt 50^\circ $), this 1-year effective exposure corresponds to a total observational period of approximately 3 calendar years.

Figure 2.  (color online) Discrimination curves for pulsar halos achievable with CTA-North and LHAASO-KM2A. We use $ \dot{E}/\dot{E}_G $ instead of $ \dot{E} $ on the vertical axis, where $ \dot{E}_G $ denotes the luminosity of the reference source, Geminga. Sources located above the curves can be morphologically distinguished as diffusion halos. Filled circles denote confirmed pulsar halos or previously reported candidates, while open circles denote other pulsar-associated sources in the LHAASO catalog. The competing spatial models and statistical decision methods used in the four panels, from the top left to the bottom right, are as follows: a Gaussian model with the AIC, a disk model with the AIC, a Gaussian model with the $ \chi^2 $ test, and a disk model with the $ \chi^2 $ test.

The positions of confirmed or potential pulsar halos are overlaid for reference, including the Geminga and Monogem halos [4, 15], LHAASO J0621+3755 [25], LHAASO J0249+6022 [16], HESS J1831-098 [26, 27], and the source around PSR J0359+5414 [28]. The characteristic parameters of the pulsars associated with these candidates are summarized in Table 1. In addition, we consider pulsar-associated sources from the LHAASO catalog with characteristic ages greater than $ 20\; \mathrm{kyr} $, excluding younger systems that may still be dominated by emission from their initial pulsar wind nebulae. Halo candidates located outside the LHAASO field of view are not shown [29, 30]. Note that the performance metrics for LHAASO-KM2A presented here are derived under optimal observation conditions, effectively representing the sensitivity for sources at $ \delta \approx 30^\circ $. For sources at other declinations, the detection performance decreases, requiring longer exposure times to achieve comparable statistical significance. For instance, a source such as HESS J1831-098, located at a less favorable declination, would require at least twice the effective observation time to reach the same level of sensitivity.

A general trend is that the AIC-based method yields less restrictive discrimination conditions than the $ \chi^{2} $ test. For the same spatial template and instrumental configuration, the AIC identifies a larger region of the d–$ \dot{E} $ parameter space as distinguishable, allowing more halo candidates to be resolved by both CTA and LHAASO-KM2A. This behavior is expected, because the AIC is designed for relative model comparison between non-nested hypotheses by explicitly accounting for differences in model complexity, whereas the $ \chi^{2} $ test primarily evaluates the absolute goodness of fit of a given model based on 1D profiles and is therefore more conservative.

Comparing the results obtained with different competing spatial templates, it is apparent that distinguishing diffusion templates from disk templates is easier than differentiating them from Gaussian ones. This reflects the fact that disk models, characterized by flat central areas and sharp boundaries, differ more significantly from the smooth radial profiles produced by diffusion. On the other hand, a Gaussian template may partially mimic the central-peaked signature of diffusion halos, especially under conditions of restricted angular resolution or photon statistics.

The complementary performance of CTA and LHAASO-KM2A as a function of pulsar distance is clearly demonstrated. For relatively nearby objects ($ d\lesssim1.5 $ kpc), where the intrinsic angular extension is large, both experiments are in principle capable of resolving the morphology. In this regime, the discrimination power is mainly driven by photon statistics rather than angular resolution. Owing to its nearly continuous operation, LHAASO-KM2A accumulates substantially higher event statistics, resulting in a stronger discrimination capability than CTA for relatively nearby halos. For more distant pulsar halos, the apparent angular size becomes smaller, and the discrimination power is increasingly limited by the PSF. In this regime, CTA outperforms LHAASO-KM2A despite its shorter exposure time, owing to its superior angular resolution. Although LHAASO-KM2A benefits from higher photon statistics, its broader PSF limits its ability to resolve subtle morphological differences at large distances. Additionally, CTA's effectiveness for nearby halos is constrained by its limited field of view (FoV). For extremely extended sources like the Geminga halo, off-axis acceptance and the lack of gamma-ray-free regions for background estimation introduce significant challenges. Thus, our CTA forecasts are more reliable for compact or distant halos, while LHAASO-KM2A remains superior for nearby objects with large angular extensions owing to its wider FoV.

The impact of enhanced instrumental performance is also clearly demonstrated. Improving the angular resolution of LHAASO-KM2A yields only marginal gains for nearby sources, as their morphology is already well resolved and the available photon statistics remain unchanged. In contrast, at larger distances, the improved PSF substantially extends the parameter space within which halos can be discerned, enabling the detection of sources that would otherwise remain unresolved. Increasing the CTA exposure time from 50 h to $ 200 $ h systematically enhances discrimination capability across both nearby and distant sources, because its angular resolution is already sufficient and the additional exposure primarily boosts the photon statistics.

Consistent with current observational results, LHAASO-KM2A can presently perform unambiguous morphological discrimination only for the nearby Geminga and Monogem halos. However, with a 40% improvement in angular resolution, LHAASO-KM2A would be able to identify several major pulsar halo candidates, namely the γ-ray halos surrounding pulsars J1831-0952, J0248+6021, and J0359+5414. With its nominal design performance and a typical exposure of 50 h, CTA is already capable of distinguishing the morphology of most candidate pulsar halos. The most challenging among the current candidates is the halo associated with LHAASO J0621+3755, because its pulsar has a spin-down luminosity comparable to that of Geminga but lies at a distance approximately six times greater. With an increased exposure of 200 h, CTA is expected to achieve morphological discrimination between the diffusion and disk templates for LHAASO J0621+3755. However, distinguishing the diffusion template from a Gaussian model would require even greater photon statistics. Furthermore, a higher pair-conversion efficiency from spin-down power leads to increased data accumulation, thereby enhancing the experimental resolving power.

Additionally, observational evidence for deviations from spherical symmetry in pulsar halos has been steadily accumulating. Specifically, an asymmetric morphology for the Geminga halo was reported in the ICRC 2023 proceedings [31], followed by a 2024 HAWC measurement [15] suggesting a marginal asymmetry at approximately the $ 2\sigma $ level. Furthermore, a significantly elongated γ-ray structure has recently been identified near the millisecond pulsar MSP J0218+4232 [32]. On the theoretical front, such asymmetries are predicted to arise from anisotropic diffusion aligned with local magnetic fields [33, 34] or spatially varying diffusion coefficients induced by the local environment [35]. To assess how these complex morphologies might influence our results, we performed a sensitivity test in Appendix A using mock data generated from an asymmetric template. Our findings indicate that the discrimination power actually increases when the underlying source is asymmetric, primarily due to the more pronounced structural mismatch between the asymmetric diffusion morphology and the symmetric Gaussian assumption used in the fit. Therefore, the possible asymmetry of pulsar halos does not necessarily undermine the power of the morphological discrimination, demonstrating that our method is robust.

We have investigated the capability of CTA and LHAASO-KM2A to identify pulsar halo candidates through morphological discrimination. Using realistic mock observations derived from diffusion-driven halo templates and instrument response functions, we compared the diffusion morphology against alternative, non-nested spatial templates, including Gaussian and disk models, employing both the AIC and the $ \chi^{2} $ test for statistical assessment. The AIC generally provides more optimistic results than the $ \chi^{2} $ test, allowing a larger fraction of halos to be distinguished from simplified morphologies. Diffusion halos are more easily separated from disk templates than from Gaussian ones, reflecting a stronger morphological contrast.

A clear complementarity between CTA and LHAASO-KM2A is evident. For relatively nearby halos with large angular extensions, LHAASO-KM2A benefits from the higher photon statistics enabled by its large effective area and continuous operation. For more distant halos with smaller apparent sizes, CTA achieves superior discrimination owing to its significantly finer angular resolution. Furthermore, improving the angular resolution can significantly enhance LHAASO-KM2A's ability to identify distant pulsar halos, such as LHAASO J0249+6022, HESS J1831-098, and the candidate halo around PSR J0359+5414. For CTA, adopting longer observation times could, in principle, resolve the diffusion signatures of all candidate pulsar halos within 5 kpc, including LHAASO J0621+3755.

The discussion in this work is based on the ideal scenario in which the object is an isolated source. If the target is contaminated by other nearby γ-ray sources, morphological discrimination would be more challenging. Furthermore, identifying diffusion signatures is a necessary but insufficient condition for confirming a pulsar halo. Other criteria for confirming γ-ray pulsar halos include the offset between the halo and the pulsar, the relationship between the energy required to produce the halo and the pulsar spin-down power, the relative sizes of the γ-ray halo and the corresponding X-ray PWN [5], and a comparison between the electron energy density required to explain the halo and the typical energy density of the ISM [36]. In practice, multi-wavelength observations are also required to determine whether a pulsar halo is the most plausible interpretation.

Observational evidence for deviations from spherical symmetry in pulsar halos has been steadily accumulating. Specifically, an asymmetric morphology for the Geminga halo was reported in the ICRC 2023 proceedings [31], and a subsequent 2024 HAWC measurement [15] suggested a marginal asymmetry at approximately the $ 2\sigma $ level. Furthermore, a significantly elongated γ-ray structure has recently been identified near the millisecond pulsar MSP J0218+4232 [32]. Moreover, theoretical studies have proposed that such asymmetries can naturally arise from anisotropic diffusion due to the alignment of local magnetic fields [33, 34] or from spatially varying diffusion coefficients induced by the local environment [35]. To assess the impact of such asymmetries on our findings, we repeated our analysis using mock halo signals with asymmetric spatial templates to verify the robustness of our results.

The asymmetric halo morphology we adopted follows recent observations by HAWC [15, 37]. In this model, the characteristic diffusion extension $ \theta_d $ varies with azimuth. Consistent with the HAWC measurement, the halo is divided into four right-angled spherical sectors centered on the pulsar. In each sector, the γ-ray surface brightness is still described by the form of Eq. (1), with $ \theta_d $ taking different values in the four sectors. To ensure model consistency at the center and the proper normalization of the morphological component, additional multiplicative factors $ \lambda_i \, (i=1,2,3,4) $ are applied. The factors $ \lambda_i $ are constrained by the following relation:

We obtain the expression in Eq. 6.

For consistency, a symmetric diffusion model is still employed alongside a symmetric Gaussian template to fit the mock data. $ \Delta \mathrm{AIC} \gt 6 $ is adopted as the criterion here. The new discrimination curves are presented in Fig. A1, indicating stronger discrimination capabilities for both LHAASO-KM2A and CTA, owing to the more pronounced deviation of the symmetric Gaussian template from the mock data, which were generated using an asymmetric diffusion model. The assumption of an asymmetric halo does not affect the main conclusions drawn in the main text, and the complementary performance between CTA and LHAASO-KM2A remains.

Figure A1.  (color online) Morphological discrimination curves for pulsar halos based on both symmetric and asymmetric diffusion models. These curves represent the sensitivities of CTA-North and LHAASO-KM2A. The analysis employs a Gaussian spatial template as the competing model, using the AIC as the statistical decision metric.