PHECT: A lightweight computation tool for pulsar halo emission

$γ$-ray pulsar halos, most likely formed by inverse Compton scattering of electrons and positrons propagating in the pulsar-surrounding interstellar medium with background photons, serve as an ideal probe for Galactic cosmic-ray propagation on small scales (typically tens of parsecs). While the associated electron and positron propagation is often modeled using homogeneous and isotropic diffusion, termed here as standard diffusion, the actual transport process is expected to be more complex. This work introduces the Pulsar Halo Emission Computation Tool (PHECT), a lightweight software designed for modeling pulsar halo emission. PHECT incorporates multiple transport models extending beyond standard diffusion, accounting for different possible origins of pulsar halos. Users can conduct necessary computations simply by configuring a YAML file without manual code edits. Furthermore, the tool adopts finite-volume discretizations that remain stable on non-uniform grids and in the presence of discontinuous diffusion coefficients. PHECT is ready for the increasingly precise observational data and the rapidly growing sample of pulsar halos.

💡 Research Summary

This paper presents PHECT (Pulsar Halo Emission Computation Tool), a lightweight, open‑source software package designed to model the γ‑ray emission from pulsar halos. Pulsar halos are extended γ‑ray structures produced when high‑energy electrons and positrons, escaped from a pulsar wind nebula (PWN), inverse‑Compton scatter ambient photon fields (cosmic microwave background, infrared dust emission, and starlight) in the surrounding interstellar medium (ISM). Because the γ‑ray surface brightness directly traces the spatial distribution of the parent leptons, halos provide a unique probe of cosmic‑ray (CR) propagation on scales of tens of parsecs.

Scientific Context
Historically, halo modeling has relied on a homogeneous, isotropic diffusion‑loss equation (the “standard diffusion” model) with a diffusion coefficient D(E)∝E^δ and energy‑loss rate b(E)∝E² (synchrotron + inverse‑Compton). Recent high‑resolution observations, however, reveal that the apparent slow diffusion (D ≈ 10⁻³–10⁻⁴ kpc² Myr⁻¹, two orders of magnitude below the Galactic average) cannot be explained by a single uniform diffusion coefficient. The authors discuss several physically motivated alternatives:

1. SNR‑induced diffusion – Turbulence generated by the host supernova remnant (SNR) or the progenitor’s stellar wind reduces the diffusion coefficient inside the SNR (D_in) relative to the ambient ISM (D_out). The model treats D as a step function with a discontinuity at the SNR boundary.

2. Anisotropic diffusion – If the mean magnetic field near the pulsar is nearly aligned with the line of sight, the diffusion coefficient parallel to the field (D_∥) can remain Galactic‑average while the perpendicular component (D_⊥) is strongly suppressed. Projection effects then mimic a slow‑diffusion halo.

3. Superluminal correction – For very fast diffusion (large D) and short propagation times, the classical diffusion equation predicts superluminal particle fronts. PHECT implements a correction term that limits the effective propagation speed to c, which is important for accurate modeling of bright, compact γ‑ray cores.

4. Superdiffusion and magnetic‑mirror diffusion – Early escape phases may be governed by Lévy‑flight‑like superdiffusion, while magnetic mirrors in the ISM can produce a “mirror diffusion” regime. PHECT currently includes a superdiffusive transport module and plans to add a time‑dependent transition to standard diffusion.

Software Architecture
PHECT is written in C with only two external dependencies: the GNU Scientific Library (GSL) for numerical integration and libyaml for configuration handling. All model parameters—pulsar spin‑down power, age, initial spin‑down timescale τ₀, injection spectral index p, cutoff energy E_c, conversion efficiency η, magnetic field strength B, background photon temperatures and energy densities, diffusion coefficients, geometry flags, etc.—are specified in a human‑readable YAML file. No source‑code modification is required to switch between models.

The computation proceeds in three stages:

1. Electron transport – The diffusion‑loss equation ∂N/∂t = ∇·(D∇N) + ∂(bN)/∂E + Q is solved using a finite‑volume (FV) discretization. The FV scheme conserves particle number across cell faces and remains stable on non‑uniform radial grids and across discontinuous D. Time integration employs a semi‑implicit Crank‑Nicolson method, allowing large time steps without sacrificing stability.

2. Line‑of‑sight integration – The three‑dimensional electron density N(E, r, t) is projected onto the sky to obtain the surface density Σ(E, θ). For cylindrically symmetric (anisotropic) models, an additional angular dimension is retained.

3. Inverse‑Compton emission – The γ‑ray surface brightness I_γ(θ, E_γ) is calculated by convolving Σ with the full Klein‑Nishina cross‑section. The background photon fields are modeled as a blackbody (CMB) plus gray‑body spectra for infrared dust and starlight. Energy‑loss rates b_syn and b_ICS are expressed as b = (b₀,syn + b₀,ics) E², with b₀,ics containing a polynomial approximation Y(x) of the Klein‑Nishina correction.

All integrals over electron energy are performed with GSL’s adaptive CQUAD routine, guaranteeing high accuracy even when the integration limits span several decades (E_min ≈ 1 GeV to E_max ≈ 50 E_c).

Validation and Applications
The authors benchmark PHECT against analytic approximations (e.g., the 1/r² surface‑brightness law) and find agreement within 1 %. They also demonstrate numerical convergence when refining the radial grid or reducing the time step, even in the presence of sharp D jumps. As a proof‑of‑concept, they model the Geminga and Monogem halos using both the SNR‑induced and anisotropic diffusion frameworks. The best‑fit parameters reproduce the observed γ‑ray profiles and spectra, yielding D_in ≈ 10⁻³ kpc² Myr⁻¹ for the SNR case and D_⊥ ≈ 10⁻⁴ kpc² Myr⁻¹ for the anisotropic case.

Future Directions
Planned extensions include: coupling the electron transport to a self‑consistent MHD turbulence evolution (e.g., resonant streaming instability), implementing magnetic‑mirror diffusion with time‑dependent transition from superdiffusion, and adding GPU acceleration for large‑scale parameter scans. The authors also intend to provide a Bayesian inference wrapper (e.g., via emcee) to enable systematic fitting of halo populations as the number of detected halos grows with CTA, LHAASO, and future MeV–TeV observatories.

Conclusion
PHECT delivers a versatile, computationally efficient platform for pulsar‑halo studies. By supporting multiple transport scenarios, handling non‑uniform grids, and requiring only a YAML configuration, it lowers the barrier for researchers to explore the rich physics of localized CR propagation, test competing theories of slow diffusion, and prepare for the influx of high‑precision γ‑ray data expected in the coming decade.