© The Authors 2026

1. Introduction

Extended gamma-ray emission around pulsars has emerged as a phenomenon distinct from compact pulsar wind nebulae (PWNe), which are typically confined within a few parsecs (Gaensler & Slane 2006; Torres et al. 2014). These extended sources, termed pulsar halos (or sometimes teraelectronvolt halos), span degree scales without clear boundaries; they were first revealed by the High Altitude Water Cherenkov (HAWC) around Geminga and Monogem (Abeysekara et al. 2017; Albert et al. 2024), and were interpreted as relativistic e^± diffusing into the interstellar medium (ISM) and producing gamma-ray emission via inverse-Compton (IC) scattering of the background radiation field (Sudoh et al. 2019). The bright, extended teraelectronvolt (TeV) morphologies of these halos (tens of parsecs) indicate suppressed diffusion coefficients (Abeysekara et al. 2017) relative to the Galactic average extrapolated from the boron-to-carbon (B/C) ratio (Génolini et al. 2019; De Sarkar et al. 2021). However, alternative explanations are available (see, e.g., Liu et al. 2019; Aloisio et al. 2009, with detailed reviews in Fang 2022; Liu 2022).

These halos may help explain the long-standing cosmic-ray positron excess problem observed by Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (PAMELA), Alpha Magnetic Spectrometer (AMS-02), and Fermi Large Area Telescope (LAT) (Adriani et al. 2009; Aguilar et al. 2013; Ackermann et al. 2012), as nearby sources like Geminga and Monogem could dominate the local e^± flux (Hooper et al. 2017; Manconi et al. 2020), though dark matter annihilation, and secondaries produced in supernova remnants (SNRs) or molecular clouds (MCs) remain viable scenarios (Hooper & Kribs 2004; Cholis et al. 2009; Cholis & Hooper 2013; Mertsch & Sarkar 2014; De Sarkar et al. 2021). The presence of halos also suggests that many unidentified extended TeV sources are halos rather than shell-type SNRs or PWNe (Giacinti et al. 2020; Martin et al. 2022a,b; López-Coto & Giacinti 2018), as halos are potentially common around middle-aged (> 20 kyr) pulsars (Albert et al. 2025). It is indeed a possibility that new halo candidates will continue to emerge in Fermi-LAT (Di Mauro et al. 2019, 2021), High Energy Stereoscopic System (H.E.S.S.) (H.E.S.S. Collaboration 2018), HAWC (Albert et al. 2020), and Large High Altitude Air Shower Observatory (LHAASO) catalogs (Cao et al. 2024). LHAASO detections above 100 TeV further highlight the prevalence of Galactic “PeVatrons”, traditionally attributed to SNR-MC systems or PWNe (see, e.g., Kar & Gupta 2022; De Sarkar & Gupta 2022; De Sarkar et al. 2022; De Sarkar 2023; De Sarkar & Majumdar 2024; Joshi et al. 2023; Wu et al. 2023; Mitchell 2024; Mitchell & Celli 2024; Abe et al. 2023; de la Fuente et al. 2023), but pulsar halos themselves may also qualify if source confusion in the crowded inner Galaxy can be resolved (Mestre et al. 2022).

Within this context, PSR J1813-1246 and its associated TeV source HESS J1813-126 offer a clean case study due to their off-plane location (b ≈ 2.5°). PSR J1813-1246 was discovered with Fermi-LAT (Abdo et al. 2009) and is the fastest-spinning, second-most energetic radio-quiet pulsar, with period P ≃ 48.1 ms, period derivative Ṗ ≃ 1.76 × 10⁻¹⁴ s s⁻¹, current spin-down energy-loss rate Ė ≃ 6.24 × 10³⁶ erg s⁻¹, and characteristic age τ_c ∼ 43 kyr (Abdo et al. 2013). The distance to the pulsar was estimated to be between a minimum of 2.5 kpc and a maximum of 24.7 kpc (Abdo et al. 2013; Marelli et al. 2014). Despite its energetics, no radio pulsations have been detected (Ray et al. 2011; Abdo et al. 2013). In X-rays, XMM-Newton and Chandra revealed a pulsed counterpart (Marelli et al. 2014). A possible detection of the pulsar in the 30–500 keV range was also reported using the International Gamma-Ray Astrophysics Laboratory (INTEGRAL) data, which is consistent with the extrapolation between the soft X-ray and gamma-ray spectra (Marelli et al. 2014).

The pulsar is positionally coincident with HESS J1813-126, a hard-spectrum (E⁻²), extended (∼0.21°) TeV source from the H.E.S.S. Galactic Plane Survey (HGPS; H.E.S.S. Collaboration 2018). No gigaelectronvolt (GeV) or X-ray PWN associated with the pulsar has been detected (Ackermann et al. 2011; Marelli et al. 2014; Guevel et al. 2025), thus supporting a halo origin. HAWC detected the source as 3HWC J1813-125 (Albert et al. 2020), while LHAASO detected > 100 TeV emission from the region as 1LHAASO J1813-1245 (Cao et al. 2024), confirming its extended morphology (≲0.32°) and PeVatron potential. Together, these data establish HESS J1813-126 as a bright, spatially extended TeV gamma-ray emitter extending into the ultra-high-energy regime and naturally linked to PSR J1813-1246, thus making it a possible pulsar halo candidate.

In this paper we study the combined PSR J1813-1246 – HESS J1813-126 system. Section 2 presents and discusses the results of this work, and Sect. 3 concludes. Appendix A presents the details of the Fermi-LAT standard and phased analyses of the GeV gamma-ray data from the pulsar and the search for a potential off-pulse emission. Appendix B discusses the synchro-curvature (SC) framework used to explain the high-energy emission of the pulsar, and Appendix C provides the formalisms of pulsar halo particle transport under the two-zone isotropic suppressed diffusion, ballistic-to-diffusion, and anisotropic diffusion scenarios.

2. Results and discussion

The method for the standard Fermi-LAT data analysis is discussed in Appendix A. PSR J1813-1246 is detected as 4FGL J1813.4-1246 in the updated Fourth Fermi Gamma-ray LAT Source Catalog – Data Release 4 (4FGL-DR4). The best-fit localization of the GeV counterpart was found at RA = 273.3471° ±0.0025° and Dec = −12.7667° ±0.0025°, with a detection significance of ∼136σ. The best-fit position of the GeV counterpart is offset from the position of the pulsar by 0.08°, and is consistent with the centroid of the HESS J1813-126 (RA = 273.340° and Dec = −12.688°), as well as being situated well within the containment radii of the HAWC and LHAASO detected counterparts. The spectral analysis of the 4FGL source reveals a significantly curved spectral energy distribution (SED). The spectrum is best described by a PLSuperExpCutoff4 function¹, which is favored over a simple power law function by 2log(ℒ_PLSEC4/ℒ_PL)≈1140.2, and log-parabola function by 2log(ℒ_PLSEC4/ℒ_LP)≈982.8. This function is particularly relevant for modeling the gamma-ray emission of pulsars. The best-fit parameters of the PLSuperExpCutoff4 function for this source are γ_{0,  PLSEC4} = −2.459 ± 0.001, E_0, PLSEC4 = 1440 MeV, d_PLSEC4 = 0.667 ± 0.001 and b_PLSEC4 = 0.426, where γ_{0,  PLSEC4} and d_PLSEC4 are the spectral index and local curvature at reference energy E_{0,  PLSEC4}. The index b_PLSEC4 and reference energy E_{0,  PLSEC4} were kept fixed at the 4FGL-DR4 catalog values. The integrated energy flux between 1 and 500 GeV is (8.304 ± 0.095)×10⁻¹¹ erg cm⁻² s⁻¹. We further examined the source extension (68% containment radius) using both RadialDisk and RadialGaussian spatial models. We found that RadialDisk model gives a best-fit extension of $0.{0515}_{-0.0109}^{+0.0098}$ deg with TS_ext ≈ 7.169 (∼2.7σ), where for RadialGaussian, these values are $0.{0474}_{-0.0099}^{+0.0092}$ deg with TS_ext ≈ 7.307 (∼2.7σ). The source spatial model clearly agrees with the point-like hypotheses and disfavors the GeV source being extended. The curved SED and point-like morphology therefore confirm that 4FGL J1813.4-1246 is indeed the phase-averaged GeV counterpart of PSR J1813-1246. Figure 1 shows the significance map of the phase-averaged 4FGL J1813.4-1246, along with H.E.S.S., HAWC, and LHAASO source extension. HESS J1813-126 significance map, as observed by H.E.S.S., is also shown in Fig. 1. Note that no previously undetected source was found within H.E.S.S. extent and in the vicinity of the 4FGL source through the source finding algorithm in the phase-averaged standard analysis.

Fig. 1. Significance maps of the HESS J1813-126 region. Left: Significance map of 4FGL J1813.4-1246 from the phase-averaged analysis. Middle: Same but for PS J1813.3-1246, detected during the off-pulse phase of the 4FGL source. Right: H.E.S.S. significance map. The color bar in each plot provides the significance level. The source extent for H.E.S.S., HAWC, and LHAASO (KM2A and WCDA) is shown as colored circles.

Next, we discuss the result from the pulsar phased analysis, the method of which has also been discussed in Appendix A. Before the analysis, we constructed weighted pulse phase histograms and applied the Bayesian blocks algorithm (Jackson et al. 2005; Scargle et al. 2013), which adaptively identifies significant changes in photon count rate. This yielded an off-pulse interval (0.74–0.99) where pulsar contamination should be minimal. Note that our off-peak definition closely matches that given in Abdo et al. (2013, 0.74–0.98); however, it is somewhat different in Ackermann et al. (2011, 0.72–0.84). Further note that Abdo et al. (2013) has already cautioned that systematic uncertainties in discerning magnetospheric and diffuse background emission may complicate the detection of faint off-peak GeV emission, and additionally, some residual low-altitude magnetospheric emission may still persist. Using a larger dataset, we reexamined the source to search for any residual off-peak emission that could potentially be crucial evidence in favor of HESS J1813-126 being a GeV PWN. The weighted pulse phase histogram, along with the off-pulse phase selections from this work and previous literature (Ackermann et al. 2011; Abdo et al. 2013), are shown in Fig. 2, where the model weight has been considered from Bruel (2019). Note that the narrower phase interval of 0.83–0.99, identified through the Bayesian blocks analysis, was not adopted as the off-pulse phase region, since the limited photon statistics in this interval were insufficient to yield a statistically significant off-pulse emission.

Fig. 2. Weighted pulse profile of PSR J1813-1246 (4FGL J1813.4-1246) over one rotational cycle (pulse phase 0-1). It was constructed using 100 uniform phase bins per period from events between 0.1 and 500 GeV with the model weight from Bruel (2019). The blue histogram shows the weighted photon counts per phase bin, while the solid red line indicates the Bayesian block decomposition of statistically significant structures in the weighted pulse profile. The red shaded region signifies the off-pulse interval estimated from this work, whereas the magenta and green shaded regions signify the same obtained in Abdo et al. (2013) and Ackermann et al. (2011), respectively.

During the off-pulse phase of the GeV counterpart of the pulsar, 4FGL J1813.4-1246, a new source, PS J1813.3-1246, was detected at a minor offset of 0.011° from the position of the pulsar. The source was detected at the best-fit localization of RA = 273.3459° ±0.0122° and Dec = −12.7781° ±0.0120° with a detection significance of 23.8σ. The PS source is consistent with the containment radii of the sources detected by H.E.S.S., HAWC, and LHAASO. The spectrum of the source is explained by a simple power law, which has best-fit parameters of α_PL = −3.123 ± 0.001 at a reference energy scale of E_0, PL = 1000 MeV. The energy flux integrated within 1-500 GeV was found to be (2.419 ± 0.122)×10⁻¹¹ erg cm⁻² s⁻¹. We further examined the spatial extension of the source and found that the RadialDisk model yields a best-fit extension of $0.{0872}_{-0.0316}^{+0.0236}\mathrm{deg}$ with TS_ext ≈ 2.8 (∼1.7σ), while the RadialGaussian model gives $0.{0752}_{-0.0292}^{+0.0268}\mathrm{deg}$ with TS_ext ≈ 2.4 (∼1.5σ), thus clearly hinting at a point-like morphology.

Given that the off-pulse source spectrum is adequately described by a simple power law and no significant extension is detected, the balance of evidence favors residual pulsar magnetospheric contamination rather than a GeV PWN. First, the source is point-like and positionally coincident with the pulsar, whereas a nebular IC component is expected to be spatially extended at GeV energies. Second, the measured power law index is relatively soft for a GeV PWN (∼1.8 − 2.2) and does not connect naturally to the TeV spectrum without introducing a break, but it is compatible with an unpulsed or tail magnetospheric component leaking into the off-pulse window due to imperfect gating or ephemeris uncertainties. Third, previous studies report no significant off-pulse GeV emission (Ackermann et al. 2011; Abdo et al. 2013; Marelli et al. 2014), consistent with our finding that any putative nebular signal would be subdominant. We therefore interpret the off-pulse SED as an upper limit and attribute the detected GeV flux to residual magnetospheric radiation rather than a GeV nebular emission. A future Fermi-LAT phased analysis with an updated ephemeris will be crucial to search for any GeV emission hidden beneath the pulsar signal. Figure 1 also shows the significance map of the off-pulse source PS J1813.3-1246.

Next, we move on to describe the gamma-ray and X-ray SEDs of PSR J1813-1246 with the SC framework. In this picture, relativistic e^± pairs are accelerated by a parallel electric field E_∥ in the outer pulsar magnetosphere and radiate while both spiraling around and sliding along curved magnetic field lines. The SC formalism naturally interpolates between the synchrotron and curvature limits and has been shown to reproduce the curved spectra observed by Fermi–LAT in many energetic pulsars. We briefly present the SC model equations in Appendix B (for a more detailed discussion, see, e.g., Cheng & Zhang 1996; Viganò et al. 2015a,b,c,d; Torres 2018; Torres et al. 2019; De Sarkar & Majumdar 2024). The gamma-ray SED corresponds to the Fermi-LAT source 4FGL J1813.4-1246 (this work), whereas the XMM-Newton X-ray spectrum and INTEGRAL X-ray upper limit are taken from Marelli et al. (2014) and Guevel et al. (2025). Assuming a minimum pulsar distance of r_s = 2.5 kpc, the SC parameters that reproduce the data are log(N₀) = 32.3, log(E_∥/V m⁻¹) = 9.25, magnetic field gradient b_sc = 3.0, and the scale length log(x₀/R_lc) =  − 3.4. The resulting model SED, together with the observational data, is shown in Fig. 3, along with the Lorentz factor Γ, the SC parameter ξ, and the pitch angle sin α as functions of the distance from the inner boundary x_in in units of the light-cylinder radius R_lc. In the SED plot of Fig. 3, the decomposed emission from different parts of the emitting region has also been shown, as a function of distance from x_in in units of R_lc. The model reproduces both bands well: the synchrotron component explains the X-ray SED, while the curvature component accounts for the GeV spectrum. The parameter ξ acts as a diagnostic of the dominant regime: if ξ ≪ 1, curvature radiation prevails, whereas ξ ≫ 1 signals synchrotron dominance. From Fig. 3 we see that the transition between the two regimes occurs at around (x − x_in)/R_lc ≈ 10⁻³, suggesting that the X-rays are produced at low altitudes where synchrotron losses dominate, while the gamma rays are emitted further out along the field lines, where curvature radiation becomes the main channel.

Fig. 3. Distance variation of different SC model parameters and the spectrum of PSR J1813-1246: Lorentz factor (Γ; first panel), SC parameter (ξ; second panel), and pitch angle (sin α; third panel). In the second panel, the dashed red line marks ξ = 1, i.e., the transition point between synchrotron- and curvature-dominated regimes. In the last panel, the model spectrum is plotted against the observed SED obtained from the Fermi-LAT data analysis of 4FGL J1813.4-1246 and the X-ray SED and upper limit from XMM-Newton and INTEGRAL data, respectively, as obtained from Marelli et al. (2014) and Guevel et al. (2025). In the last panel, decomposed emission from different parts of the emitting region is also shown; the colorbar signifies the distance from the inner boundary (xin) in units of the light-cylinder radius (Rlc).

It is to be noted that the origin of the observed quarter phase lag between the X-ray and gamma-ray light curves of PSR J1813-1246 (Marelli et al. 2014) remains a puzzling issue. Several possibilities could, in principle, account for such a lag, including emission from spatially distinct regions within the gap, such as a synchrotron-dominated inner zone and a curvature-dominated outer zone, and a presence of strong magnetic field sweepback near the light-cylinder (Contopoulos et al. 1999; Spitkovsky 2006; Kalapotharakos et al. 2012; Dyks & Rudak 2003; Dyks 2008; Dyks et al. 2004; Dyks & Harding 2004; Watters et al. 2009; Bai & Spitkovsky 2010), which can shift the formation of caustic peaks. However, quantitatively addressing these possibilities requires a fully self-consistent description of the magnetospheric electrodynamics, ideally within global resistive magnetohydrodynamic or kinetic (particle-in-cell) frameworks similar to that discussed in Marelli et al. (2014), or a more methodological approach (see, e.g., Viganò & Torres 2019; Íñiguez-Pascual et al. 2024) that incorporates realistic field topology, particle acceleration, and radiative transport in a coupled manner. Such an analysis falls beyond the scope of the present work, which is primarily focused on elucidating the pulsar halo origin of HESS J1813–126 and its connection to PSR J1813–1246. A detailed investigation of the phase-lag problem under the SC framework is therefore left to future dedicated studies.

Finally, we explored the pulsar halo origin of HESS J1813-126. To describe the transport of e^± pairs and to further discuss the SED and morphology of the resulting gamma-ray emission via the IC upscattering of the cosmic microwave background (CMB), infrared photons from interstellar dust, and optical stellar fields, in this section we consider three scenarios that have been proposed in the literature: the two-zone isotropic suppressed diffusion (2ZISD) model, the ballistic-to-diffusion (B2D) model, and the anisotropic diffusion (AD) model. These models differ mainly in how the spatial diffusion of particles is parameterized, but share the same general transport framework. We briefly discuss the formulae pertinent to the transport framework in Appendix C (for a more detailed discussion, see Fang 2025). We applied all of these three models to explore the pulsar halo origin of HESS J1813-126 using the PHECT code (Fang 2025).

In this paper we primarily focus on explaining the gamma-ray SED using IC emission, and do not delve into explaining the X-ray upper limit with synchrotron emission. Guevel et al. (2025) has reported no evidence of excess X-ray emission, and from the calculated X-ray upper limits (4.32 × 10⁻⁴ keV⁻¹ cm⁻² s⁻¹ and 5.38 × 10⁻⁴ keV⁻¹ cm⁻² s⁻¹ at 1 keV assuming E⁻² spectrum), they have posited that the magnetic field in the region should be in the range 3.5–6.5 μG, implying that the halo magnetic field is not significantly enhanced compared to the average Galactic magnetic field. In this work, we used a representative halo magnetic field of 5 μG to be consistent with Guevel et al. (2025). The current age of the system is unknown, so we considered it to be t_s = 30 kyr. We chose this value to be consistent with Albert et al. (2025), who posited that the typical age of pulsars that contain a pulsar halo is > 20 kyr. The initial spin-down timescale can be calculated as τ₀ = τ_c − t_s ≈ 13 kyr, given that we considered n = 3. The spectral index of the particle spectrum is chosen to be p = 1.6 and s = 1, making it a typical power law with exponential cutoff spectrum. Additionally, to calculate the IC emission, we considered radiation fields from CMB, interstellar dust, and starlight, with the parameters: T_CMB = 2.73 K, ϵ_CMB = 0.26 eV cm⁻³, T_dust = 20 K, ϵ_dust = 0.3 eV cm⁻³, T_starlight = 5000 K, and ϵ_starlight = 0.3 eV cm⁻³, respectively. The distance of the source was kept the same as the pulsar at 2.5 kpc. Finally, to calculate the intrinsic surface brightness profiles (SBPs) of different models, we convolved the model predictions with a symmetric 2D Gaussian point spread function (PSF; see Eq. (C.8) and the discussion in Appendix C). The adopted width σ_PSF in the figures is illustrative, serving only to demonstrate the effect of instrumental smearing, and should not be taken as corresponding to any specific instrument response. In practice, σ_PSF can be set to the appropriate value for a given telescope. For this work, as a representative case, we considered the width of the 2D Gaussian PSF as 0.35° (in the 8–40 TeV range), 0.25° (in the 40–100 TeV range), and 0.18° (in the 100–150 TeV range). The primary free parameters varied to explain the gamma-ray data for the three models are the conversion efficiency (η) and the cutoff energy (E_e, c). Furthermore, the inner zone diffusion coefficient in the 2ZISD model, and different combinations of M_A and ζ in the AD model have been explored. In the B2D model, the diffusion coefficient is fixed at the standard ISM value, i.e., D_ISM, 0 = 3.16 × 10²⁹ cm² s⁻¹ at 100 TeV, and Kolmogorov-like δ = 0.3. Figures 4–7 show the modeled SEDs and SBPs of HESS J1813-126 under the three transport models introduced in Appendix C. In the SED panels of all figures, the primary curve corresponds to an aperture of 0.34° (the upper limit of integration in Eq. (C.6)), matching the spectral extraction region measured by H.E.S.S. The additional curves show progressively larger apertures, which allow us to track how the integrated flux depends on the chosen extraction region. The considered apertures have also been indicated in the SBPs using vertical lines.

Fig. 4. Exploration of SEDs and SBPs of HESS J1813-126 in the 2ZISD framework. Left: SEDs corresponding to different choices of inner zone diffusion coefficients. The D0 in the plot signifies the normalization of the inner zone diffusion coefficient. Middle: Angle-integrated SEDs with increasing aperture sizes starting from 0.34°, assuming D0 = 5 × 1027 cm2 s−1. In these two panels, the brown squares and green triangles (along with the upper limits) correspond to the off-pulse emission from Ackermann et al. (2011) and this work (PS J1813.3-1246), respectively. The blue data points and butterfly plot correspond to the H.E.S.S. data, whereas the purple, pink, and sky blue butterfly plots correspond to the HAWC, LHAASO WCDA, and KM2A observations, respectively. We also provide the sensitivity curves for the CTA-North (gray), CTA-South (yellow-green), and SWGO (red) for comparison. Right: SBPs for three different energy ranges scaled to start from the same point. Vertical lines correspond to the aperture sizes and have the same color scheme as in the middle panel.

In Fig. 4 we show the results of the 2ZISD model. Given that the diffusion coefficient of the inner zone is unknown and cannot yet be constrained due to a lack of observed SBP data, we considered five different values of it that aptly cover the range of values typically considered in pulsar halo models. The SEDs, corresponding to these five choices, and integrated within 0.34°, reproduce the H.E.S.S., HAWC, and LHAASO data well. The SEDs also do not exceed the upper limit provided by Fermi-LAT analysis. The required values of free parameters are given in the Table 1. Then, we considered a representative inner diffusion coefficient value of 5 × 10²⁷ cm² s⁻¹, consistent with that typically considered (see, e.g., Xi et al. 2019), and enlarged the aperture in multiple steps, finding that the SED primarily increases at the lower energies. This is understandable since the lower-energy electrons, having longer cooling timescales, can propagate farther from the pulsar before losing their energy, making the IC emission more spatially extended at lower energies. As a result, when a larger aperture is considered, more of this extended, lower-energy emission from the outer regions of the halo is collected, so the integrated SED increases significantly. From Fig. 4, it can be seen that SEDs corresponding to enlarged apertures are consistent with the gamma-ray spectra posited by HAWC and LHAASO, before getting saturated. Enlarging the aperture even further does not increase the SED any more, since the steep decline of the particle density at the bubble boundary prevents much additional contribution. The corresponding SBPs, shown in the right panel of Fig. 4, reinforce this picture. At 8–40 TeV, the profile exhibits a bright central peak followed by a radial decline, indicating that the bulk of the emission is confined within the inner region. As the photon energy increases, the SBPs become progressively narrower, reflecting the shorter cooling timescales of higher-energy electrons that prevent them from diffusing to large distances. The contrast between the extended low-energy profile and the compact high-energy profiles demonstrates the energy-dependent morphology that naturally arises in the 2ZISD framework. In practice, this implies that any aperture growth primarily enhances the flux in the lower-energy bands, while the high-energy emission essentially gets saturated. From the same logic, future instruments observing the source in the TeV range, such as Cherenkov Telescope Array (CTA), are expected to observe even larger source extent and comparatively elevated gamma-ray emission from the HESS J1813-126 source region.

The results of the B2D model are shown in Fig. 5. Within the 0.34° aperture, the predicted flux reproduces the H.E.S.S. data reasonably well. The required parameters to match the spectrum are an efficiency factor of η ≈ 0.22 and a cutoff energy of E_e, c ≈ 80 TeV. However, as the aperture is increased beyond the H.E.S.S. extraction region, the integrated SED curves shift upward and begin to overshoot the HAWC and LHAASO measurements. This behavior is a direct consequence of the quasi-ballistic propagation of pairs, which produces flatter SBPs extending to larger radii. In contrast to the sharp drop-off seen in the 2ZISD case, the B2D profiles decline gradually with angle, so that a significant fraction of the flux resides outside the 0.34° aperture. As larger integration regions are considered, additional emission from the outer halo is collected. While the effect is most pronounced at lower energies – where long-lived electrons produce highly extended IC emission – the flatter SBPs also allow some higher-energy particles to contribute at larger radii, clearly enhancing the high-energy part of the SED as well. Nevertheless, these high-energy electrons are still dominated by radiative losses, so the SBPs continue to retain an energy dependence, i.e., higher-energy profiles remain more compact than those at lower energies. The predicted flux increasingly overshoots the observed HAWC and LHAASO spectra, as well as Fermi-LAT upper limits once apertures larger than 0.34° are considered. Although the required efficiency is not unphysically larger than unity (as found in some earlier works; see Wu et al. 2024), it is nevertheless high (≳10%), which is considerably above typical expectations for pulsar halos. Thus, while the B2D framework can reproduce the H.E.S.S. data, the combination of a high conversion efficiency and the systematic overprediction of the HAWC and LHAASO fluxes indicates that this model provides a poorer overall description of HESS J1813-126.

Fig. 5. Exploration of the SEDs and SBPs of HESS J1813-126 in the B2D framework. Left: Angle-integrated SEDs, calculated for increasing aperture sizes, together with the same data points, upper limits, butterfly plots, and sensitivity curves as in Fig. 4, using the same color scheme. Right: SBPs for three different energy ranges scaled to start from the same point. Vertical lines correspond to the aperture sizes and follow the same color scheme as in the left panel.

Finally, we turn to the results of the AD model. Because the predicted emission is strongly geometry-dependent, we explored the parameter space by considering M_A = 0.2, 0.24, 0.35 and 0.43, and ζ = 0° and 10°. In the AD scenario, the parallel diffusion coefficient D_zz must remain consistent with the average Galactic value, so M_A has been chosen in such a way that we can explore different perpendicular diffusion coefficients D_rr. The two values of ζ have been chosen to probe different symmetries. As before, assuming an aperture of 0.34°, the SEDs for different combinations of M_A and ζ are shown in Fig. 6, with the corresponding parameter values summarized in Table 2. From the table, one can see the efficiency factor η (as well as E_e, c) in the AD model is comparable with those in the 2ZISD model. Additionally, note that for a given M_A, increasing ζ reduces the fraction of projected emission recovered within the chosen aperture and therefore requires even higher efficiency (η) to match the observed flux. For the present work, we restricted ourselves to the representative cases of M_A = 0.35, ζ = 0°, and M_A = 0.35, ζ = 10°, and progressively enlarged the aperture (see Fig. 6). In the first scenario, the SED rises primarily at lower energies as the aperture is increased, eventually becoming consistent with the observed HAWC, LHAASO, and Fermi-LAT fluxes before getting saturated. On the second scenario, although the primary increase in the SEDs remains in the lower energy range, the model SEDs eventually overshoot the HAWC, LHAASO, and Fermi-LAT fluxes before saturation, a behavior that can be further understood by examining the model SBPs. Figure 7 shows the SBPs and gamma-ray morphologies predicted by the AD model. For ζ = 0°, where the line of sight is aligned with the ordered magnetic field, the projected gamma-ray emission map is nearly isotropic and circular, and the SBPs extend to large radii with a steep decline, yielding a single, centrally peaked profile in each energy band. In contrast, for ζ = 10°, the projected morphology becomes slightly elongated along the magnetic field: the emission extends further in the field direction but falls off more rapidly across it. This broken symmetry introduces distinct profiles, producing multiple curves in the SBPs and a more complex overall structure. From Fig. 7, it can be seen that the integrated SBP curve corresponding to the azimuthal interval around φ ∼ 90° lies below the others, whereas the same for the azimuthal intervals of φ ≪ 90° and φ ≫ 90° produce comparatively broader, flatter profiles. This happens as the viewing angle increases, the line of sight intersects a longer segment of the field-aligned diffusion region, causing the integrated SBPs in φ ≪ 90° and φ ≫ 90° cases to appear flatter even though the same for φ ∼ 90° case remains steeper. As a result, for ζ = 0° the integrated SED quickly saturates with aperture since most of the flux is already enclosed within the compact, symmetric region, while for ζ = 10° the flatter SBPs allow additional extended emission to contribute at larger radii, leading to a gradual flux increase and a mild overshoot before eventual saturation. The energy dependence further reflects the effect of radiative cooling: in the 8–40 TeV band, the emission is broad and extended, at 40–100 TeV, the profile starts to contract, and at 100–150 TeV, the emission is comparatively more suppressed at large radii, yielding a more centrally peaked profile.

Fig. 6. Exploration of SEDs of HESS J1813-126 in the AD framework. Left: Model SEDs of HESS J1813-126 for different MA − ζ combinations. Middle and right: SEDs for the combinations MA = 0.35 and ζ = 0°, and MA = 0.35 and ζ = 10°, respectively. In both panels, the primary SEDs correspond to an aperture of 0.34°, while additional SEDs show progressively larger apertures. The same data points, upper limits, butterfly plots, and sensitivity curves are provided as in Fig. 4, using the same color scheme.

Fig. 7. Exploration of SBPs and gamma-ray maps of HESS J1813-126 in the AD framework. The model SBPs are shown for two combinations of MA and ζ and for three energy ranges, scaled to start from the same point, and along with vertical lines, corresponding to the same aperture sizes and with the same color scheme as in Fig. 6. The SBP curves are divided across different azimuthal intervals, visualized as circular pie charts in the corresponding plots. The figure also shows gamma-ray maps in three energy ranges for the two considered combinations, to emphasize the anisotropy of the AD model. Top two panels: Results for MA = 0.35 and ζ = 0°. Bottom two panels: Results for MA = 0.35 and ζ = 10°.

In summary, our exploration of the three transport models underscores both shared features and key distinctions. The 2ZISD scenario provides the most consistent description of the multi-instrument spectra (H.E.S.S., HAWC, LHAASO, Fermi-LAT), reproducing the data with reasonable efficiencies and exhibiting an aperture-dependent growth confined primarily to lower energies, without overshooting the observed data. The B2D model can match the H.E.S.S. spectrum within 0.34°, but its flatter profiles lead the integrated SED to exceed the HAWC and LHAASO fluxes at larger apertures, while also requiring comparatively high conversion efficiencies. The AD model naturally captures the geometric imprint of anisotropic diffusion, predicting characteristic SBPs and morphologies testable with future imaging atmospheric Cherenkov telescopes (IACTs). The field-aligned configuration (ζ = 0°) explains the H.E.S.S. data well and shows aperture-dependent behavior consistent with multi-instrument observations, whereas the field-misaligned case (ζ = 10°) overshoots the fluxes at wider apertures, suggesting that the halo is unlikely to be strongly asymmetric and will likely display spherical morphology. Overall, the 2ZISD framework and the field-aligned AD model – though grounded in different underlying physics – both provide a coherent and self-consistent interpretation of HESS J1813-126. Future high-resolution IACT observations will be decisive in distinguishing between these scenarios and constraining particle transport in this system.

A useful way to relate the pulsar central engine PSR J1813-1246 to the surrounding pulsar halo HESS J1813-126 is to compare the high-energy pair budgets implied by the magnetospheric SC emission and by the halo models. In both cases, the model normalizations can be translated into particle injection rates, which we express as dimensionless multiplicities normalized to the Goldreich–Julian (GJ) outflow rate. Note that these multiplicities refer strictly to the particular high-energy subset of the pair population relevant for SC emission and pulsar halo formation in our modeling and they do not represent the global pair yield of the magnetosphere.

The spatial distribution of SC-emitting particles along the gap is given by Eq. (B.12). Interpreting this as a steady-state occupancy distribution within the emission region, extending from x_in to x_out with Δ ≡ x_out − x_in, the mean residence length measured from the inner boundary is

$$\begin{array}{c}\hfill x_{\mathrm{res}}\equiv \frac{{\displaystyle {\int }_{x_{\mathrm{in}}}^{x_{\mathrm{out}}}(x-x_{\mathrm{in}})\frac{\mathit{dN}}{\mathit{dx}}dx}}{{\displaystyle {\int }_{x_{\mathrm{in}}}^{x_{\mathrm{out}}}\frac{\mathit{dN}}{\mathit{dx}}dx}}=x_0-\frac{\mathrm{\Delta }}{e^{\mathrm{\Delta }/x_0}-1}.\end{array}$$(1)

This expression has the appropriate limiting behaviors: x_res → x₀ when the distribution decays steeply (x₀ ≪ Δ), and x_res → Δ/2 when the distribution is nearly uniform (x₀ ≫ Δ).

Adopting the standard assumption of ultra-relativistic streaming at v ≃ c along open field lines, the mean residence time is t_res = x_res/c. The associated SC-emitting particle throughput is therefore

$$\begin{array}{c}\hfill {\dot{N}}_{\mathrm{mag}}=\frac{N_0}{t_{\mathrm{res}}}=\frac{c{N}_0}{x_0-\frac{\mathrm{\Delta }}{e^{\mathrm{\Delta }/x_0}-1}}.\end{array}$$(2)

This quantity represents only the rate of high-energy particles participating in the SC channel and is not intended to describe the full magnetospheric outflow.

For comparison, we normalized Ṅ_mag to the canonical GJ outflow rate (see, e.g., Amato 2024),

$$\begin{array}{c}\hfill {\dot{N}}_{\mathrm{GJ}}\simeq \frac{B_0{R}_{\ast }^3{\mathrm{\Omega }}^2}{2ec},\end{array}$$(3)

with B₀ the surface polar magnetic field, R_* the stellar radius, and Ω = 2π/P. We therefore defined a high-energy magnetospheric multiplicity:

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{mag}}\equiv \frac{{\dot{N}}_{\mathrm{mag}}}{{\dot{N}}_{\mathrm{GJ}}}.\end{array}$$(4)

On the halo side, we adopted an injection power (η Ė) carried by relativistic e^± with a spectrum given by Eq. (C.3). The corresponding particle injection rate is

$$\begin{array}{c}\hfill {\dot{N}}_{\mathrm{halo}}=\frac{\eta \dot{E}}{⟨E_e⟩},\end{array}$$(5)

where the mean injected energy is

$$\begin{array}{c}\hfill ⟨E_e⟩=\frac{{\int }_{E_{e,min}}^{E_{e,max}}{E}_e{q}_E(E_e)d{E}_e}{{\int }_{E_{e,min}}^{E_{e,max}}{q}_E(E_e)d{E}_e}.\end{array}$$(6)

This rate refers only to e^± above the adopted minimum energy threshold (E_e, min∼ GeV in this case), i.e., the portion of the particles relevant for powering the pulsar halo. Normalizing to the same Ṅ_GJ gives the halo multiplicity,

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{halo}}\equiv \frac{{\dot{N}}_{\mathrm{halo}}}{{\dot{N}}_{\mathrm{GJ}}}=\frac{\eta \dot{E}}{⟨E_e⟩{\dot{N}}_{\mathrm{GJ}}}.\end{array}$$(7)

For our model parameters, we find κ_mag ≈ 5.9 × 10⁴ and κ_halo ≈ 2.2 × 10² (corresponding to the 2ZISD case with D₀ = 5 × 10²⁷ cm² s⁻¹). The estimated magnetospheric multiplicity is consistent with the typically considered range of 10³ − 10⁵ in pulsars (Timokhin & Harding 2019; Amato 2024; Spencer & Mitchell 2025). The halo multiplicity is comfortably below the very large total multiplicities inferred for classical PWNe (κ_PWN ∼ 10⁵–10⁷), which arise primarily from low-energy pairs irrelevant for pulsar halo production. Restricted to the high-energy component, the approximate values of κ_mag and κ_halo indicate that the magnetosphere does not suffer from a deficit of energetic pairs. This comparison does not, by itself, imply that the SC-emitting region supplies the halo directly; alternative high-energy escape channels such as the equatorial wind or the current sheet may also contribute. Rather, the comparison only shows that the global high-energy pair budget in our model is energetically self-consistent without a need for additional pair-creation channels or exotic accelerators outside the magnetosphere.

3. Concluding remarks

A compelling case can be made that HESS J1813-126 is indeed a pulsar halo powered by PSR J1813-1246. The multiwavelength phenomenology strongly favors a halo interpretation: no radio or X-ray PWN has been detected despite deep observations (Ackermann et al. 2011; Marelli et al. 2014; Guevel et al. 2025), which disfavors a confined nebula scenario. The TeV emission is extended (∼0.2°–0.3°), with a hard spectrum, and with no shell-like morphology to suggest a SNR or associated MCs (Marelli et al. 2014). The GeV counterpart detected by Fermi-LAT is point-like and coincident with the pulsar, consistent with being dominated by magnetospheric emission. Our phased analysis shows that any residual off-pulse GeV emission is soft and compact, and is more plausibly explained by the imperfect gating of magnetospheric photons rather than an underlying GeV PWN. The extended TeV source HESS J1813-126 has been independently detected by H.E.S.S., HAWC, and LHAASO, all of which confirm its extended morphology and reveal emission up to > 100 TeV, establishing it as a bright Galactic PeVatron candidate. The positional coincidence with PSR J1813-1246, an energetic, middle-aged, radio-quiet pulsar located off the Galactic plane, makes a chance alignment highly unlikely. Our transport modeling shows that the observed H.E.S.S. SED can be consistently explained within the pulsar halo framework under multiple plausible transport scenarios (2ZISD, B2D, and AD models). Each yields distinct predictions for the SBPs and the aperture-dependent evolution of the SED that provide clear observational diagnostics for future IACTs, yet all converge on the idea that PSR J1813-1246 has efficiently injected e^± pairs into the ISM over its lifetime. Taken together, the absence of a classical PWN, the extended and multi-TeV nature of the emission, the spatial coincidence with a powerful pulsar, and the successful reproduction of the data by halo models provide strong evidence that HESS J1813-126 is best interpreted as a pulsar halo. Recent reanalyses of several extended TeV sources by the H.E.S.S. collaboration have shown that some objects previously reported as moderately extended in the HGPS can exhibit substantially larger and more complex morphologies when studied with improved reconstruction techniques and deeper exposures. Notable examples include HESS J1809-193 and HESS J1813-178, for which updated analyses (H.E.S.S. Collaboration 2023, 2024) revealed broad, multi-degree emission structures. These results have transformed our understanding of evolved PWNe, highlighting scenarios in which a significant fraction of pairs escape the nebula and form very extended pulsar halos, as also discussed for HESS J1809-193 (H.E.S.S. Collaboration 2023; Martin et al. 2024). In this broader context, a dedicated, high-statistics H.E.S.S. reanalysis of HESS J1813-126 would be particularly valuable, as it may uncover similarly extended emission components and provide decisive constraints on the transport physics operating in this system.

In addition to the halo modeling, we also modeled the high-energy pulsar emission with the SC framework, which successfully reproduces the combined X-ray and gamma-ray spectra of PSR J1813-1246. The approximate high-energy pair multiplicities found for the magnetospheric and halo populations in our model suggest that the global energetic requirements for both the SC emission and the pulsar halo emission can be met without invoking additional pair-creation channels. This reinforces the association between the pulsar and its surrounding pulsar halo, as both the magnetospheric and halo phenomenology can be explained within a coherent and unified picture.

Acknowledgments

A.D.S. thanks the anonymous reviewer for helpful suggestions and constructive criticism, which vastly improved the quality of this work. This work is part of the grant Juan de la Cierva JDC2023-052168-I, funded by MCIU/AEI/10.13039/501100011033 and by the ESF+. This work has been supported by the grant PID2024-155316NB-I00 funded by MICIU/AEI/10.13039/501100011033/FEDER, UE and CSIC PIE 202350E189. This work was also supported by the Spanish program Unidad de Excelencia María de Maeztu CEX2020-001058-M and also supported by MCIN with funding from European Union NextGeneration EU (PRTR-C17.I1).

References

Appendix A: Fermi-LAT data analysis

To investigate the gamma-ray emission from PSR J1813-1246, we analyzed ∼16.75 years of P8R3 SOURCE class Fermi-LAT data (evclass = 128, evtype = 3; Atwood et al. 2013; Bruel et al. 2018), covering 2008 August 4 (Mission Elapsed Time (MET) 239557417) to 2025 May 2 (MET 767880664) in the 1–500 GeV energy range. We adopted a minimum energy threshold of 1 GeV to reduce the impact of the bright diffuse emission and to exploit the substantially improved Fermi-LAT PSF at higher energies. The analysis was performed with the Fermipy² package (v1.4.0; Wood et al. 2017), which interfaces with the official Fermitools³. We adopted the instrument response function P8R3_SOURCE_V3, excluded events with zenith angle > 105° to minimize Earth-limb contamination, and used the standard Galactic (galdiff) and isotropic (isodiff) diffuse models, gll_iem_v07.fits and iso_P8R3_SOURCE_V3_v1.txt, respectively. For the source modeling, we adopted the latest Fermi-LAT source catalog, 4FGL-DR4 (Abdollahi et al. 2020; Ballet et al. 2023).

The region of interest was defined as a circle with 20° radius centered on PSR J1813-1246, with a spatial bin size of 0.1° and 8 energy bins per decade. A binned maximum-likelihood analysis was carried out over a 15° ×15° region centered on the pulsar. We first optimized the baseline model by adjusting the normalizations of diffuse components and 4FGL sources within the analysis region to establish a stable starting point for subsequent likelihood fitting. Then, we freed the normalizations of the 4FGL cataloged sources within 5° of the pulsar, along with the diffuse components, while fixing distant sources to their 4FGL-DR4 values. To search for un-modeled excesses, we applied Fermipy’s source-finding algorithm, requiring a test statistic (TS) > 16 and separation of at least 0.2° from known sources, followed by a global refit by repeating the previous step. From this procedure, the spectral and spatial parameters of all relevant sources were evaluated, and subsequently, the SED, extension, and localization of the target were computed. In Fig. A.1, we show the smoothed data counts map, the residual counts map (with all sources marked), and the $\sqrt{\mathit{TS}}$ significance map of the 10° ×10° region around PSR J1813-1246 to demonstrate that all signal has been properly modeled.

Fig. A.1. Fermi-LAT analysis diagnostic plots. Smoothed data counts map (left panel) and residual (data−model) counts map (middle panel) of a 10° ×10° region centered on PSR J1813-1246 are shown. The residual map includes the positions of all the 4FGL sources in the region (black plus signs). The color bars in both panels indicate the photon counts. The right panel shows the $\sqrt{\mathit{TS}}$ significance map of the same region, where the color bar indicates the $\sqrt{\mathit{TS}}$ value of the region. The position of PSR J1813-1246 is marked with a green cross in all panels.

To specifically probe possible GeV PWN emission hidden beneath the pulsar contribution, we performed a pulsar phased analysis using ∼11 years of Fermi-LAT data, from 2008 August 11 (MET 240162626) to 2019 June 17 (MET 582479577), in the energy range 1–500 GeV. We considered ∼ 11 years of data in this analysis, corresponding to the time span for which a valid ephemeris was available. Event files with pulse phases were obtained from the Fermi-LAT third pulsar catalog (3PC)⁴, with spacecraft files from the Fermi-LAT data server⁵. The analysis was performed in a 15° ×15° region centered on the pulsar position. We followed the same event selection and diffuse models as above. The phased analysis proceeded in three main steps. First, we optimized the baseline model as before. Then, the initial fit was performed by freeing the normalizations of all 4FGL sources within 5° of the pulsar having a TS > 9, together with the normalizations of the diffuse templates galdiff and isodiff. Second, we applied the Fermipy source-finding algorithm to search for additional excesses in the analysis region, requiring candidate sources to exceed TS = 16 and to be separated by at least 0.2° from cataloged sources to avoid confusion. Third, we refitted the updated source model, performing a global maximum-likelihood optimization to obtain stable convergence. During the off-pulse interval, this procedure revealed a new source at the pulsar location, which was then added to the model. Its SED was extracted using the standard Fermipy SED tool, and we tested for possible spatial extension with RadialDisk and RadialGaussian morphologies. We used likelihood ratio tests to refine its localization and to compare point-like and extended-source hypotheses. Note that throughout this work, we have assumed that the extended hypothesis is true only if TS_ext ≥ 25. The properties of this new source were subsequently examined as potential evidence of faint GeV emission associated with an underlying PWN, emerging once the dominant pulsed magnetospheric emission was gated off.

Appendix B: The synchro-curvature model

For a particle of Lorentz factor Γ and pitch angle α moving along a field line of curvature radius r_c in a local magnetic field B, the gyroradius and SC parameter are

$$\begin{array}{cc}\hfill r_{\mathrm{gyr}}& =\frac{m_e{c}^2\mathrm{\Gamma }sin\alpha }{\mathit{eB}},\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill \xi & =\frac{r_{\mathrm{c}}}{r_{\mathrm{gyr}}}\frac{{sin}^2\alpha }{{cos}^2\alpha }.\hfill \end{array}$$(B.2)

The corresponding characteristic photon energy is

$$\begin{array}{c}\hfill E_{\mathrm{c}}(\mathrm{\Gamma },r_{\mathrm{c}},r_{\mathrm{gyr}},\alpha )=\frac{3}{2}ħc{Q}_2{\mathrm{\Gamma }}^3,\end{array}$$(B.3)

with

$$\begin{array}{c}\hfill Q_2^2=\frac{{cos}^4\alpha }{r_{\mathrm{c}}^2}(1+3\xi +{\xi }^2+\frac{r_{\mathrm{gyr}}}{r_{\mathrm{c}}}).\end{array}$$(B.4)

The single-particle spectral power per unit energy at a given position then reads

$$\begin{array}{c}\hfill \frac{d{P}_{\mathrm{sc}}}{\mathit{dE}}=\frac{\sqrt{3}{e}^2\mathrm{\Gamma }y}{4\pi ħr_{\mathrm{eff}}}[(1+z)F(y)-(1-z)K_{2/3}(y)],\end{array}$$(B.5)

where

$$\begin{array}{c}\hfill y\equiv \frac{E}{E_{\mathrm{c}}},z={(Q_2{r}_{\mathrm{eff}})}^{-2},r_{\mathrm{eff}}=\frac{r_{\mathrm{c}}}{{cos}^2\alpha }{(1+\xi +\frac{r_{\mathrm{gyr}}}{r_{\mathrm{c}}})}^{-1},\end{array}$$(B.6)

and

$$\begin{array}{c}\hfill F(y)={\displaystyle {\int }_y^{\infty }}{K}_{5/3}(y^{\prime })d{y}^{\prime },\end{array}$$(B.7)

with K_n being the modified Bessel functions of the second kind of index n. The total SC power radiated by a single particle is

$$\begin{array}{c}\hfill P_{\mathrm{sc}}=\frac{2e^{2}c{\mathrm{\Gamma }}^4}{3r_{\mathrm{c}}^2}{g}_r,\mathrm{where}{g}_r=\frac{r_{\mathrm{c}}^2}{r_{\mathrm{eff}}^2}\frac{1+7{(Q_2{r}_{\mathrm{eff}})}^{-2}}{8{(Q_2{r}_{\mathrm{eff}})}^{-1}}.\end{array}$$(B.8)

The evolution of Γ, α, and ξ is obtained by solving the coupled differential equations

$$\begin{array}{cc}\hfill \frac{d(psin\alpha )}{\mathit{dt}}& =-\frac{P_{\mathrm{sc}}}{c}sin\alpha ,\hfill \end{array}$$(B.9)

$$\begin{array}{cc}\hfill \frac{d(pcos\alpha )}{\mathit{dt}}& =e{E}_{\Vert }-\frac{P_{\mathrm{sc}}}{c}cos\alpha ,\hfill \end{array}$$(B.10)

with p(=Γm_ec) being the relativistic momentum and x (dx = c dt) being the gap coordinate. We adopted Γ_in = 10³ and α_in = 45° as initial conditions. The average SC spectrum throughout the trajectory is then

$$\begin{array}{c}\hfill \frac{d{P}_{\mathrm{tot}}}{\mathit{dE}}={\displaystyle {\int }_{x_{\mathrm{in}}}^{x_{\mathrm{out}}}}\frac{d{P}_{\mathrm{sc}}}{\mathit{dE}}(x)\frac{\mathit{dN}}{\mathit{dx}}(x)dx,\end{array}$$(B.11)

with an effective particle distribution

$$\begin{array}{c}\hfill \frac{\mathit{dN}}{\mathit{dx}}(x)=\frac{N_0{e}^{-(x-x_{\mathrm{in}})/x_0}}{x_0(1-e^{-(x_{\mathrm{out}}-x_{\mathrm{in}})/x_0})}.\end{array}$$(B.12)

In our setup, the inner boundary (x_in) is at 0.5 R_lc and the outer boundary (x_out) is at 1.5 R_lc of the emission region, and x₀ is the length scale, which was considered a free parameter. The magnetic field and curvature radius along the path are parameterized as

$$\begin{array}{c}\hfill B(x)=B_0{\left(\frac{R_{\ast }}{x}\right)}^{b_{\mathit{sc}}},r_{\mathrm{c}}(x)=R_{\mathrm{lc}}{\left(\frac{x}{R_{\mathrm{lc}}}\right)}^{{\eta }_{\mathit{sc}}},\end{array}$$(B.13)

where R_* (≈ 10 km) is the pulsar radius, $B_0(\approx 6.4\times {10}^{19}\sqrt{P\dot{P}}\mathrm{G})$ is the dipolar magnetic field at the polar surface, R_lc (=Pc/2π) is the light cylinder radius and η_sc = 0.5. The magnetic field gradient b_sc, which determines how quickly the magnetic field decreases along particle trajectories, along with the $N_0(={\displaystyle {\int }_{x_{\min }}^{x_{\max }}}\frac{\mathit{dN}}{\mathit{dx}}dx)$, denoting the total number of particles and used to set the normalization of the spectrum, were also kept free.

Appendix C: The pulsar halo particle transport models

The evolution of the lepton distribution N(E_e, r, t) is governed by the time-dependent, diffusion-loss equation

$$\begin{array}{c}\hfill \frac{\partial N}{\partial t}=\mathrm{\nabla }·(D(E_e,r)\mathrm{\nabla }N)+\frac{\partial }{\partial E_e}(b(E_e)N)+Q(E_e,r,t),\end{array}$$(C.1)

where D is the diffusion coefficient, b(E_e) = |dE_e/dt| is the absolute energy-loss rate, and Q is the source term. For relativistic leptons above GeV energies, the losses are quadratic in energy, b(E_e)≈b₀E_e², with b₀ determined by the sum of synchrotron and IC contributions. The IC term is evaluated for each photon field with the full Klein–Nishina cross section.

The source term Q can be divided into spatial (q_r), temporal (q_t), and energy components (q_E). The spatial component of injection is taken to be point-like in space at the pulsar position, i.e., q_r(r) = δ(r). The time dependence follows the pulsar spin-down power in the dipole case (n = 3):

$$\begin{array}{c}\hfill q_t(t)=\frac{{(1+t/{\tau }_0)}^{-2}}{{(1+t_s/{\tau }_0)}^{-2}},t\ge 0,\end{array}$$(C.2)

where τ₀ is the initial spin-down timescale and t_s the source age. The spectrum of injected particles at the current epoch is assumed as

$$\begin{array}{c}\hfill q_E(E_e)=q_{E,0}{E}_e^{-p}exp[-{\left(\frac{E_e}{E_{e,c}}\right)}^s],\end{array}$$(C.3)

within E_e, min ≈ 1 GeV and E_e, max ≈ 1 PeV with normalization set by an efficiency factor η through η Ė = ∫q_E(E_e) E_e dE_e.

Once N(E_e, r) is obtained by solving Eq. C.1, the observable gamma-ray signal is calculated. For the spherically symmetric case, the surface distribution is calculated by performing line-of-sight (LOS) integration on the electron number density:

$$\begin{array}{c}\hfill S_e(E_e,\theta )={\displaystyle {\int }_0^{\infty }N(E_e,r)dl,}\end{array}$$(C.4)

with $r=\sqrt{r_s^2+l^2-2r_{s}lcos\theta }$, where r_s is the distance to the pulsar, l is the distance from a point to the observer along LOS, and θ is the angle between LOS and the direction from the observer to the pulsar. The IC surface brightness is then

$$\begin{array}{c}\hfill s_{\gamma }(E_{\gamma },\theta )=\frac{1}{4\pi }{\displaystyle \int d}{E}_e{\displaystyle \int d}ϵS_e(E_e,\theta )f(E_e,ϵ,E_{\gamma }),\end{array}$$(C.5)

where f signifies the production rate of emitted gamma-ray photons with energy E_γ from the interaction between electrons of energy E_e and target photons of energy ϵ. The energy spectrum within an angular range can be calculated using

$$\begin{array}{c}\hfill F(E_{\gamma })={\displaystyle {\int }_{{\theta }_{\gamma ,1}}^{{\theta }_{\gamma ,2}}}{s}_{\gamma }(E_{\gamma },\theta )2\pi \theta d\theta ,\end{array}$$(C.6)

whereas the energy band-integrated SBP is

$$\begin{array}{c}\hfill S_{\gamma }(\theta )={\displaystyle {\int }_{E_{\gamma ,1}}^{E_{\gamma ,2}}}{s}_{\gamma }(E_{\gamma },\theta )w(E_{\gamma })d{E}_{\gamma },\end{array}$$(C.7)

with w(E_γ) = 1 for photon counts or w(E_γ) = E_γ for energy flux. Finally, this intrinsic profile must be convolved with the instrument’s PSF:

$$\begin{array}{c}\hfill {\stackrel{\sim }{S}}_{\gamma }(\theta )={\displaystyle {\int }_0^{\infty }{\int }_0^{2\pi }}{S}_{\gamma }({\theta }^{\prime })g(\theta ,{\theta }^{\prime },\varphi ){\theta }^{\prime }d\varphi d{\theta }^{\prime },\end{array}$$(C.8)

where (θ, θ′,ϕ) denote polar coordinates in the plane of the sky, with θ and θ′ being the observed and intrinsic angular distances from the pulsar, respectively, and ϕ being the azimuthal angle. g is the PSF. In this work, a symmetric 2D Gaussian PSF with width σ_PSF is considered. For a cylindrical model, such as the AD model, the LOS integration takes the form

$$\begin{array}{c}\hfill S_e(E_e,\varphi ,\theta )={\displaystyle {\int }_0^{\infty }}N(E_e,r,z)dl,\end{array}$$(C.9)

where the mapping of (l, θ, ϕ)→(r, z) has been done following Liu et al. (2019). We briefly discuss the underlying idea and corresponding equations of the three models below.

C.1. Two-zone isotropic suppressed diffusion model

In isotropic cases, the diffusion-loss equation reduces to spherical symmetry:

$$\begin{array}{c}\hfill \frac{\partial N(E_e,r,t)}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^{2}D(E_e,r)\frac{\partial N}{\partial r}\right)+\frac{\partial }{\partial E_e}(b(E_e)N)+Q(E_e,t)\delta (r).\end{array}$$(C.10)

A single-zone suppressed diffusion model is generally inconsistent with the B/C ratio and cosmic-ray electron data (Aguilar et al. 2016; H.E.S.S. Collaboration 2017; Hooper & Linden 2018), motivating the two-zone approach (Hooper et al. 2017; Fang et al. 2018; Profumo et al. 2018; Tang & Piran 2019; Jóhannesson et al. 2019). This introduces a slow-diffusion bubble of radius r_2z (tens of parsecs) around the pulsar,

$$\begin{array}{c}\hfill D(E_e,r)=\{\begin{array}{cc}{D}_{\mathrm{in}}(E_e)={\mathrm{\Xi }}_{2\mathrm{z}}{D}_{\mathrm{ISM},0}{\left(\frac{E_e}{100\mathrm{TeV}}\right)}^{\delta },\hfill & r\le r_{2\mathrm{z}},\hfill \\ D_{\mathrm{out}}(E_e)=D_{\mathrm{ISM}}=D_{\mathrm{ISM},0}{\left(\frac{E_e}{100\mathrm{TeV}}\right)}^{\delta },\hfill & r>r_{2\mathrm{z}},\hfill \end{array}\end{array}$$(C.11)

with Ξ_2z ≪ 1. In this work, we adopted r_2z = 50 pc (Xi et al. 2019; Wu et al. 2024). The diffusion coefficient D_ISM, 0 = 3.16 × 10²⁹ cm² s⁻¹ at 100 TeV, and Kolmogorov-like δ = 0.3 were used unless otherwise stated. This configuration reproduces the compact halos observed by HAWC around Geminga and Monogem (Abeysekara et al. 2017; Hooper et al. 2017). Unlike the B2D or AD models, it enforces a sharp radial break in D(E_e, r) with suppressed diffusion within r_2z. This suppressed diffusion ensures particles cool before escaping, yielding a bright inner core and a drop in brightness beyond r_2z. Cosmic-ray self-generated turbulence, via resonant and nonresonant streaming instabilities, has been proposed as the origin of this slow-diffusion region (Skilling 1971, 1975; Bykov et al. 2013; Evoli et al. 2018a,b; Mukhopadhyay & Linden 2022; Fang et al. 2019; Bell 2004; Gupta et al. 2021; Schroer et al. 2021). Although semi-analytical solutions are present for the 2ZISD model (Osipov et al. 2020), we used the finite-volume discretization method (Fang et al. 2018) to solve Eq. C.10 with discontinuous diffusion coefficients in Eq. C.11.

C.2. Ballistic-to-diffusion model

The one-zone normal diffusion model with spherical transport (Eq. C.10) suffers from an unphysical superluminal problem. With mean free path λ(E_e) = 3D(E_e)/c, freshly injected particles cannot diffuse for t ≤ λ/c; otherwise, their average displacement, $\sqrt{4Dt}$, would exceed that at light speed (Aloisio et al. 2009). A more physical picture is ballistic propagation at r ≪ λ and diffusion at r ≫ λ, implying a transition between regimes.

This can be treated by introducing the Jüttner propagator (Aloisio et al. 2009), which preserves causality and interpolates between ballistic and diffusive transport. The diffusion-loss solution then reads

$$\begin{array}{cc}\hfill N(E_e,r)& ={\displaystyle \int d}{t}_0\frac{q_t(t_0)q_E(E_{e,★})}{b(E_e)}\frac{b(E_{e,★})}{4\pi {[c(t-t_0)]}^3}\hfill \\ \hfill & \times \frac{H(1-{\xi }_1){\kappa }_1}{{(1-{\xi }_1^2)}^2{K}_1({\kappa }_1)}exp[-\frac{{\kappa }_1}{\sqrt{1-{\xi }_1^2}}],\hfill \end{array}$$(C.12)

where ξ₁ = r/[c(t − t₀)], κ₁ = [c(t − t₀)]²/(2λ_★), λ_★ = ∫_Ee^E_e, ★D(E′_e) dE′_e/b(E′_e), and K₁ is the first-order modified Bessel function of the second kind.

During quasi-ballistic propagation, the velocity distribution is anisotropic, so IC emission assuming isotropy (Eq. C.4) is invalid. A correction is included in the LOS integration:

$$\begin{array}{c}\hfill S_e(E_e,\theta )={\displaystyle {\int }_0^{\infty }}2N(E_e,r)M(x,\mu )dl,\end{array}$$(C.13)

with (Aharonian et al. 2010; Prosekin et al. 2015)

$$\begin{array}{c}\hfill M(x,\mu )=\frac{1}{Z(x)}exp[-\frac{3(1-\mu )}{x}],\end{array}$$(C.14)

where x = rc/D(E_e), $Z(x)=\frac{x}{3}[1-exp(-6/x)]$, and μ = (r_s cos θ − l)/r. The IC emission was then evaluated with Eq. C.5.

Recchia et al. (2021) demonstrated that the B2D framework can reproduce the multi-TeV SBPs of pulsar halos using a standard ISM diffusion coefficient, albeit at the cost of invoking larger pair conversion efficiencies compared to the 2ZISD model. In this scenario, the halo can appear more extended, with a possible shrinkage at high energies due to the limited ballistic propagation length and radiative cooling effects (Prosekin et al. 2015; Yang & Liu 2022). In the B2D model, the apparent halo size is mainly set by the ballistic propagation length of pairs before they isotropize and enter the diffusive regime, whereas in the 2ZISD model it is controlled by the combined effects of radiative cooling and the diffusion distance of the pairs.

C.3. Anisotropic diffusion model

Given that interstellar turbulence is sub-Alfvénic (Alfvénic Mach number M_A < 1), and the magnetic field maintains a mean direction within one coherence length of 10 - 100 pc, transport is no longer isotropic within that length, since particle diffusion along the mean magnetic field direction is much faster than that in the perpendicular direction. In this scenario, the perpendicular diffusion coefficient D_rr is suppressed by a factor M_A⁴ with respect to the parallel diffusion coefficient D_zz (Yan & Lazarian 2008), assuming the mean magnetic field is oriented along the z axis. For this model, the diffusion–loss equation then takes the cylindrical form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial N(E_e,r,z,t)}{\partial t}& =\frac{1}{r}\frac{\partial }{\partial r}\left(r{D}_{\mathit{rr}}(E_e)\frac{\partial N}{\partial r}\right)+D_{\mathit{zz}}(E_e)\frac{{\partial }^{2}N}{\partial z^2}\hfill \\ \hfill & +\frac{\partial }{\partial E_e}(b(E_e)N)+Q(E_e,t)\delta (r)\delta (z),\hfill \end{array}\end{array}$$(C.15)

with

$$\begin{array}{c}\hfill D_{\mathit{rr}}(E_e)=D_{rr,0}{\left(\frac{E_e}{100\mathrm{TeV}}\right)}^{\delta },D_{\mathit{zz}}(E_e)=M_A^{-4}{D}_{\mathit{rr}}(E_e).\end{array}$$(C.16)

After rescaling z′=M_A²z, the diffusion becomes isotropic in (r, z′) coordinates, giving

$$\begin{array}{c}\hfill N(E_e,r,z)={\displaystyle \int d}{t}_0\frac{q_t(t_0)q_E(E_{e,★})}{b(E_e)}\frac{b(E_{e,★})M_A^2}{{(4\pi {\lambda }_{★})}^{3/2}}exp[-\frac{r^2+M_A^4{z}^2}{4{\lambda }_{★}}],\end{array}$$(C.17)

where λ_★ = ∫_Ee^E_e, ★D_rr(E′_e) dE′_e/b(E′_e). The AD model predicts a distinctive observational signature of elongated, asymmetric halo morphology that depends on the viewing angle ζ relative to the ordered field. Liu et al. (2019) demonstrated that the multi-TeV SBP of Geminga can be reproduced if ζ ≤ 5° with respect to the mean magnetic field. The AD model removes the need for suppression via turbulence amplification as present in the 2ZISD model. AD and 2ZISD models typically yield comparable conversion efficiencies, as both produce steep SBP declining at small radii due to reduced particle diffusion, though the physical origin of slow diffusion differs. Nevertheless, it is to be noted that the consistency of the AD model in light of the general pulsar halo observations is still debatable (De La Torre Luque et al. 2022).

All Tables

All Figures

Fig. 1. Significance maps of the HESS J1813-126 region. Left: Significance map of 4FGL J1813.4-1246 from the phase-averaged analysis. Middle: Same but for PS J1813.3-1246, detected during the off-pulse phase of the 4FGL source. Right: H.E.S.S. significance map. The color bar in each plot provides the significance level. The source extent for H.E.S.S., HAWC, and LHAASO (KM2A and WCDA) is shown as colored circles.

In the text

Fig. 2. Weighted pulse profile of PSR J1813-1246 (4FGL J1813.4-1246) over one rotational cycle (pulse phase 0-1). It was constructed using 100 uniform phase bins per period from events between 0.1 and 500 GeV with the model weight from Bruel (2019). The blue histogram shows the weighted photon counts per phase bin, while the solid red line indicates the Bayesian block decomposition of statistically significant structures in the weighted pulse profile. The red shaded region signifies the off-pulse interval estimated from this work, whereas the magenta and green shaded regions signify the same obtained in Abdo et al. (2013) and Ackermann et al. (2011), respectively.

In the text

Fig. 3. Distance variation of different SC model parameters and the spectrum of PSR J1813-1246: Lorentz factor (Γ; first panel), SC parameter (ξ; second panel), and pitch angle (sin α; third panel). In the second panel, the dashed red line marks ξ = 1, i.e., the transition point between synchrotron- and curvature-dominated regimes. In the last panel, the model spectrum is plotted against the observed SED obtained from the Fermi-LAT data analysis of 4FGL J1813.4-1246 and the X-ray SED and upper limit from XMM-Newton and INTEGRAL data, respectively, as obtained from Marelli et al. (2014) and Guevel et al. (2025). In the last panel, decomposed emission from different parts of the emitting region is also shown; the colorbar signifies the distance from the inner boundary (xin) in units of the light-cylinder radius (Rlc).

In the text

Fig. 4. Exploration of SEDs and SBPs of HESS J1813-126 in the 2ZISD framework. Left: SEDs corresponding to different choices of inner zone diffusion coefficients. The D0 in the plot signifies the normalization of the inner zone diffusion coefficient. Middle: Angle-integrated SEDs with increasing aperture sizes starting from 0.34°, assuming D0 = 5 × 1027 cm2 s−1. In these two panels, the brown squares and green triangles (along with the upper limits) correspond to the off-pulse emission from Ackermann et al. (2011) and this work (PS J1813.3-1246), respectively. The blue data points and butterfly plot correspond to the H.E.S.S. data, whereas the purple, pink, and sky blue butterfly plots correspond to the HAWC, LHAASO WCDA, and KM2A observations, respectively. We also provide the sensitivity curves for the CTA-North (gray), CTA-South (yellow-green), and SWGO (red) for comparison. Right: SBPs for three different energy ranges scaled to start from the same point. Vertical lines correspond to the aperture sizes and have the same color scheme as in the middle panel.

In the text

Fig. 5. Exploration of the SEDs and SBPs of HESS J1813-126 in the B2D framework. Left: Angle-integrated SEDs, calculated for increasing aperture sizes, together with the same data points, upper limits, butterfly plots, and sensitivity curves as in Fig. 4, using the same color scheme. Right: SBPs for three different energy ranges scaled to start from the same point. Vertical lines correspond to the aperture sizes and follow the same color scheme as in the left panel.

In the text

Fig. 6. Exploration of SEDs of HESS J1813-126 in the AD framework. Left: Model SEDs of HESS J1813-126 for different MA − ζ combinations. Middle and right: SEDs for the combinations MA = 0.35 and ζ = 0°, and MA = 0.35 and ζ = 10°, respectively. In both panels, the primary SEDs correspond to an aperture of 0.34°, while additional SEDs show progressively larger apertures. The same data points, upper limits, butterfly plots, and sensitivity curves are provided as in Fig. 4, using the same color scheme.

In the text

Fig. 7. Exploration of SBPs and gamma-ray maps of HESS J1813-126 in the AD framework. The model SBPs are shown for two combinations of MA and ζ and for three energy ranges, scaled to start from the same point, and along with vertical lines, corresponding to the same aperture sizes and with the same color scheme as in Fig. 6. The SBP curves are divided across different azimuthal intervals, visualized as circular pie charts in the corresponding plots. The figure also shows gamma-ray maps in three energy ranges for the two considered combinations, to emphasize the anisotropy of the AD model. Top two panels: Results for MA = 0.35 and ζ = 0°. Bottom two panels: Results for MA = 0.35 and ζ = 10°.

In the text

Fig. A.1. Fermi-LAT analysis diagnostic plots. Smoothed data counts map (left panel) and residual (data−model) counts map (middle panel) of a 10° ×10° region centered on PSR J1813-1246 are shown. The residual map includes the positions of all the 4FGL sources in the region (black plus signs). The color bars in both panels indicate the photon counts. The right panel shows the $\sqrt{\mathit{TS}}$ significance map of the same region, where the color bar indicates the $\sqrt{\mathit{TS}}$ value of the region. The position of PSR J1813-1246 is marked with a green cross in all panels.

In the text