Analytical modeling of the gamma-ray morphology of TeV halos around moving pulsars

Introduction

The first TeV halos were discovered in 2017 by the High-Altitude Water Cherenkov Observatory (HAWC), which reported [1] extended very-high-energy (VHE, $>0.1$ TeV) gamma-ray emission, at photon energies of approximately 8–40 TeV, around the nearby middle-aged pulsars Geminga and PSR B0656+14 (Monogem). Both sources are located at distances of a few hundred parsecs from Earth, and their gamma-ray emission extends over approximately 20–30 pc. Subsequent observations by the Large High Altitude Air Shower Observatory (LHAASO) have revealed additional extended gamma-ray sources associated with other middle-aged pulsars [2, 3], and a number of analyzes have been devoted to search for TeV-halo-like structures around this class of astrophysical accelerators [17, 18]. The natural interpretation of this emission is inverse-Compton scattering (ICS) by electrons and positrons released by the pulsars and propagating in the surrounding interstellar medium (ISM). At multi-TeV electron energies, the Klein–Nishina suppression makes the cosmic microwave background (CMB) the dominant target radiation field, unless a sufficiently intense local infrared field is present [see e.g. 4], while at lower energies, other components of the interstellar radiation field (ISRF) may also contribute [5]. Since the CMB and diffuse ISRF provide uniform target photon fields, the gamma-ray morphology at a given photon energy largely traces the spatial distribution of the leptons within the characteristic energy range that mostly contributes to the corresponding inverse-Compton emission.

One of the most striking results of the HAWC observations was the inference of a diffusion coefficient approximately two to three orders of magnitude smaller than the average Galactic value inferred from cosmic-ray (CR) measurements [29], most notably the boron-to-carbon ratio, and extrapolated to multi-TeV energies. Indeed, within the interpretation of the emission in terms of isotropic, loss-limited diffusion, the observed extension of TeV halos requires $D(100\,{\rm TeV})\sim10^{27}$–$10^{28}\,{\rm cm^2\,s^{-1}}$. The possible presence of regions of suppressed diffusion around CR accelerators has important implications for particle transport and may influence the locally measured CR fluxes [31]. Such regions may also affect the diffuse gamma-ray emission and the contribution of pulsars to the CR positron excess [see 8, 7 for reviews].

Several explanations have been proposed for these findings. Some studies have focused on explaining suppressed diffusion in terms of the properties of the pulsar environment, which may be more turbulent than the Galactic average [9, 10], or in terms of the amplification of magnetic turbulence by streaming instabilities driven by the escaping electron–positron pairs [15, 16]. Other studies have questioned the assumption of isotropic suppressed diffusion, considering, for instance, the possible impact on gamma-ray emission of the transition from quasi-ballistic to diffusive propagation inherent to any diffusion scenario [11, 12], or, more importantly, the effect of anisotropic diffusion with respect to the large-scale magnetic field in the source region [9, 14, 13]. In this scenario the parallel diffusion coefficient can remain close to typical Galactic values, while diffusion perpendicular to the local field is strongly suppressed. This possibility has recently received renewed attention in light of new analyses of the Geminga and Monogem halos [25].

The implications of suppressed diffusion around pulsars for the CR positron excess have also been extensively investigated, often adopting two-zone diffusion models in which a region of inhibited diffusion surrounds the source, while particle propagation at larger distances is governed by the typical Galactic diffusion coefficient [22, 23, 21, 19, 20, 26].

An additional ingredient that may substantially affect the TeV-halo morphology is the pulsar proper motion. Due to natal kick velocities that can exceed several hundred km/s [27], pulsars can travel appreciable distances during their lifetimes, but also during the time over which the emitting particles lose their energy through synchrotron and inverse-Compton losses. Consequently, the observed emission results from the superposition of particle populations injected at different positions along the pulsar trajectory [21], potentially producing asymmetric halos with emission extending behind the moving pulsar along its trajectory. This effect has already been investigated in the case of Geminga, including analyses of its lower-energy gamma-ray emission with Fermi-LAT [19], and in subsequent studies of the morphology of gamma-ray halos around moving pulsars [24].

The possible relevance of the pulsar proper motion can be appreciated by comparing the characteristic diffusion distance of the particles with the distance traveled by the pulsar during their radiative loss time. For particles of energy $E$, these scales are approximately

where $D(E)$ is the diffusion coefficient, $v_{\rm psr}$ is the pulsar velocity, and $\tau_{\rm loss}(E)$ is the energy-loss time. For particles with energies $E \approx 10$–$200$ TeV, the latter is dominated by ICS losses on the CMB and synchrotron losses in the background magnetic field, giving a loss time [6]:

Assuming a Geminga-like diffusion coefficient, $D(100\,{\rm TeV})\sim 5\times10^{27}\,{\rm cm^2\,s^{-1}}$, electrons with energies $E\sim100$ TeV would diffuse over a distance of $\sim 20 \,\rm pc$, while a pulsar moving at $v_{\rm psr}\sim 100\,{\rm km\,s^{-1}}$ would travel about 1 pc. Larger pulsar velocities, which are not uncommon, can appreciably increase this displacement, as in the case analysed in [3], where a proper-motion velocity of $\sim 500 \, \rm km/s$ was reported. On the other hand, for a diffusion coefficient that increases with energy as $D(E)\propto E^\delta$, with typical values $0.3 \lesssim\delta \lesssim 0.7$ [see e.g. 29], the ratio $r_{\rm psr}/r_d$ tends to increase at lower energies. A similar trend is expected for particles whose loss time exceeds the pulsar age, in which case $\tau_{\rm loss}(E)$ is replaced by $t_{\rm age}$. Thus, depending on the pulsar velocity and the diffusion coefficient, proper motion can significantly affect the halo morphology, producing asymmetries that could be detectable in gamma-ray observations [24].

Although the morphology of gamma-ray halos around moving pulsars has already been investigated [19, 24], the inclusion of proper motion in gamma-ray analyses generally requires evaluating the particle distribution resulting from continuous injection along the pulsar trajectory [21], integrating it along the line of sight (LOS), and convolving the resulting projected distribution with the ICS kernel to obtain the gamma-ray surface brightness. This procedure can be computationally expensive, potentially limiting the inclusion of proper motion in extensive studies of TeV halos with current and future observations. The aim of this work is to derive an analytical prescription for the gamma-ray morphology of TeV halos around moving pulsars, avoiding computationally expensive numerical integrations. Our goal is to provide a practical means of constructing spatial templates for gamma-ray observations, extending the analytical approaches commonly adopted for stationary TeV halos [see, e.g., 1, 28].

We begin by re-deriving the three-dimensional lepton distribution around a moving source, following [21]. This solution involves a time integral over the particle injection history. We then show that the integration along the line of sight can be performed exactly, retaining the full geometrical dependence on the source trajectory and observer position. The projected lepton distribution is thus reduced to a single time integral, which can be evaluated numerically for arbitrary injection spectra and energy-loss rates. The expression can be further simplified in the far-field limit, when the source distance is much larger than the characteristic size of the halo.

Finally, adopting the far-field limit, we assume a constant injection luminosity and a simplified treatment of radiative energy losses, following the prescription introduced by [6] and commonly adopted in TeV-halo analyses [1]. Under these assumptions, we show that the remaining time integral can also be approximated analytically, yielding a fully analytical expression for the projected lepton distribution. The resulting solution naturally factorizes into an energy-dependent spectral term, containing the injection spectrum and diffusion coefficient, and a dimensionless spatial term describing the energy-dependent morphology and its dependence on the diffusion coefficient and pulsar velocity.

To facilitate their application, we provide a Python implementation of the exact and approximate solutions, together with an illustrative example, as supplementary material.

The paper is organized as follows. In Sect. § 2, we derive the three-dimensional lepton distribution around a moving pulsar. In Sect. § 3, we perform the exact line-of-sight integration, introduce the far-field approximation, and derive a fully analytical expression for the projected lepton distribution under the simplifying assumptions discussed above. In Sect. § 4, we assess the accuracy of the approximations and illustrate their application to the gamma-ray morphology of TeV halos, using LHAASO J0248+6021, associated with PSR J0248+6021 [3], as a reference case. Finally, we summarize our results in Sect. § 5.

Electron-positron distribution: continuous injection from a moving source with diffusion and losses

The distribution function for leptons instantaneously released at time $t_{\rm inj}$ from a point source located at position $\mathbf{x}_s$ and observed at time $t_{\rm obs}$ reads [6]

where $b(E)=-dE/dt$ is the positive energy-loss rate, $E_{\tau}$ is the initial energy of particles cooling down to energy $E$ over a time $\tau=t_{\rm obs}-t_{\rm inj}$, and $r_d$ is the diffusion-loss radius, defined by

We model a pulsar of age $t_{\rm age}$, moving at a uniform velocity $\mathbf{v}$, as a sequence of burst-like sources injecting particles at times $t_i$ and positions $\mathbf{x}_{s,i} = \mathbf{x}_{\rm age} + \mathbf{v}(t_i - t_{\rm age})$, where in the latter equality we assumed that the pulsar is located at position $\mathbf{x}_{s, \rm age}$ at present time.

The resulting distribution is the sum of this individual contributions,

where $\tau_i=t_{\rm age}-t_i$. In addition, we assume that the source injects particles with a time-dependent spectrum $q(E,t)$, that can be normalized to the injected particle luminosity $L_e(t)$ according to $\int q(E,t)E\,dE=L_e(t)$. For each injection interval $\Delta t$, the corresponding burst-like injection spectrum is $Q(E_{\tau_i},t_i)=q(E_{\tau_i},t_i)\Delta t$.

Substituting this expression into Eq. (2) and taking the limit $\Delta t\to 0$, the sum becomes an integral over the injection time. We then change the integration variable from $t$ to $\tau=t_{\rm age}-t$ and obtain

The energy loss rate due to synchrotron and inverse-Compton losses, with the addition of ionization, Coulomb and Bremsstrahlung losses, is calculated following [6], while the Klein-Nishina effect in ICS is taken into account following [5].

Derivation of the gamma-ray morphology

The gamma-ray surface brightness is obtained by integrating the particle distribution of Eq. (3) over the radiative production rate and along the LOS:

where $\mathbf{x}_{\rm los}=\mathbf{x}_{\rm los}(s,\mathbf{x}_{\perp})$ denotes a position along the LOS, $\mathbf{x}_{\perp}$ is the projected position on the plane of the sky at the source distance, and $\mathrm{PR}_{\gamma}$ is the differential $\gamma$-ray production rate.

For inverse-Compton emission, the target radiation fields, particularly the CMB, can be approximated as spatially uniform. Under this assumption, the spatial dependence of the emission is entirely determined by the particle distribution, and the surface brightness can be written as

where $\Sigma (t_{\rm age},\mathbf{x}_{\perp},E)$ is the LOS-integrated particle distribution:

with

To evaluate the LOS-integrated particle distribution, we consider the geometry illustrated in Figure 1 and derive three expressions, whose detailed derivations are presented in Appendix Figure 6:

(i) Exact LOS integration: an exact expression for the LOS integral, that retains the full dependence on the pulsar and observer positions, the direction of pulsar motion, the injection history, and radiative energy losses. The resulting projected distribution requires only a one-dimensional numerical integration over the particle injection history.

(ii) Far-field approximation: a simplified expression for the projected distribution, valid when the source distance is much larger than the characteristic halo extension, while retaining the full treatment of injection history and energy losses.

(iii) Fully analytical approximation: an approximate expression for the remaining time integral, valid in the far-field limit under the additional assumptions of constant injection luminosity and a simplified treatment of energy losses.

Exact expression for the line-of-sight integrated particle distribution

As shown in Appendix Figure 6, the integration along the LOS (variable $s$) in Eq. (6) is analytic and gives

with

Here, $\mathbf{n}$ is the unit vector along the LOS, pointing from the observer towards the position $\mathbf{x}_{\perp}$ on the plane of the sky.

This expression is exact and valid for a generic position of the observer $\mathbf{x}_{\rm obs}$ relative to the current position of the pulsar $\mathbf{x}_{\rm age}$. It can be combined with Eq. (5) to compute the $\gamma$-ray surface brightness for a generic injection spectrum and time-dependent injection luminosity $L_e(t)$, using packages such as NAIMA [30] to evaluate the inverse-Compton emission.

Far-field approximation for the line-of-sight integrated particle distribution

Without loss of generality, and referring to the geometry illustrated in Figure 1, we place the observer along the $z$-axis, namely $\mathbf{x}_{\rm obs} = d\,\hat e_z$, and assume that the pulsar moves along the $x$-axis and is currently located at the origin, namely $\mathbf{x}_{\rm age} = \mathbf{0}$ and $\mathbf{v}\tau = v\tau\,\hat e_x$. We also introduce polar coordinates to identify a position in the plane of the sky with respect to the current pulsar position

With these choices and assuming the far-field limit, $\rho\ll d$, $v\tau\ll d$ and $r_d \ll d$, one can show (see Appendix Figure 6) that the LOS-integrated particle distribution of Eq. (8) simplifies to

Analytical approximation for constant injection luminosity and simplified energy losses

We now seek a fully analytical approximation to the LOS-integrated particle distribution. Following the treatment of radiative losses proposed by [6], we assume that particles of energy $E$ contribute to the distribution over a maximum time

where $\kappa\simeq0.75$. Within this interval, we approximate $E_\tau\simeq E$ and $r_d^2(E,E_\tau)\simeq4D(E)\tau$. Note that the same assumption leads to the well-known approximation of [6] for the electron distribution around a stationary source with constant injection luminosity and isotropic diffusion:

This expression has been widely adopted in analyses of TeV halos [1]. Assuming a constant injection luminosity, one gets $A(E, \tau)\simeq q(E)$ and the far-field expression of Eq. (11) reduces to

where we have introduced the characteristic timescales

These represent, respectively, the diffusion time over a distance $\rho$, the time required for the pulsar to travel the same distance, and the characteristic timescale associated with the competition between diffusion and source motion. As shown in Appendix Figure 6, the time integral in Eq. (14) can be computed analytically in the cases of negligible proper motion $v\to0$ and of infinite $\tau_{\rm max}$, giving

and

where $E_1$ is the first-order exponential integral and $K_0$ is the modified Bessel function of the second kind.

Note that the limit $v \to 0$ recovers the projected particle distribution of a spherically symmetric TeV halo centered on a stationary source, providing an analytical parameterization that can also be used to model such systems. For a finite $\tau_{\rm max}$ and a generic $v$, in Appendix Figure 6 we show that the time integral can be approximated analytically, giving

where

Note that the result of Eq. (18) naturally splits into a spectral component $\mathcal{S}(E)$ and a dimensionless spatial factor $\mathcal{T}(\rho,\phi,E)$ that depends on energy through the diffusion coefficient $D(E)$ and $\tau_{\rm max}(E)$:

Equilibrium lepton spectrum and total gamma-ray flux

The spatially-integrated lepton spectrum is obtained by integrating the projected lepton density over the entire plane of the sky. This corresponds to the equilibrium lepton spectrum, as shaped by injection and energy losses. As shown in Appendix app:equilibrium-spectrum, this yields

The result is independent of the diffusion coefficient and pulsar velocity, since spatial transport redistributes the particles without changing their total number at a given energy. In the simplified treatment, this reduces to

The total gamma-ray flux produced by inverse-Compton scattering at distance $d$,

can therefore be calculated directly from the equilibrium lepton spectrum $N_{\rm eq}(E)$, for instance using NAIMA.

Validation of the approximations

In this section, we assess the accuracy of the analytical approximations introduced in the previous sections by comparing their predictions with the exact calculation of the projected electron distribution reported of Eq. (8).

As a representative case, we consider the extended multi-TeV gamma-ray source LHAASO J0248+6021 [3], associated with PSR J0248+6021 (spin-down power $\dot{E} \approx 2\times10^{35}\,\mathrm{erg\,s^{-1}}$, age $t_{\rm age} \approx 60\,{\rm kyr}$, distance $d\sim 2\,{\rm kpc}$). This system is particularly interesting for investigating the effects of pulsar proper motion, owing to the relatively high pulsar transverse velocity, $v_{\rm psr} \sim 500\,\mathrm{km\,s^{-1}}$, and the diffusion coefficient inferred from the LHAASO observations, $D(E) \approx 2\times10^{28}(E/160\,\mathrm{TeV})^{0.5} \,\mathrm{cm^2\,s^{-1}}$. The combination of substantial pulsar velocity and relatively slow particle diffusion makes the effects of proper motion potentially relevant for the spatial distribution of the gamma-ray emission.

Energy losses include synchrotron emission in a magnetic field of $3\,\mu\mathrm{G}$ and inverse-Compton scattering on the CMB and on the interstellar radiation field components considered in [5].

The corresponding LHAASO measurements are included as an observational reference for the characteristic spectral and spatial properties of the emission. We emphasize that the adopted parameters are not obtained through a fit to the data, and that the purpose of this exercise is to test the approximations rather than to construct a detailed model of the source. We assume continuous injection of electron–positron pairs with a power-law spectrum and an exponential cutoff, adopting a constant injection luminosity of $L_e = 2\times10^{34}\,\mathrm{erg\,s^{-1}}$, corresponding to approximately $10\%$ of the current pulsar spin-down power, a spectral index of $1.2$, and a cutoff energy of $E_{\rm cut}=50\,\mathrm{TeV}$. These parameters are chosen to approximately reproduce the overall spectral energy distribution (SED) reported by LHAASO, following the discussion of Sect. § 3.4 and without performing a detailed fit. Figure 2 compares the gamma-ray SEDs obtained from the exact and approximate calculations, together with the best-fit spectral model reported by LHAASO for a log-parabola function, shown for reference. The two calculations are in good overall agreement, although small differences emerge around the spectral cutoff. These differences arise from the simplified treatment of energy losses in the approximate calculation, which does not fully reproduce the equilibrium electron spectrum in this energy range.

In Figure 3, we compare the projected electron distributions obtained with the exact calculation, the far-field approximation, and the analytical approximation. The comparison is performed at representative electron energies between $0.1$ and $100\,\mathrm{TeV}$, and along different azimuthal directions relative to the pulsar proper motion.

As discussed in Sect. Figure 6, the accuracy of the analytical approximation of the time integral depends on the parameter

Good agreement is found for $\eta < 10$, with particularly high accuracy for $\eta < 0.1$. Although deviations increase for larger values of $\eta$, they remain below 10% up to $\eta \sim 100$, as shown in Sect. Figure 6. Larger discrepancies occur only in regions where the integral is already close to zero, and therefore have little impact on the resulting particle distribution. For the reference values considered here, $\eta$ can significantly exceed unity, particularly at small projected distances and high electron energies. However, even in these regimes, the approximation generally remains within the range of good accuracy identified above, with larger deviations limited to regions where the projected particle density is strongly suppressed.

As expected, the projected electron distribution is asymmetric, extending to larger distances behind the moving pulsar. The three calculations show remarkably good agreement at energies below $\sim 1\,\mathrm{TeV}$ across the entire considered spatial range. At higher energies, deviations become more pronounced, particularly at projected distances exceeding $\sim 10\,\mathrm{pc}$.

In general, the analytical approximation becomes less accurate at large distances and at energies approaching the cutoff of the injected spectrum. As already illustrated by the comparison of the total gamma-ray SEDs, these discrepancies are primarily associated with the simplified treatment of electron energy losses, which becomes less reliable near the spectral cutoff.

In Figure 4, we show the inverse-Compton gamma-ray surface-brightness maps at representative photon energies between $0.1$ and $100\,\mathrm{TeV}$, calculated using the exact LOS-integrated electron distribution. The emission is computed using the same radiation fields adopted in the calculation of electron energy losses. The maps clearly shows an asymmetric morphology, with the emission extending preferentially behind the moving pulsar. This asymmetry is more pronounced at lower photon energies, where longer electron cooling times and smaller diffusion coefficients allow particles injected at earlier stages of the pulsar evolution to retain a stronger spatial imprint of its past positions.

To further investigate the spatial distribution of the gamma-ray emission, we consider the surface brightness integrated over the photon-energy range $25$–$200\,\mathrm{TeV}$. In Figure 5, we compare the resulting radial profiles with the best-fit one-dimensional gamma-ray flux distribution reported by LHAASO [3] for $E_{\gamma}> 25~\rm TeV$, obtained using their diffusion model. Also in this case, the LHAASO profile is shown for reference only, and no attempt is made to fit the observational data.

For the illustrative parameters adopted here, the exact calculation provides a reasonable description of the normalization and angular dependence of the LHAASO reference profile. The analytical approximation exhibits larger deviations, consistent with the discrepancies already identified in the total gamma-ray SED at energies approaching the spectral cutoff. Interestingly, these differences primarily affect the overall normalization, while the radial dependence is reasonably well reproduced. This is illustrated by the additional curve obtained by dividing the integrated surface brightness predicted by the analytical approximation by a factor of $\sim 1.6$, which brings it into closer agreement with the exact calculation. This suggests that, in the energy range considered, the main limitation of the analytical approximation lies in the simplified treatment of electron energy losses rather than in the description of the spatial morphology. Such a discrepancy could introduce systematic uncertainties in the inferred particle injection parameters, particularly in the lepton injection efficiency relative to the pulsar spin-down power.

Finally, we note that the predicted profiles have not been convolved with the instrumental point-spread function (PSF). Consequently, the agreement with the LHAASO reference profile should be interpreted as illustrative, and not as a quantitative assessment of the source morphology.

Conclusions

We have derived analytical expressions for the projected distribution of electrons and positrons continuously injected by a moving pulsar and propagating through the surrounding medium under the effects of energy-dependent diffusion and radiative losses.

In particular, we have shown that the integration of the particle density along the line of sight can be performed exactly, without approximations on the observer and source positions, and on the particle injection history and radiative energy losses. This result reduces the calculation of the projected particle distribution to a single time integral, eliminating the need for an additional, computationally expensive spatial integration. It therefore provides an efficient and accurate way to calculate the gamma-ray surface brightness of extended sources associated with moving pulsars, while retaining the full effects of their proper motion and injection history.

We have further simplified the exact solution by introducing the far-field approximation, valid when the source distance is much larger than the characteristic size of the halo and the pulsar displacement. Under the additional assumptions of constant injection luminosity and a simplified treatment of radiative energy losses, we have derived a fully analytical expression for the projected lepton distribution, eliminating the remaining time integral. The resulting expression factorizes into an energy-dependent spectral term, containing the injection spectrum and diffusion coefficient, and a dimensionless spatial term describing the energy-dependent morphology and its dependence on the diffusion coefficient and pulsar velocity. This provides a practical analytical parameterization of the asymmetric emission around moving pulsars. In the limit of vanishing pulsar velocity, the expression also provides an analytical parameterization of the projected lepton distribution of spherically symmetric TeV halos around stationary pulsars.

By comparing the exact and approximate solutions, we have assessed the accuracy of the approximations and identified their main limitations, particularly at large distances from the source, where the emission is generally fainter, and at high energies, especially near the spectral cutoff, where the simplified treatment of radiative losses becomes less reliable. We have also shown that the latter primarily affects the inferred lepton injection spectrum, particularly its normalization, when the model is compared with gamma-ray observations.

We have illustrated the results through the calculation of the inverse-Compton gamma-ray emission from a representative TeV-halo source, using parameters inspired by LHAASO J0248+6021, associated to PSR J0248+6021. The comparison confirms that the exact line-of-sight integration provides a practical means of calculating energy-dependent gamma-ray morphologies, while the analytical approximations reproduce their main spatial features with reasonable accuracy. These results are particularly relevant for interpreting extended gamma-ray sources associated with pulsars, where proper motion can introduce significant departures from spherical symmetry.

To facilitate the application of these results, we provide a Python implementation of the exact and approximate solutions, together with an illustrative example, as supplementary material.