pulsar or generic particle transport behaviors. The latter scenario predicts that TeV halos would be a universal phenomena for all pulsars. We searched for extended gamma-ray emission around 36 isolated middle-aged pulsars identified by radio and gamma-ray facilities using 2321 days of data from the High-Altitude Water Cherenkov (HAWC) Observatory. Through a stacking analysis comparing TeV flux models against a background-only hypothesis, we identified TeV halolike emission at a significance level of 5.10σ. Our results imply that extended TeV gamma-ray halos may commonly exist around middle-aged pulsars. This reveals a previously unknown feature about pulsars and opens a new window to identify the pulsar population that is invisible to radio, x-ray, and GeV gamma-ray observations. DOI: 10.1103/PhysRevLett. 134.171005 Gamma-ray emission with angular size of a few degrees was first detected by the High-Altitude Water Cherenkov (HAWC) Observatory around two nearby middle-aged energetic pulsars between 8 and 40 TeV in 2017 [1]. The TeV emission seen around those pulsars is much more extended than the pulsar wind nebulae whose sizes are no more than a few arc minutes [2][3][4][5]. These extended sources were named TeV halos [6]. Since the discovery of the first two TeV halos, more similar sources have been identified by HAWC [7] and LHAASO [8]. To date, the online catalog for TeV Astronomy (TeVCat [9]) reports a total of about ten TeV halos and TeV halo candidates. Counterparts of TeV halos at lower energies have been found in Fermi-LAT [5,10] and H.E.S.S. [4] data.

The origin of the TeV halos is largely unknown [11]. The emission cannot be explained by electron diffusion in the interstellar medium (ISM) as it requires a suppression of the cosmic-ray diffusivity by 100-1000 times [1]. Several explanations have been proposed, but all have weaknesses. Models can be divided into two main classes, attributing the slow diffusion to a turbulence that universally exists around the central source [12,13] or preexisting regions of nonstandard relativistic electron transport [14,15]. A related debate is whether the halos commonly exist around middleaged pulsars [16] or most middle-aged pulsars do not develop halos [17].

The several known TeV halos present similar physical extensions [6], which are related to the cooling time of very-high-energy (0.1-100 TeV) electrons in the cosmic microwave background (CMB) and interstellar radiation field. The flux of the TeV emission seems to be related to the pulsar's emissivity [18]. Motivated by the similarity of the known sources, we perform the first systematic search for TeV halos in the HAWC data around pulsar populations that are identified by radio and GeV gamma-ray observations.

Samples of isolated middle-aged pulsars-Our analysis uses radio pulsars from the Australia Telescope National Facility (ATNF) catalog [19,20] and gamma-ray pulsars from the Third LAT Pulsar catalog (3PC) [21]. We apply the following selection criteria: (1) a pulsar is in HAWC's field of view (-26°< δ < 64°), (2) not in a binary system, (3) middle aged, defined as the spin down age P=ð2 ṖÞ > 20 kyr, where P and Ṗ are the barycentric period of the pulsar and its time derivative, respectively, and (4) with the spin-down flux above the HAWC sensitivity. The spin-down flux is defined as Ė=ð4πd 2 Þ, where Ė and d are the spin-down power and distance, respectively. A total of 1304 sources pass the first three selection criteria. Since the distance to a pulsar is needed to estimate its detectability and the angular size of its halo, we remove 3 pulsars without reported distances from the list. After applying the fourth criterion, the selection yields a total of 86 sources, including 49 gamma-ray pulsars and 72 radio-loud pulsars. Among the gamma-ray pulsars, 14 are radio quiet and mutually exclusive from the radio-loud sample.

Some of the sources are located in highly populated regions of TeV gamma-ray sources. To avoid source confusion, we further removed pulsars within 2°:5 of the TeVCat sources that are not classified as a TeV halo or TeV halo candidate. The choice of 2°∶5 makes sure that the contamination from most TeV sources is avoided while the sample remains sufficiently large. Our final sample has 36 isolated middle-aged pulsars, including 24 gamma-ray pulsars and 24 radio-loud pulsars (12 pulsars are both radio-loud and gamma-ray bright). Among them, 28 pulsars have not been associated with a known TeV halo or TeV halo candidate. Figure 1 shows a HAWC significance map of the entire sky above 300 GeV in Galactic coordinates, overlaying the positions of the pulsars used in the analysis. Detailed properties of the pulsars are presented in Supplemental Material [22].

HAWC data-HAWC is a gamma-ray and cosmic-ray observatory located near Puebla, Mexico at an altitude of 4,100 meters. This analysis uses the HAWC Data Pass 5 [23] with 2321 days of exposure. The gamma-ray energy is estimated using the fraction of the photomultiplier tubes (PMTs) that are triggered in each shower event, called f Hit bins [24].

Stacking analysis-Motivated by the similarity of the observed TeV halos, we combine the observations of the pulsar population accessible to HAWC to constrain the common factors that impact the gamma-ray production in TeV halos.

We model the TeV halo emission of the ith pulsar in the source list with a spatial model and a spectral model. The spatial model is a Gaussian template with an extension σ PHYSICAL REVIEW LETTERS 134, 171005 (2025) scaled to the extension of the Geminga TeV halo, σ i ¼ σ Geminga ðd Geminga =d i Þ, where σ Geminga ¼ 2° [6], d Geminga ¼ 0.25 kpc and d i is the distance to the pulsar in kpc. The value of σ i is fixed in the fit. The spectral model describes the flux as a simple power law, dN=dEdAdt ∝ ðE=E piv Þ -α , where E piv is the pivot energy and α is the spectral index. We fit such a model to data in both the full energy range, 0.316 TeV-100 TeV, and in half-decade bins. In the full energy range the pivot energy is fixed to 30 TeV, while in the fit in each energy bin, the pivot energy is chosen to be the median energy of each energy bin. We fix the spectral index to α ¼ 2.7 and evaluate the uncertainties associated with the choices of spatial extension and spectral index later in a systematic study.

We assume that all TeV halos in the population may be described by the same physics, where the differential energy flux of a TeV halo, F i ≡ E 2 ðdN=dEdAdtÞ i , is scaled to the available power of a pulsar, A i , through a common efficiency, η: F i ¼ ηA i . The efficiency η represents the fraction of the pulsar's power that is converted into the differential TeV halo luminosity at a given energy. We consider two scenarios for A i , namely, the spin-down flux scenario and the GeV flux scenario. In the former model,

i is the current-day spin-down power with I ¼ 10 45 g cm -2 being the moment of inertia of a neutron star which is fixed to be the same for all pulsars. In this scenario, the closer and more energetic pulsars have a higher gamma-ray flux. In the GeV flux scenario, we assume A i ¼ F 0.1-100 GeV , where F 0.1-100 GeV is the integrated flux of the pulsar between 100 MeV and 100 GeV provided by 3PC. In this second scenario, brighter GeV pulsars would produce brighter very-high energy (VHE) gamma-ray halos. Independent of the distance, the spindown flux and the GeV flux are correlated, though the slope has a large dispersion as shown in Ref. [21]. We therefore treat them as two separate scenarios in the stacking analysis.

We fit each source in the pulsar samples using the spectral and spatial model described above. We then add their log-likelihood profiles to constrain the common efficiency η in each power scenario and energy bin.

We compare the combined signal with a null hypothesis, where the flux of extended TeV halo emission around pulsars is consistent with background fluctuations. We first construct a test statistic (TS), which is defined to be twice the logarithm of the ratio of the likelihoods when fitting the data with and without TeV halos around the selected pulsars, TS ≡ 2 ln½LðηÞ=Lðη ¼ 0Þ. To understand the TS distribution of the null hypothesis, we perform a Monte Carlo simulation to generate random positions based on the spatial distribution of Galactic pulsars (see Supplemental Material for details regarding the background source generation). These positions are treated as fake sources and we stack their likelihood profiles in the actual data. We use the same selection criteria as for the real sources to form a fake source sample, that is, a fake source needs to be in the HAWC field of view and located at least 2°:5 away from TeVCat sources. The black curve in the left panel of Fig. 2 presents the distribution of the TS under the null hypothesis. Out of 101 965 trials of stackings of 10 fake sources, the highest TS is found to be TS bg ¼ 69.5. Therefore, a TS much higher than TS bg can safely reject the background hypothesis at least at the level of 4.27σ. We verified that the background TS distribution is not sensitively impacted by the number of stacked sources (see Supplemental Material). The background TS distribution departs from a Chi-square over two (χ 2 =2) distribution that is expected from Wilks' theorem [25]. This is caused by the fake sources that land at the outskirts of TeV sources that are more extended than 2°:5, the separation threshold adopted in our source selection criteria.

Results-We perform the stacking analysis using radio and gamma-ray pulsar samples and the full energy range. We adopt the spin-down flux scenario for radio pulsars, considering that not all radio pulsars are detected in GeV gamma rays. For gamma-ray pulsars, we consider both weighting scenarios. This led to a total of 3 trials of stacking analyses. To avoid bias caused by bright TeV halos, we first consider the subsample that is not associated with known TeV halo or TeV halo candidates, which includes 17 gamma-ray and 20 radio-loud pulsars. The highest TS is found to be TS ¼ 108, obtained from the trial using gamma-ray pulsars and the spin-down flux scenario. Based on an extrapolation of the simulated background TS distribution (gray curve in left panel of Fig. 2), this TS value corresponds to a p value p ¼ 5.5 × 10 -8 , based on an extrapolation of the TS distribution of the background trials.

The other two trials result in TS ¼ 92 and 96, respectively, as also indicated in the figure. We apply a penalty trial factor to account for the fact that multiple trials may increase the possibility of finding a significant event.

Although the three signal hypotheses are correlated, we account for a conservative trial factor of 3 when evaluating a post-trials p value. Specifically, we multiply the lowest p value by 3 to obtain a post-trial p value and convert the p value to a Gaussian significance of 5.10σ. This positive detection of extended TeV gamma-ray emission around pulsar populations not previously associated with TeV To ensure that the stacked TS is not dominantly contributed by a few bright sources, we gradually remove the sources with the highest TS values from a pulsar sample and examine the evolution of the stacked TS values. We use all pulsars including those associated with known TeV halos in this test. The right panel of Fig. 2 shows that in all weighting scenarios and pulsar groups, the stacked TS gradually declines as the number of stacked sources decreases. The trend cannot be explained by the background which is indicated by the gray shaded region. The evidence of an excess halo emission against the background comes from the entire pulsar population as opposed to just the brightest sources.

Finally, to examine whether the emission is physically extended, we perform the stacking analysis using a pointlike spatial template. Similar to the extended-model analysis, we confine the efficiency parameter using a combined likelihood. We compare the two spatial models based on the TS values and the Bayesian information criterion (BIC) from the analysis using the full energy range. We find that in all weighting scenarios and pulsar groups, the extended model is significantly preferred [26] over the pointlike model with TS extended -TS point > 25 and BIC point -BIC extended > 10.

The halo efficiencies computed using all pulsars are generally consistent with those from pulsars not associated with known TeV halos within uncertainties as shown in Fig. 6 in Supplemental Material. The systematic uncertainties include the effects due to spectral models, extension models, as well as uncertainties in our detector configurations. See Supplemental Material for details regarding the evaluation of the systematic errors. The efficiencies decrease as a function of the photon energy because the gamma-ray spectrum follows a power-law distribution. On average, the differential energy flux in TeV halos is at the level of (0.01-1)% of the spin-down flux. For gamma-ray pulsars, the energy deposited in TeV halos ranges from (0.1-10)% of that in the GeV emission of the pulsars.

Implications of commonly existing TeV halos-We have used a stacking technique that combines the likelihood profiles of radio and gamma-ray pulsars accessible to HAWC to test the hypothesis that the flux of extended TeV halo emission around pulsars is consistent with background fluctuations. Our analysis has assumed that the TeV flux of a pulsar halo scales with the spin-down flux or the GeV gamma-ray flux of the pulsar. Using pulsars that have not been associated with a known TeV halo, the null hypothesis is rejected at the level of 5.10σ.

Our analysis implies an existence of ubiquitous TeV halos around isolated middle-aged pulsars. This establishes TeV halos as a new pulsar phenomenon that has a distinct origin from pulsar wind nebulae, as the latter are usually orders of magnitude more compact for pulsars at this age.

A common existence of TeV halos suggests that the production mechanism of TeV halos is more likely related to the transport of relativistic particles in the vicinity of a pulsar than the local environment of the observed halos. It also suggests a leptonic origin of the TeV gamma-ray emission since hadronic emission highly depends on the gas density (e.g., Ref. [27]) and would vary from pulsar to pulsar.

The longer confinement of relativistic electrons by the halos than the average ISM has strong implications for indirect searches for dark matter. While these halos provide an alternative way to explain the GeV excess in the Galactic center [28], they also make it harder to explain the cosmicray positron excess [29][30][31] with nearby pulsars [1,27,32], compared to scenarios involving exotic physics. A common existence of TeV halos also impacts the view of Galactic gamma-ray astronomy. These halos may contribute to the diffuse Galactic plane emission observed by Tibet ASγ [33], HAWC [34], and LHAASO [35] experiments and explain the hardening of the diffuse γ-ray emission around 1 TeV [36][37][38][39].

Our Letter suggests that a conversion of the spin-down power to the TeV halo emission could be as efficient as ∼0.1% and that the power carried by a TeV halo is about 1% of that by the GeV magnetospheric emission at 10 TeV (see Fig. 6 in Supplemental Material). The spin-down flux model resulted in a higher TS than the GeV flux model, though the difference is not statistically significant. This could be related to the fact that the observation of GeV magnetospheric emission may be biased by orientation effects. We further extended this analysis to constrain the diffusion radius of electrons around our sample of TeV halos. From this, we derived a diffusion coefficient of D 0 ≈ 2.0 × 10 27 cm 2 s -1 at a reference energy of E 0 ¼ 10 GeV, consistent with measurements around Geminga [1] (see Supplemental Material and Table IV for more details regarding the diffusion coefficient study). These results reinforce the idea that electron propagation in TeV halos occurs in a suppressed diffusion environment, significantly lower than the typical Galactic diffusion coefficient.

Future observations of the spatial profiles of individual halos over a wide energy range would be needed to understand how electrons diffuse in TeV halos. Wide-field gamma-ray detectors with better sensitivity and covering the southern sky can carry out such population studies [40]. Finally, while pulsars are typically identified via their pulsed emission in radio, x-ray, and GeV gamma-ray wavelengths, TeV halos open a new window on pulsar observations. In particular, they offer a unique way to probe the "invisible" pulsars [6] with beaming angles misaligned with observers and middle-aged and old pulsars with a surface density too low to be detected in other wavelengths.