© The Authors 2026

1. Introduction

Since the discovery of TeV haloes around the Geminga and B0656+14 pulsars by the HAWC experiment in 2017 (Abeysekara et al. 2017), and around J0622+3749 by LHAASO in 2021 (Aharonian et al. 2021), theoretical efforts have been dedicated to understanding this source class. These very high-energy (≳100 GeV) gamma-ray emissions are produced by the inverse-Compton scattering of leptons accelerated in the pulsar wind nebula of the pulsar with the cosmic microwave background and the interstellar photon fields (Di Mauro et al. 2019). Known haloes originate from middle-aged pulsars (∼100 kyr) expected to have escaped their parent supernova remnant (SNR) into the interstellar medium (ISM) (Abeysekara et al. 2017; Fang 2022; Amato & Recchia 2024). Surprisingly, the gamma-ray emission hints at a stronger confinement – with a suppression factor ≳100 with respect to average conditions in the Galaxy – of the high energy (≳100 TeV) particles around a region of radius 20 − 30 pc around the pulsars.

The debate regarding the number of detectable TeV haloes has remained active since their discovery. Several studies expected the number of detected pulsar haloes to grow rapidly (Linden et al. 2017; Linden & Buckman 2018). However, after almost 10 years of searches, the number of candidate TeV haloes grew to just a dozen with most requiring further spectral and morphological studies for confirmation. A high occurrence rate of TeV haloes around middle-aged pulsars also seems disfavoured by the local positron flux and the properties of the known population of nearby pulsars (Martin et al. 2022).

This discovery has consequences not only in gamma-ray astronomy but also in cosmic-ray physics. For example, before 2017, an effort was made to understand whether pulsars or dark matter could explain the anomalous positron excess in cosmic rays (Yuan et al. 2015; Boudaud et al. 2015). The detection of TeV haloes led to the idea that the large-scale cosmic-ray diffusion coefficient of the Galaxy – deduced from secondary and primary cosmic-ray nuclei fluxes, mainly probing the large diffusive halo of the Milky Way – is not representative of transport conditions in the disc. As a consequence of this lower diffusion coefficient, cosmic-ray leptons would take more time to diffuse between nearby sources and the Earth, suffer stronger energy losses particularly at the highest energies, which would make nearby pulsars unable to explain the positron excess (Abeysekara et al. 2017). Later, more refined studies showed that pulsars could still explain the positron excess if propagation around them is described in a two-zone framework, with suppressed diffusion within tens of parsec of the source – the halo – and normal standard diffusion at larger distances (Tang & Piran 2019; Schroer et al. 2023).

The theoretical efforts to explain the suppressed diffusion can be divided into three categories, depending on whether the suppression is actually necessary or considered necessary but either linked to the pulsar progenitor or not. Belonging to the first category, for example, are models of anisotropic transport along the magnetic fields surrounding the pulsar (Liu et al. 2019b; De La Torre Luque et al. 2022), or models that take into account non purely diffusive transport (Recchia et al. 2021). Belonging to the second category, are models where the suppressed diffusion is attributed to turbulence self-generated by the lepton pairs emitted by the pulsar (Evoli et al. 2018; Mukhopadhyay & Linden 2022) through the excitation of plasma instabilities. Belonging to the third category – the one relevant for the present work – include, for instance, the work by Fang et al. (2019), which proposes that pulsars exhibiting a TeV halo remain inside their parent SNRs, whose interiors could be turbulent, and the work by Schroer et al. (2022), which suggests that the turbulence responsible for the formation of the TeV halo may be self-generated by protons accelerated by the SNR and escaping the system.

The lepton pairs creating TeV haloes are also expected to emit in other wavelengths, such as radio or X-ray through synchrotron emission. Both should in principle allow us to gain new information on the environment of the pulsar, especially its magnetic field strength or structure (Liu et al. 2019a), and the emitting electron and positron population (Khokhriakova et al. 2024). In the case of Geminga, no X-ray synchrotron halo could be detected, and Manconi et al. (2024) gave an upper limit on the magnetic field strength around the pulsar.

Despite intense interest in TeV haloes, there remains no clear consensus on the environments that host TeV-halo pulsars – critically limiting our ability to pinpoint the physical mechanisms behind the observed suppression of cosmic-ray diffusion. This uncertainty is exemplified by Geminga, whose putative SNR remains undetected, and by PSR B0656+14, which appears projected within the Monogem ring, itself a candidate parent SNR. Since the nature and level of turbulence – and thus particle propagation – can vary dramatically depending on the local environment, resolving this ambiguity is essential. In this work, we take a decisive step forward: by reconstructing the past trajectory and environment of each pulsar, we directly investigate the medium traversed throughout its life, providing crucial insights into the origin and properties of TeV haloes.

van der Swaluw et al. (2003) computed the escape time t_esc at which a pulsar with speed v_kick reaches the ISM when escaping from a SNR of energy E_SN propagating in the uniform ISM of density n_ISM and found

$$\begin{array}{c}\hfill t_{\mathrm{esc}}=14{\left(\frac{E_{\mathrm{SN}}}{{10}^{51}\mathrm{erg}}\right)}^{1/3}{\left(\frac{v_{\mathrm{kick}}}{1000\mathrm{km}/\mathrm{s}}\right)}^{-5/3}{\left(\frac{n_{\mathrm{ISM}}}{1{\mathrm{cm}}^{-3}}\right)}^{-1/3}\mathrm{kyr}.\end{array}$$(1)

This leads to the hypothesis that the middle-aged pulsars creating the haloes are all in the ISM.

However, pulsars are born from massive stars that, whether isolated or in groups, blow powerful winds during their lives. The stellar winds carve rarefied, hot, and turbulent bubbles that can extend up to several hundreds of parsecs (Weaver et al. 1977; Mac Low & McCray 1988; Vieu et al. 2022). Then, the supernova explosion following the birth of the pulsar ejects stellar material inside the bubble, which further modifies the surrounding medium. Accounting for this peculiar environment challenges the usual understanding of the environment of pulsars. Stellar bubbles and their interaction with the SNR were never studied in the literature relative to TeV haloes.

In this paper, we present in Sect. 2 a model for the environment probed by pulsars following their formation and the impulse they receive in core-collapse supernova explosions of massive-star progenitors. We provide in Sect. 3 the statistical description of the environment of pulsars as a function of time, for a fully synthetic Galactic population. In Sect. 4 we apply our model to a selection of known pulsars to statistically determine the nature of the medium they are most likely to reside in, and we discuss in more detail a number of individual objects including the current sample of established TeV haloes. Finally, we conclude in Sect. 5.

2. Model

2.1. Pulsar kinematics

Pulsar kick velocities result from asymmetries in the SN explosion happening at the end of the life of the progenitor massive star (Lambiase & Poddar 2024). Assuming a uniform and rectilinear trajectory, the pulsar position depends on kick velocity v_kick and true age t_age as r_PSR = v_kickt_age.

A recent work using parallax and proper-motion measurements of radio pulsar by Igoshev (2020) suggests using single and bimodal Maxwellian distributions for the kick velocities of old and young pulsars, respectively (young meaning a spin-down characteristic age τ_c < 3 Myr). Since the middle-aged pulsars that we consider in the present work have a characteristic age comprised in the range of 100 − 300 kyr, we adopted a bimodal Maxwellian distribution f_2ℳ:

$$\begin{array}{c}\hfill f_{2\mathcal{M}}(v)\mathrm{d}v=w\mathcal{M}(v,{\sigma }_1)\mathrm{d}v+(1-w)\mathcal{M}(v,{\sigma }_2)\mathrm{d}v,\end{array}$$(2)

where v is short for v_kick, ℳ(v) is the Maxwellian distribution, and the optimal parameters are w = 0.2, σ₁ = 56 km/s, and σ₂ = 336 km/s¹.

We assessed the impact of using different kick velocity distributions on a subset of our results by considering the total sample of pulsars in Igoshev (2020) or previous distributions, such as those from Faucher-Giguère & Kaspi (2006) and Hobbs et al. (2005), adopted in the literature.

2.2. Progenitor stars

We assumed that pulsars are born from massive stars of zero-age main-sequence (ZAMS) masses M_ZAMS between 8 M_⊙ < M_ZAMS < 20 M_⊙ (Sukhbold et al. 2016). Stars in this mass range follow a power-law initial-mass function (IMF) of index −2.35 (Kroupa & Weidner 2003). Stars with ZAMS masses in the interval 2 M_⊙ < M_ZAMS < 16 M_⊙ are considered B stars, while stars with higher masses are O stars. After reaching the end of the main-sequence (MS) phase, massive stars enter the red supergiant (RSG) and Wolf-Rayet (WR) phases. As progenitor properties, we used the results of Seo et al. (2018) and references therein, which are based on stellar evolution calculations by Ekström et al. (2012) for non-rotating stars.

The MS and RSG phase durations, τ_MS and τ_RSG, depend on M_ZAMS and are given in megayears as (Zakhozhay 2013; Seo et al. 2018)

$$\begin{array}{cc}\hfill log{\tau }_{\mathrm{MS}}& \approx 9.96-3.32log{M}_{\mathrm{ZAMS}}+0.63{(log{M}_{\mathrm{ZAMS}})}^2\hfill \\ \hfill & +0.19{(log{M}_{\mathrm{ZAMS}})}^3-0.057(log{M}_{\mathrm{ZAMS}}),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill log{\tau }_{\mathrm{RSG}}& \approx -2.76log{M}_{\mathrm{ZAMS}}+9.38.\hfill \end{array}$$(4)

The fraction of the lifetime spent as RSG decreases strongly with the mass and amounts to 17% for a 8 M_⊙ B star or 7% for a 20 M_⊙ O star.

We did not account for the WR phase as it only affects stars with masses ≥25 M_⊙.

The total mass lost ΔM during each phase is given in M_⊙ by

$$\begin{array}{cc}\hfill log\mathrm{\Delta }{M}_{\mathrm{MS}}& \approx 2.36log{M}_{\mathrm{ZAMS}}-3.15,\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill log\mathrm{\Delta }{M}_{\mathrm{RSG}}& \approx 0.85log{M}_{\mathrm{ZAMS}}-0.24.\hfill \end{array}$$(6)

For the instantaneous wind mass-loss rate Ṁ (in solar masses per year) and velocity u_w (in kilometres per second) in each phase, we followed the prescription of Seo et al. (2018):

$$\begin{array}{cc}\hfill log{\dot{M}}_{\mathrm{MS}}& \approx -3.38{(log{M}_{\mathrm{MS}})}^2+14.59log{M}_{\mathrm{MS}}-20.84,\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill log{\dot{M}}_{\mathrm{RSG}}& \approx 4.16log{M}_{\mathrm{RSG}}-10.27,\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill log{u}_{\mathrm{w},\mathrm{B}}& \approx 0.21log{M}_{\mathrm{MS}}+2.85,\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill log{u}_{\mathrm{w},\mathrm{O}}& \approx 0.08log{M}_{\mathrm{MS}}+3.28,\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill u_{\mathrm{w},\mathrm{RSG}}& \approx 1.9·M_{\mathrm{RSG}}^{0.85}\mathrm{km}/\mathrm{s},\hfill \end{array}$$(11)

where the quantities are given as a function of the mass of the start at the beginning of each phase, so M_MS = M_ZAMS and M_RSG = M_ZAMS − ΔM_MS. The wind luminosity L_w in ergs per second derives from the above quantities as $L_{\mathrm{w}}=\dot{M}{u}_{\mathrm{w}}^2/2$, which yields the following for the MS:

$$\begin{array}{c}\hfill log{L}_{\mathrm{w},\mathrm{B}}\approx -3.38{(log{M}_{\mathrm{MS}})}^2+15.02log{M}_{\mathrm{MS}}+20.36,\end{array}$$(12)

$$\begin{array}{c}\hfill log{L}_{\mathrm{w},\mathrm{O}}\approx -3.38{(log{M}_{\mathrm{MS}})}^2+14.77log{M}_{\mathrm{MS}}+21.21.\end{array}$$(13)

The formula for L_w, RSG is not given as it is not directly used. As explained below, the RSG wind properties are used to determine the pressure balance in the bubble of an isolated star, which requires only u_w, RSG and Ṁ_RSG.

2.3. Structure of the surrounding medium

The majority of O and B stars are formed in star clusters (SCs) and OB associations (Wright 2020; Portegies Zwart et al. 2010). A number of massive stars are found in relative isolation in the Galactic field, either because they were created this way, or because they were ejected from clusters due to gravitational interactions, thus becoming runaway stars. In the following, we handle two different families of progenitors: isolated massive stars and SCs. The latter category gathers all types of groupings irrespective of size or compactness, and it is where most massive stars in our Galaxy are to be found. Zinnecker & Yorke (2007) state that only 2% of B stars and 10 − 25% of O stars are found isolated in the Galactic field. Motivated by the subsequent analysis of Carretero-Castrillo et al. (2023), we adopted the high-end value of 25% for runaway O stars. Due to the IMF favouring low-mass stars, the impact of the number of runaway O stars is modest, and most runaways are actually B stars. Combining these probabilities we find that 10% of the massive stars in the mass range we considered are found in isolation, while the rest is found in SCs and OB associations.

As a consequence, we considered two different types of surrounding medium for the pulsar progenitors: wind-blown bubbles (WBB) and superbubbles (SB), depending on whether the progenitor star is isolated or belongs to a cluster. In both cases, for simplicity, we assumed spherical symmetry for the environment around a star or cluster located at the centre. We therefore neglected effects such as the proper motions of isolated massive stars or the substructures in stellar clusters, both potentially giving rise to asymmetric developments of the environment.

2.3.1. Wind-blown bubbles around isolated massive stars

The interaction of supersonic stellar winds from isolated massive stars with a typically warm ISM of density n_ISM generates what we call here a WBB. This results in a stratified bubble-like structure (Weaver et al. 1977). In the immediate surroundings of the star, freely expanding wind of speed u_w ends in a wind termination shock (WTS) at radius r_w. After the WTS, one finds a bubble of hot rarefied plasma with density n_b and temperature T_b made of shocked wind and evaporated interstellar material that extends up to radius r_b and occupies most of the volume of the structure. The bubble is bounded by a thin and dense shell of density n_shell found at radius r_shell, which consists of all the ISM material that was swept by the bubble during its expansion.

At the end of the MS phase, massive stars enter the RSG phase for ≳100 kyr and lose most of their mass in powerful and dense but slow winds. These slow winds and the short duration of the RSG phase are such that they do not impact the extent of the bubble, but they set a new pressure balance in its innermost regions. We therefore assumed that the size of the bubble at the time of the supernova is set from the MS phase only, while the pressure balance is determined from the RSG wind.

The WTS radius r_w was obtained by equating the ram pressure of the wind during the RSG phase and the thermal pressure of the bubble interior at the end of the MS phase:

$$\begin{array}{c}\hfill r_{\mathrm{w}}={\left(\frac{{\dot{M}}_{\mathrm{RSG}}{u}_{\mathrm{RSG}}}{4\pi \mu n_{\mathrm{b}}{\mathrm{k}}_{\mathrm{B}}{T}_{\mathrm{b}}}\right)}^{1/2}.\end{array}$$(14)

In general, a star of a certain mass in our range exhibits more powerful winds in the RSG than the MS phase, so the WTS positions we obtained (of the order of 10 pc) are several times the WTS positions obtained from MS winds (of order 1 pc) as in a number of other similar studies (e.g. Dwarkadas 2005). The bubble density n_b and the bubble temperature T_b were taken from Castor et al. (1975) and Weaver et al. (1977), who model the evaporation of the shell matter into the hot bubble:

$$\begin{array}{cc}\hfill n_{\mathrm{b}}& =8\times {10}^{-3}{\mathrm{cm}}^{-3}{\left(\frac{{\zeta }_{\mathrm{L}}}{0.22}\right)}^{6/35}{\left(\frac{L_{\mathrm{w},\mathrm{MS}}}{{10}^{36}\mathrm{erg}/\mathrm{s}}\right)}^{6/35}\hfill \\ \hfill & \times {\left(\frac{n_{\mathrm{ISM}}}{{\mathrm{cm}}^{-3}}\right)}^{19/35}{\left(\frac{{\tau }_{\mathrm{w}}}{\mathrm{Myr}}\right)}^{-22/35},\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill T_{\mathrm{b}}& =1.46\times {10}^6\mathrm{K}{\left(\frac{{\zeta }_{\mathrm{L}}}{0.22}\right)}^{8/35}{\left(\frac{L_{\mathrm{w},\mathrm{MS}}}{{10}^{36}\mathrm{erg}/\mathrm{s}}\right)}^{8/35}\hfill \\ \hfill & \times {\left(\frac{n_{\mathrm{ISM}}}{{\mathrm{cm}}^{-3}}\right)}^{2/35}{\left(\frac{{\tau }_{\mathrm{w}}}{\mathrm{Myr}}\right)}^{-6/35}.\hfill \end{array}$$(16)

The wind time τ_w contributing to the expansion of the bubble is defined as τ_w = τ_MS + τ_RSG.

We added a corrective factor ζ_L on the luminosity to account for the effect of radiative losses, possibly in relation to mixing and turbulence at the bubble-shell interface. Radiative losses are not included in Weaver et al. (1977) for the evolution of WBBs, leading to bubble sizes bigger than observed (Yadav et al. 2017). Lacking a value specific to WBBs and expecting the behaviour to be similar to SBs, we adopted the value of ζ_L = 22% proposed in Härer et al. (2023) and based on the observation of SBs.

The bubble radius r_b is given in parsecs from Weaver et al. (1977):

$$\begin{array}{c}\hfill r_{\mathrm{b}}=21\mathrm{pc}{\left(\frac{{\zeta }_{\mathrm{L}}}{0.22}\right)}^{1/5}{\left(\frac{L_{\mathrm{w},\mathrm{MS}}}{{10}^{36}\mathrm{erg}/\mathrm{s}}\right)}^{1/5}{\left(\frac{{\tau }_{\mathrm{w},\mathrm{MS}}}{1\mathrm{Myr}}\right)}^{3/5}{\left(\frac{n_{\mathrm{ISM}}}{1{\mathrm{cm}}^{-3}}\right)}^{-1/5}.\end{array}$$(17)

The wind density profile is obtained by ${\rho }_{\mathrm{w}}=\dot{M}/4\pi u_{\mathrm{w}}{r}^2$. We recall that the mass ρ and number n densities, used interchangeably throughout this paper, are linked by ρ = μm_pn with μ the molecular fraction taken at 1.4 and m_p the proton mass.

After the bubble, the shell contains all the ISM matter that was in the volume of the bubble at the star’s birth and that was swept up by the expanding bubble. The shell width is assumed to be Δr_shell = 1 pc, and the volume of the shell lies between r_b and r_shell = r_b + Δr_shell:

$$\begin{array}{c}\hfill n_{\mathrm{shell}}=n_{\mathrm{ISM}}\frac{r_{\mathrm{b}}^3}{r_{\mathrm{shell}}^3-r_{\mathrm{b}}^3}·\end{array}$$(18)

The shell and ISM are both considered a warm medium, so from Recchia et al. (2022), we take the corresponding temperatures T_shell = T_ISM = 8000 K. We provide in Table 1 the properties of the WBB at the time of supernova for two stellar progenitors of initial mass M_ZAMS = 8 M_⊙ and M_ZAMS = 20 M_⊙.

2.3.2. Superbubbles around massive SCs

Superbubbles can be described in the same theoretical framework as WBBs (Weaver et al. 1977; Härer et al. 2023), replacing the wind properties of a single star by the collective mechanical power of all stars of the cluster. As a result, SBs extend over hundreds of parsecs instead of tens of parsecs for WBBs.

In this work, we neglected any structure in the stellar cluster or its spatial extent and assumed that SBs are the result of a mechanical outflow emanating from a central point. For each SB, we first sampled the total cluster mass M_cl between 10³ M_⊙ ≲ M_cl ≲ 10⁵ M_⊙ from the following young massive SC IMF (Portegies Zwart et al. 2010),

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}{M}_{\mathrm{cl}}}\propto M_{\mathrm{cl}}^{-2}exp(-M_{\mathrm{cl}}/M_{\mathrm{cl},\mathrm{cut}}),\end{array}$$(19)

with the cutoff mass M_cl, cut = 2 × 10⁵ M_⊙ for spiral galaxies. We then populated the SC with stars of mass M_ZAMS taken from a Salpeter IMF with index α = 2.35,

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}{M}_{\mathrm{ZAMS}}}\propto M_{\mathrm{ZAMS}}^{-\alpha }exp(-M_{\mathrm{ZAMS},\mathrm{cut}}/M_{\mathrm{ZAMS}}),\end{array}$$(20)

for masses between 0.1 M_⊙ < M_ZAMS < 150 M_⊙, and adopted a cutoff mass of M_ZAMS, cut = 0.35 M_⊙ (Larson 1998). We sampled from this distribution until the cluster’s mass was reached.

The SC properties needed for modelling the SB are its mechanical luminosity L_w, SC, its collective wind mass loss rate Ṁ_SC, and speed u_w, SC. The first two quantities are simple sums of the luminosity L_w and mass loss rate Ṁ of individual massive stars i:

$$\begin{array}{cc}\hfill L_{\mathrm{w},\mathrm{SC}}& ={\displaystyle \sum _i}{L}_{\mathrm{w},i},\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill {\dot{M}}_{\mathrm{SC}}& ={\displaystyle \sum _i}{\dot{M}}_i.\hfill \end{array}$$(22)

In contrast, the SC collective wind speed u_w, SC is linked to the wind speeds u_w and mass loss rates of individual stars using momentum conservation as in Menchiari et al. (2024):

$$\begin{array}{c}\hfill u_{\mathrm{w},\mathrm{SC}}=\frac{{\sum }_i{u}_{\mathrm{w},i}{\dot{M}}_i}{{\dot{M}}_{\mathrm{SC}}}·\end{array}$$(23)

In 1D, all the stellar winds contribute fully to the wind luminosity, while in reality some energy is lost through wind interactions. The losses remain marginal when considering the bubble size. This was shown in Vieu et al. (2024) where the bubble still evolves by following the Weaver et al. (1977) model, even if they show an absence of WTS around the OB association. As stellar wind properties of individual stars, we adopted a different approach than that followed in Sect. 2.3.1 for the WBBs, for the sake of simplicity (because the superposition of different winds from a high number of stars at different evolutionary stages is a complex problem). First, we considered O and B stars with masses > 2 M_⊙ (and not just for 8 M_⊙ < M_ZAMS < 20 M_⊙, which was motivated by the ability to form a pulsar), as this broader mass range yields a better account of the total energetics of the cluster. Second, we restricted the wind properties to those of the MS stage, which we picked from Seo et al. (2018). Last, we assumed that the collective mechanical power of the SC remains constant and equal to its initial value, when all its member stars are in the MS stage. Thus, we did not account for any change in luminosity when stars of the cluster leave the MS and enter RSG or WR stages, and we also missed the decrease in luminosity when the massive stars exploded. Conversely, we did not consider the energy input from SN explosions. Nevertheless, this contribution was assessed in Vieu et al. (2022), and the authors showed that the mechanical power of a cluster is roughly constant over the first tens of millions of years, which suggests that our assumption of a constant luminosity is reasonable. A more realistic description of a SB including these effects would have to rely on hydrodynamical simulations of the time-dependent outflows from an ageing star population, but this complex modelling goes beyond the scope of our statistical work.

Finally, Parizot et al. (2004) found that SCs are born in giant molecular clouds (GMCs) with a typical gas number density n_GMC ∼ 100 cm⁻³. We took this value for the density of the ambient medium surrounding the SB, thus setting n_ISM = n_GMC.

Out of all the stars able to produce pulsars (we recall 8 M_⊙ < M_ZAMS < 20 M_⊙) in the cluster, we randomly picked one and computed its MS time. Contrary to the WBB case, when the selected massive star that produces a pulsar dies, the rest of the cluster can still act as an engine for the SB; therefore, we evolved the system after the SN of the selected star until the pulsar reached the shell and escaped.

2.4. Supernova remnant evolution

At the end of its life, the progenitor massive star explodes as a core-collapse SN, which leads to the formation of a pulsar with a natal kick velocity and the ejection of stellar material at supersonic speeds. A SNR results from the expansion and interaction of the latter in the aforementioned parent environments. We here describe how the SNR dynamics was computed in our model.

The SNR parameters are the SN energy, for which we used the fiducial value of E_SN = 10⁵¹ erg, and the ejecta mass m_ej. It could be possible to add a dispersion on E_SN, but observationally this parameter is degenerate with the medium density, resulting in loose constraints on the parameter. Therefore, we kept it at its commonly assumed value. The ejecta mass was computed as the progenitor mass minus the mass lost in winds and the remaining pulsar mass m_PSR = 1.4 M_⊙:

$$\begin{array}{c}\hfill m_{\mathrm{ej}}=M_{\mathrm{ZAMS}}-(\mathrm{\Delta }{M}_{\mathrm{MS}}+\mathrm{\Delta }{M}_{\mathrm{RSG}}+m_{\mathrm{PSR}}).\end{array}$$(24)

There are no comprehensive analytical models for the evolution of a SNR inside a stratified environment similar to a bubble that also incorporates the radiative phase and merger. We used the thin-shell approximation for the evolution of the SNR blast wave radius and speed in the stratified environment to which we added a criterion for the merger with the medium surrounding the SNR.

The model is based on the method given in Appendix A of Ptuskin & Zirakashvili (2005) for the computation of the SNR position r as a function of time t. We recall the equations used here, beginning with the equation for the SNR mass M(r) = M_ej + M_sw(r), with M_sw(r) representing the swept-up mass until the shock position r:

$$\begin{array}{c}\hfill M(r)=M_{\mathrm{ej}}+4\pi {\displaystyle {\int }_0^r}\mathrm{d}r\prime {r\prime }^2\rho (r\prime ),\end{array}$$(25)

where ρ(r) corresponds to the density profile.

The SNR mass enters in the shock velocity u_SNR profile expression:

$$\begin{array}{c}\hfill u_{\mathrm{SNR}}(r)=\frac{\gamma +1}{2}{\left[\frac{2\alpha E_{\mathrm{SN}}}{M^2(r)r^{\alpha }}{\displaystyle \underset{I(r)}{\underbrace{{\int }_0^r\mathrm{d}r\prime {r\prime }^{\alpha -1}M(r\prime )}}}\right]}^{1/2},\end{array}$$(26)

with γ = 5/3 being the adiabatic index for a monatomic gas in the medium and α = 6(γ − 1)/(γ + 1).

Finally, the evolution of the shock position as a function of time is given by

$$\begin{array}{c}\hfill t(r)={\displaystyle {\int }_0^r}\frac{\mathrm{d}r\prime }{u_{\mathrm{SNR}}(r\prime )}·\end{array}$$(27)

The SNR speed profile can be solved analytically by assuming a density profile. We recall that the density profile we use follows:

$$\begin{array}{c}\hfill \rho (r)=\{\begin{array}{cc}{\rho }_{\mathrm{w}}(r)\hfill & r<r_{\mathrm{w}},\hfill \\ {\rho }_{\mathrm{b}}\hfill & r_{\mathrm{w}}<r<r_{\mathrm{b}},\hfill \\ {\rho }_{\mathrm{shell}}\hfill & r_{\mathrm{b}}<r<r_{\mathrm{shell}},\hfill \\ {\rho }_{\mathrm{ISM}}\hfill & r>r_{\mathrm{shell}}.\hfill \end{array}\end{array}$$(28)

We solved Eq. (26) analytically in this density profile and Eq. (27) numerically to finally find the SNR’s age as a function of its position. When and if the SNR reaches the bubble shell, it sweeps the mass emitted by the progenitor star in the form of MS and RSG winds, as well as some interstellar matter evaporated from the shell as a result of thermal conduction (Weaver et al. 1977). Yet, the shell made up of swept-up ISM material is much heavier than the SNR and effectively acts as a solid wall. As a result, the SNR collides and very rapidly loses most of its energy radiatively, as discussed both theoretically and observationally in Vink (2020), to merge with the shell. In our model framework, the SNR therefore has an impact on the bubble interior but never reaches the ISM. In case the SNR does not encounter the shell, it can merge with its surrounding medium at the position r_merge when it stops being a strong shock:

$$\begin{array}{c}\hfill u_{\mathrm{SNR}}(r_{\mathrm{merge}})=\beta c_{\mathrm{s}}(r_{\mathrm{merge}}),\end{array}$$(29)

where β = 2 (Cioffi et al. 1988) is an arbitrary parameter that gives the minimal shock strength and $c_{\mathrm{s}}(r)=\sqrt{\gamma k_{\mathrm{B}}T(r)/\mathrm{\mu }{\mathrm{m}}_{\mathrm{p}}}$ is the speed of sound in the medium.

In this analytical framework, the radiative losses of the SNR were not considered. This modelling can only reproduce correctly the ejecta dominated and Sedov-Taylor phases of the SNR. This means that the position of the SNR as a function of time is always lightly overestimated in these calculations, especially at the later times.

3. Synthetic population of pulsars

In this section, we present the results of our model on the type of environment crossed by a moving pulsar as a function of time, up to pulsar ages of a few megayears. In the following, we use the concept of ‘parent environment’ to refer to the structured medium that surrounds a pulsar at birth. This immediate-vicinity environment is determined by the nature of the progenitor star (isolated or in a cluster), as well as by the properties of the SNR resulting from the core-collapse explosion, whose expansion modifies the surrounding medium. We first introduce the properties of this parent environment in relation to our model parameters, before presenting the statistical implications for the medium experienced by pulsars along their trajectory.

3.1. Properties of the parent environment

Figure 1 shows the evolution of SNRs for the two extreme values of progenitor star mass in the range we considered, 8 M_⊙ and 20 M_⊙, along with the corresponding wind-bubble structure at the time of supernova. In both cases, the expansion of the SNR is halted by its encounter with the bubble shell, at 3 kyr and 15 kyr, respectively. This is due to the mass of the bounding shell being much higher than that of the expanding remnant, by more than one order of magnitude (see Table 1). In the WBB case, a typical pulsar with a kick velocity of 350 km/s escapes from the bubble shell and enters the ISM at 64 kyr and 130 kyr, respectively.

Fig. 1. Evolution of the SNR as a function of time in the WBB scenario, for two values of the stellar progenitor initial mass. For each case, we overplotted the positions of the wind-termination shock (dotted lines, on top of each other, see Table 1) and outer bubble radius (dashed lines) at the time of explosion. The purple line indicates the trajectory of a pulsar with kick velocity 350 km/s.

Figure 2 shows the evolution of SNRs for a progenitor star of mass 8 M_⊙ in a SB, for three different values of the cluster mass, 10³ M_⊙, 10⁴ M_⊙, and 10⁵ M_⊙. As expected, the higher the cluster mass, the larger the SB radius because of the more powerful collective wind. In these examples, rather large SB are obtained because the chosen progenitor star has a long lifetime, during which the bubble could expand to large distances. Since the parameters of the explosion are the same, the trajectory changes because of the density structure inside the bubble cavity, which is primarily determined by the cluster luminosity and age, and by the exterior density. We see in all three cases that the SNR evolves from an ejecta-dominated phase to a Sedov-Taylor phase, before merging inside the SB without even reaching the bubble shell. When this happens, a typical pulsar with a kick velocity of 350 km/s remains inside the remnant.

Fig. 2. Evolution of the SNR as a function of time in the SB scenario, for a stellar progenitor of initial mass 8 M⊙ and three different masses of the stellar cluster it belongs to. Graphical elements are similar to Fig. 1. The SNR curves turn flat when the remnant merges with the cavity interior and disappears.

We notice that the SB radius in the most massive case of 10⁵ M_⊙ is 320 pc, which is larger than the typical scale height of the gas disc of our Galaxy (Langer et al. 2014). A proper treatment would require us to model the evolution of the SB in at least 2D with a gas density gradient, as well as a proper modelling of the stellar population in the cluster but this is beyond the scope of our study. In addition, because of the cluster IMF, this case is statistically rare, so we simply neglect these considerations.

In the adopted model framework, the WBB and SB shells are the relevant boundaries marking the entry into the ISM, and not the forward shock of the SNR (as would be the case when not taking into account a modified circumstellar medium). Before a pulsar reaches the ISM, it can be either inside its parent SNR, the WBB, or the SB (owing to their small widths, we neglect the short time spent in the bounding shells of the bubbles). When the pulsar resides in the bubble medium, it is interesting to distinguish two different setups: one in which the SNR is still alive and expanding, and one in which it has dissolved in the WBB or SB interior. Yet, in our Monte-Carlo simulations and given the adopted pulsar kick velocities, the latter case appears to be more common. Given the range of progenitor masses, in the WBB case, the SNR always reaches the bubble shell before the pulsar, namely the pulsar can either be found inside the SNR or in the ISM. In the SB case, the prospects are more interesting. A pulsar can escape a still-expanding SNR and be found in the SB interior for a cluster mass 10³ M_⊙ – the most favoured by the cluster IMF – but only if it has a kick velocity typically exceeding 800 km/s and a massive progenitor of 20 M_⊙. However, the most probable situation for the SB case is that the SNR merges with the bubble cavity while the pulsar was still moving within. Observationally, such a setup could have relevant consequences, particularly in the context of pulsar haloes, since it implies the possibility of finding a young or middle-aged pulsar in what is likely a turbulent medium (the interior of the bubble), even though the parent SNR can no longer be detected. Geminga may very well be in such a situation.

3.2. Statistical description of a pulsar environment

For each family of progenitor stars described above (isolated or in cluster), we generated a population of 10 000 pulsars in their parent environments using a Monte-Carlo method. This was achieved by random sampling the initial star mass for WBB systems, the cluster and the progenitor star masses for SB systems, and the pulsar kick velocities in both cases.

From these mock populations, we computed the probability as a function of age that a pulsar is located in a particular medium among those two: the interior of the SNR and the interior of the WBB or SB. A null probability means that the pulsar has escaped into the ISM. We plotted the corresponding results in Fig. 3, using as pulsar kick velocity distribution that of Igoshev (2020) since it is based on the most complete pulsar sample. For comparison, we overplotted as vertical coloured lines the characteristic ages of the first four pulsar haloes discovered by HAWC (B0656+14 and Geminga) and LHAASO (J0622+3749 and J0248+6021), and as a shaded band the classical values of 40 − 60 kyr routinely found in the literature for the time at which a pulsar enters the ISM (e.g. Gaensler & Slane 2006; Amato & Recchia 2024).

Fig. 3. Probability for a pulsar to be found inside its parent remnant or bubble environment as a function of time. The red and blue lines correspond to the SB and WBB cases, respectively. The grey-shaded band corresponds to values typically found in the literature for the time at which a pulsar exits into the ISM, thus showing that they are significantly underestimated. The four vertical lines correspond to the characteristic ages of the confirmed TeV haloes.

Our population synthesis actually predicts that pulsars will enter the ISM somewhat later than usually assumed: at an age of 60 kyr, more than 52% of pulsars in the WBB case and all pulsars in the SB case still reside in their parent environment. As expected due to the smaller size of the structure, pulsars remain inside WBBs for a shorter amount of time compared to SBs. In the latter case, the transition from pulsars being mostly in to pulsar being mostly out occurs at an age of 200 − 500 kyr.

Both the WBB and SB cases can be combined, with the proper weighing for the global proportion of isolated versus clustered massive stars (see Sect. 2.3), to make predictions for a Galactic population. This is illustrated in Fig. 4, which resembles strongly the SB case in Fig. 3 since most massive stars live in clusters. Overall, it can be seen that before 30 kyr, almost all pulsars are inside their parent environment, while all escaped after 6 Myr. Most pulsars escape into the ISM at ∼100 − 1000 kyr, with half of the Galactic population being out at 300 kyr. This is precisely the age range of the so-called middle-aged pulsars presumably involved in TeV haloes.

Fig. 4. Probability for a pulsar to be found inside its parent environment as a function of time, for a Galactic population made of a proper mix of isolated and cluster massive-star progenitors (see text). The different curves correspond to different assumptions on the pulsar kick velocity distribution. Other graphical elements are similar to Fig. 3.

From our model and population synthesis, there is therefore no clear-cut interpretation of TeV haloes in terms of the medium a pulsar lies in, since middle-aged pulsars have comparable chances of being inside their (possibly turbulent) parent environment or in the (a priori more quiescent) ISM. Conversely, if we assume that the parent environment, provides enough turbulence for efficient trapping of ultrarelativistic electron-positron pairs, and that this is the main formation channel for TeV haloes, our model predicts a strongly decreasing occurrence rate of TeV haloes as pulsars get older than 300 kyr. At the characteristic age of Geminga, half of the pulsars would then be able to form haloes. Our model also provides the complementary information that about half of these halo-forming pulsars would still be located inside their parent SNR (see Fig. 3). This is qualitatively consistent with the proposal put forward in Martin et al. (2022) that pulsar haloes may not be so common around middle-aged pulsars.

As can be see in Fig. 2 in the SB case, the SNR can merge before reaching the bubble shell. A pulsar can be observed without a SNR surrounding it but still being inside the SB. Increasing the required shock strength before merger from β = 2 to β = 5 makes the SNR merge with the interior of the bubble marginally earlier, slightly increasing the time the pulsar spends inside the SB without probing the downstream of a SNR.

Finally, in Fig. 4, we illustrate the impact of the kick velocity distribution adopted in the Monte-Carlo sampling. We see a clear impact on the position of the transition – whether predominantly in or out of the parent environment – and, to a lesser extent, on its shape. As a consequence, the quantitative statements made above or below should be treated with some caution.

Fig. 5. Probability for a selection of 30 known pulsars within 2 kpc to be found inside their parent environments, as a function of their characteristic ages. This is based on a Galactic population model made of a proper mix of isolated and cluster massive-star progenitors (see text), and a kick velocity distribution from Igoshev (2020). Dashed and dotted lines mark the 80% and 20% probabilities, respectively.

4. Application to a selection of known pulsars

4.1. Sample construction

We applied our model to a sample of real pulsars in order to derive constraints on their environments, in a statistical sense. From the ATNF catalogue (Manchester et al. 2005), we extracted a subset of pulsar halo candidates with properties that are relatively similar to Geminga, the prototypical pulsar halo.

We considered pulsars with luminosity to distance ratio Ė/d ≥ 10³³ erg/s/kpc and characteristic age 40 kyr < τ_age < 1 Myr. The first criterion aims at jointly maximising the power and angular size of the halo emission, while the second is meant to discard that are too young systems and are more likely to be classical PWNe. We added a criterion to restrict our sample to relatively nearby objects. For a TeV halo to be detected as an extended source with current experiments, an upper limit needs to be set on the distance of the system. Assuming a typical halo size of r_halo ∼ 30 pc and considering that it would be detected as extended with HAWC and LHAASO if spanning at least three point-spread functions of α_PSF ∼ 0.3 deg, we derived an upper limit on the halo visibility at distance d_max = r_halo/tan(α_PSF) = 2 kpc. We complemented the sample with the four well-established pulsar TeV haloes (J0633+1746, B0656+14, J0622+3749 and J0248+6021) and three other objects not fulfilling the above criteria but mentioned in the literature as good halo candidates, mainly on the basis of gamma-ray observations (B0540+23, J0633+0632, and J0631+1036; see Celli & Peron 2024; Khokhriakova et al. 2024). Finally, our set of real pulsars was divided into those having estimated transverse velocity and those without, as the former category allows us to make a finer treatment and to derive richer constraints on their environments. The full sample is listed in Table A.1.

4.2. Model predictions

Before discussing the predictions of our model for the selected known pulsars, we recall that the characteristic age of a pulsar should actually be understood as an upper limit to the true pulsar age. In cases where a more accurate estimate of the age of the pulsar can be obtained by other methods, it appears that the characteristic age can markedly overestimate the true age, by up to one order of magnitude (Suzuki et al. 2021, see also the cases of PSR J0908-4913 and PSR J0538+2817 in Sect. 4.3). In the absence of a reliable estimate for the age of the pulsars in our sample, we use their characteristic age in the following discussion; nonetheless, the caveat above should be kept in mind.

From the age of selected pulsars without velocity information, we can plot Fig. 5 based on the combined WBB and SB case for a Galactic population. There, 13 (one) out of 30 pulsars in that sub-sample have a probability higher than 80% (below 20%) of being found inside their parent environment. Such an analysis can be useful to prioritise targets in searches for pulsar haloes, depending on the confinement scenario to be tested. The two halo candidates in this sample (J0631+1036, and B0540+23) and the established TeV halo (J0622+3749) have probabilities ≳65% of still residing in their parent environment. In the case of J0631+1036, the probability reaches ≳95%, making it a particularly good target for testing whether the parent environment provides the conditions for pair confinement.

The sample of pulsars with known transverse velocities is smaller but allows us to extract more information since transverse velocity is a crucial parameters in our model. We still consider for this analysis the combined WBB and SB case for a Galactic population. First, we recall that transverse velocity v_t is the projection on the sky of the kick velocity v_kick of the pulsar such that v_kick = v_t/cos(θ) with θ the angle between the two velocities. A measured v_t therefore results from any value v_kick ≥ v_t, each implying a different kinematics of the pulsar with respect to the parent environment. The probability of each v_kick is given by the distribution of natal kick velocities from Igoshev (2020), plus a weighing ∝cos(θ) to account for the fact that a uniform distribution of kick orientations in space statistically favours velocity vectors close to the plane of the sky (the allowed solid angle is larger for a given dθ element).

For each pulsar with a given characteristic age and transverse velocity, we randomly sampled a high number of possible values of v_kick ≥ v_t (based on the distributions of kick velocities and orientations mentioned above), and we computed for each v_kick the probability for the pulsar to be located in its parent environment. We then defined ranges for the most frequently encountered probabilities together with confidence intervals. The result is displayed in Fig. 6 in the form of a box plot. For each object in our sub-sample, the following graphical elements are provided: coloured symbols with black edges indicate the probability for the smallest possible velocity v_kick = v_t; boxes indicate the 68% containment interval, while whiskers extending above and below the boxes show the 95% containment interval; last, the horizontal segment inside each box is the median probability.

Fig. 6. Probability for a selection of 22 known pulsars within 2 kpc to be found inside their parent environments. Pulsars in this sample have proper motion information vt, which enters as a further constraint on their actual situation (see text). Coloured symbols with black edges indicate the probability for vkick = vt. Boxes indicate the 68% containment interval over possible kick velocities vkick ≥ vt, while whiskers extending above and below the boxes show the 95% containment interval. The horizontal segment inside each box is the median probability.

Since v_t is the smallest possible kick velocity and is sometimes much smaller than the actual v_kick, using it in our model leads to a high probability that most objects, even quite old pulsars, remain inside their parent environment. Applying the above calculation corrects for this bias and avoids overestimating this probability. As a result and in agreement with previous discussions, pulsars younger than τ_age ≲ 300 kyr still have a high probability of being inside their parent environment, while there is a transition from inside to outside between 300 kyr ≲ τ_age ≲ 600 kyr. All pulsars associated with a TeV halo or a halo candidate have high probabilities of still residing inside their parent environment according to our model.

4.3. Discussion of individual pulsars

We provide here a more detailed discussion of a number of specific pulsars and their most likely environment in view of our model’s predictions. We focussed on the pulsars in our sample associated with a PWN, because those were the targets of in-depth studies, which provide more information about the context of the system (for instance the association with a SNR, the evolutionary state or system age, or the system viewing geometry). We cross-matched our selection of pulsars with the list of PWNe provided in Olmi (2023), keeping the distinction made between normal and fast-moving pulsars (translated into the PWN and SPWN identifiers in Table A.1). Our selection includes nine pulsars with PWNe. We aimed at classifying these nine pulsars into the following two categories: those that can be located with respect to their putative parent remnant, and those that exhibit typical signatures of a propagation in the ISM (in principle, a pulsar can belong to the two categories, but this did not happen here).

The first category is well-defined, while the second category is less straightforward. The association of a pulsar with a bow-shock pulsar wind nebula (BSPWN) could in principle indicate that it is travelling in the ISM. Indeed, prototypical BSPWNe are expected to develop when pulsars propagate in a cold and mostly neutral ISM with a relatively low sound speed, and they are characterised by a double shock structure and optical line emissions generated through collisional excitation and charge exchange (Olmi 2023). Yet, the identification to a textbook BSPWN condition is far from simple. There are many circumstances that can give rise to a BSPWN appearance in X-rays or radio (high pulsar proper motion and/or recent interaction with the reverse shock of the remnant), while the typical signatures of a bow shock are not easily detectable or interpretable (complex geometry, low-surface-brightness extended emission, and lack of optical emission lines owing to upstream ionisation by the pulsar itself). Each case requires a dedicated investigation.

Six pulsars in our sample are listed in Kargaltsev et al. (2017) and Olmi (2023) as fast-moving pulsars with morphological features reminiscent of BSPWNe, and three are associated with more classical PWNe. Out of these, three display strong evidence for a connection with a SNR (J1809-2332, J0538+2817, and J0908-4913), to which we added B0656+14 (for reasons detailed below), while two seem to be travelling in the ISM (J0742-2822 and J1741-2054). We discard J0357+3205, J1740+1000, and J0358+5413 since none of the available information on these pulsars can be used in the present context (the first two are high-latitude objects that may be not representative of the pulsar population). Finally, we provide a specific discussion of J0633+1746, the Geminga pulsar and the prototypical halo.

4.3.1. Pulsars likely associated with a SNR

4.3.1.1. PSR J1809-2332

The pulsar is associated with SNR G7.5-1.7, detected as a partial radio shell (Roberts & Brogan 2008), and powers the Taz PWN. It lies in projection within the remnant, at about half a radius from the centre. The proper motion measured from a decade of X-ray observations points back towards the centre of the SNR and suggests a kinematic age of 48 kyr (Van Etten et al. 2012, about 30% less than the characteristic spin-down time). The picture is fully consistent with the prediction of our model, which places the pulsar inside its parent environment with a high probability of 95 − 100%.

4.3.1.2. J0538+2817

The pulsar is associated with SNR Sim 147 and lies well within it in projection, even accounting for the estimated total kick velocity (see Khabibullin et al. 2024a, and references therein). The angular displacement from the geometrical centre of the remnant suggests a kinematic age of 38 kyr, which is much smaller than the characteristic spin-down time of 618 kyr. We therefore adopted the former value in lieu of the latter. Interestingly, the properties of the SNR suggest that the explosion of the progenitor star took place in a wind-blown cavity (Khabibullin et al. 2024b), similar to our WBB case or to our SB case for a low-mass cluster. For such a small age, our model unambiguously places the pulsar inside its parent environment with a high probability.

4.3.1.3. PSR J0908-4913/B0906-49

Johnston & Lower (2021) updated the distance of the pulsar from 1 kpc in the ATNF catalogue to 3 kpc, which increases the transverse velocity from 240 km/s to 670 km/s. A SNR is associated with PSR J0908-4913, and the pulsar lies in projection on its edge, on the verge of escaping it. Our model predicts a ∼88% probability of being inside its parent environment. Yet, there are caveats in this comparison. First, the association with the SNR suggests a kinematical age of 12.7 kyr, significantly lower than the characteristic age of 110 kyr of the pulsar (that we used in Fig. 6). Reducing the age of the system to such a lower value would increase the probability of it being inside the parent environment up to almost 100%. On the other hand, the physical size of the putative SNR is 10 pc at a distance of 3 kpc, which is smaller than the values we obtain at around 10 kyr, even in the WBB case. This would tend to pull the probability to lower values. A dedicated modelling of this object with different parameters (lower explosion energy or wind mass loss rate for instance) would be required to get a more solid conclusion. But overall the model predicts probabilities of residing in the parent environment higher than 88%, which is consistent with observations of a pulsar in the downstream of the forward shock of its SNR.

4.3.1.4. PSR B0656+14

This pulsar lies in projection inside the Monogem ring, which has been suggested to be its parent SNR (Knies et al. 2018). With its characteristic age of 110 kyr, PSR B0656+14 has a median probability of 91% of still residing in its parent environment in our model, consistent with this picture. The fact that a TeV halo was detected around B0656+14 strongly suggests that the confinement of particles on the scale of several tens of parsecs has to do with the conditions inside the remnant or bubble. We notice that the radius of the SNR and ISM density are respectively r = 87 pc and n_H = 4 × 10⁻³ H cm⁻³, for a distance of 300 pc and using a Sedov-Taylor evolutionary solution. Such low densities can be expected from SB cavities, and it is therefore possible that the progenitor of B0656+14 exploded in a SB environment, which enabled the SNR to grow to such a large size without being stopped by the SB shell lying at ≳100 pc (see Fig. 2).

4.3.2. Pulsars likely travelling in the ISM

4.3.2.1. J0742-2822/B0740-28

This pulsar sports a clear optical bow shock (Jones et al. 2002) and, as such, could be considered to be in the predominantly neutral ISM, away from its parent environment. With τ_c = 157 kyr, PSR J0742-2822 is in the age range where the transition between mostly in to mostly out of the parent environment just starts to occurs. Accordingly, our model predicts a probability distribution from 85 to 90% (95% confidence interval). The observed stand-off distance favours small inclination of the kick velocity vector with respect to the plane of the sky (Jones et al. 2002), which would in turn favour high probabilities of residing in the parent environment. On the other hand, the morphology of the optical nebula suggests that it recently left a low-density medium (n_H ≲ 3 × 10⁻³ H cm⁻³) to enter a denser medium (n_H ≃ 3 × 10⁻¹ H cm⁻³), which can be interpreted as a transition from the interior of a bubble to the bounding shell of ambient ISM.

4.3.2.2. PSR J1741–2054

The pulsar displays a clear bow-shock nebula in optical lines (Mignani et al. 2016; Romani et al. 2010) and UV (Abramkin et al. 2025). This strongly suggests the pulsar propagates in the mostly neutral ISM, which is weakly consistent with our model prediction of an approximately median 60% probability that it resides outside its parent environment. Both Camilo et al. (2009) and Abramkin et al. (2025) noticed that the properties of J1741-2054 are very close to that of the ‘Three Musketeers’ (Geminga, B0656+14, and B1055-52). Since TeV haloes were discovered around the first two, and that B1055-52 is in the blind spot of HAWC, J1741-2054 appears as a promising target to search for a TeV halo and make a comparison with the most established instances of the phenomenon. If a halo were discovered around J1741-2054, our model would suggest pair confinement is linked to some properties of the ISM rather than be due to the parent environment.

4.3.3. The case of Geminga

With a characteristic age of 342 kyr and a transverse speed of 150 km/s, PSR J0633+1746 has travelled at least ∼53 pc since its birth. Because of this large distance, it is commonly assumed that Geminga is in the ISM.

However, our model statistically places the TeV haloes and halo candidates in our sample inside their parent environment, including Geminga. Since Geminga is usually taken as the prototypical pulsar halo, we provide three additional elements that further support the idea that it is not probing the ISM.

First, as summarised in Amato & Recchia (2024), there are many hints that Geminga is still in a hot environment. The non-detection of H_α lines in the near vicinity of the pulsar suggests that the medium Geminga traverses is hot and highly ionised. The ambient medium density inferred from Geminga’s X-ray jets is quite low, n_H ≲ 10⁻² H cm⁻³, which is typical of a bubble interior. This seems consistent with the fact that Geminga is located close to the Gemini H_α ring, a bubble that hosts an OB association. Using the formalism from Weaver et al. (1977) and comparing with Suzaku observations of the H_α ring, Knies et al. (2018) shows that the Gemini ring can be seen as a bubble blown by at least two stars of an OB association found inside. Fang et al. (2019) noticed that Geminga and the brightest part of its halo are found in projection inside the Gemini H_α ring.

Second, the morphology of X-ray jets around Geminga suggests that the kick velocity of the pulsar is nearly perpendicular to the line of sight (hence in the plane of the sky), especially when jets are interpreted as bent polar outflows (Posselt et al. 2017). In Fig. 6, this implies that the actual position of Geminga is closer to the star symbol, which implies a probability of residing in its parent environment above 90%.

Third and last, the discussion so far is based on the characteristic age of Geminga. As mentioned earlier, the characteristic age may well overestimate the true age by a factor of a few (a caveat valid for all objects in our sample). Since Geminga is exactly at the age range where most pulsars transitions from inside to outside of their parent environment, any decrease in age has significant consequences. This is illustrated in Fig. 7, where we display box plots similar to Fig. 6 for a variety of age assumptions for Geminga between 150 − 342 kyr. We see that if Geminga’s true age is roughly half of its characteristic age, its median probability of being inside its parent environment reaches ≳90%.

Fig. 7. Probability that Geminga resides in its parent environment for five different assumptions on its age. For each age, a box plot such as those in Fig. 6 is displayed.

5. Conclusions

In this work, we aimed at understanding the type of medium crossed by pulsars over their lifetime in order to constrain the physical mechanism leading to strong confinement of electron-positron pairs and ultimately to the formation of gamma-ray haloes observed at greater than or approximately equal to TeV photon energies around some pulsars. We developed a statistical model for the environment probed by pulsars over the first few megayears of their life, under two different scenarios for the circumstellar conditions: pulsars resulting from either isolated massive star progenitors surrounded by WBBs at the time of a supernova or massive stars belonging to clusters that eventually explode inside SBs. We also produced a mixed scenario from a weighted combination of these two cases to describe the population of our Galaxy.

Using a Monte-Carlo method to sample over progenitor properties and pulsar kick velocities, we computed the statistical distribution of the times at which pulsars escape their parent SNR first, and later their parent WBB or SB interior, to eventually enter into the undisturbed ISM. We applied this model to a fully synthetic population of pulsars and to a sample of known pulsars that either have the right properties to produce detectable haloes or were claimed to be halo candidates on the basis of gamma-ray observations. Our main conclusions are the following:

1. In the vast majority of cases, the parent remnant disappears before the pulsar can escape it. In WBBs, the SNR evolution is halted by its collision with the much more massive bubble shell, while in the SBs, SNRs merge with the hot interior.

2. Pulsars escape into the ISM much later than commonly assumed. Only half of a Galactic population of synthetic pulsars have left their parent environments at 300 kyr, compared to the typical escape times of 40 − 60 kyr frequently used in the literature.

3. In our model, middle-aged pulsars with ages 100 − 300 kyr, those presumably involved in TeV haloes, have higher chances of being inside their parent environment. While this does not favour any explanation for the origin of haloes in terms of the medium they develop into, this has implications on the occurrence rate of TeV haloes: if the parent environment, stellar-wind bubble or SB, provides the turbulence for efficient trapping of ultrarelativistic electron-positron pairs, our model predicts a strongly decreasing occurrence rate of TeV haloes as pulsars get older than 300 kyr. At the characteristic age of Geminga, a half of pulsars would then be able to form haloes.

4. From a sample of seven known pulsars with a confirmed or suspected halo, all of them have a significant probability of being inside their parent environment. Geminga has the highest median probability in the sample (only 35%) of being in the ISM at face value. Yet, including additional information on the kick velocity orientation and considering the uncertainty of the pulsar age increase the probability that Geminga resides in its parent environment to above 90%, consistent with the hot and ionised properties inferred for the medium around Geminga.

We provide the probability of residing inside the parent environment for a list of 53 known pulsars. Searching for extended gamma-ray emission around them with HAWC, LHAASO, and, in the future, the upcoming SWGO experiment could help us connect the occurrence of the pulsar halo phenomenon with the environment, at least from a statistical point of view. Together with complementary studies of the ISM or circumstellar medium (for instance searching for hot ionised gas in X-rays, or for neutral gas in optical lines), this should contribute to our knowledge of the mechanisms underlying the formation of TeV haloes.

Beyond these considerations, our model’s prediction that pulsars escape their parent environment much later than commonly assumed will have implications for their contribution to the flux of cosmic-ray positrons and electrons directly and locally measured. All else being equal, a delayed escape likely reduces the flux at Earth, as a smaller fraction of the pairs produced by pulsars will be effectively injected in the ISM. Accounting for the AMS-02 positron spectrum thus implies either higher pair production efficiency, contributions from a complementary source class, or different cosmic-ray transport assumptions near and/or far from the sources.

Acknowledgments

This work has made use of NASA’s Astrophysics Data System Bibliographic Services. Lioni-Moana Bourguinat acknowledges the hospitality and support of IRAP, where part of this work was carried out. Pierrick Martin acknowledges financial support by ANR through the GAMALO project under reference ANR-19-CE31-0014. The work of Sarah Recchia is funded by the European Union – “NextGenerationEU” RRF M4C2 1.1 under grant PRIN-MUR 2022TJW4EJ. The work of Carmelo Evoli has been partially funded by the European Union – “NextGenerationEU”, through PRIN-MUR 2022TJW4EJ and by the European Union – NextGenerationEU under the MUR National Innovation Ecosystem grant ECS00000041 – VITALITY/ASTRA – CUP D13C21000430001.

References

Appendix A: Properties of the individual pulsars studied

This table compiles the properties of all the individual pulsars that fit the conditions given in Sec. 4.1. The probability in the second column represents the probability of pulsars without proper motion information being found inside the parent environment, and the median probability (with the 1σ interval) for pulsars with proper motion information, in both cases corresponding to the combined WBB+SB scenario for a Galactic population. In the last column, we provide the identification as TeV halo or halo candidate, and as PWN for a normal and fast-moving pulsar (PWN and SPWN, respectively). The data shown come from the ATNF catalogue (Manchester et al. 2005) unless a different reference is provided: (1) Khabibullin et al. (2024b); (2) Van Etten et al. (2012); (3) Danilenko et al. (2019); (4) Johnston & Lower (2021); (5) Pletsch et al. (2012).

All Tables

All Figures

	Fig. 1. Evolution of the SNR as a function of time in the WBB scenario, for two values of the stellar progenitor initial mass. For each case, we overplotted the positions of the wind-termination shock (dotted lines, on top of each other, see Table 1) and outer bubble radius (dashed lines) at the time of explosion. The purple line indicates the trajectory of a pulsar with kick velocity 350 km/s.
In the text

	Fig. 2. Evolution of the SNR as a function of time in the SB scenario, for a stellar progenitor of initial mass 8 M⊙ and three different masses of the stellar cluster it belongs to. Graphical elements are similar to Fig. 1. The SNR curves turn flat when the remnant merges with the cavity interior and disappears.
In the text

Fig. 3. Probability for a pulsar to be found inside its parent remnant or bubble environment as a function of time. The red and blue lines correspond to the SB and WBB cases, respectively. The grey-shaded band corresponds to values typically found in the literature for the time at which a pulsar exits into the ISM, thus showing that they are significantly underestimated. The four vertical lines correspond to the characteristic ages of the confirmed TeV haloes.

In the text

	Fig. 4. Probability for a pulsar to be found inside its parent environment as a function of time, for a Galactic population made of a proper mix of isolated and cluster massive-star progenitors (see text). The different curves correspond to different assumptions on the pulsar kick velocity distribution. Other graphical elements are similar to Fig. 3.
In the text

Fig. 5. Probability for a selection of 30 known pulsars within 2 kpc to be found inside their parent environments, as a function of their characteristic ages. This is based on a Galactic population model made of a proper mix of isolated and cluster massive-star progenitors (see text), and a kick velocity distribution from Igoshev (2020). Dashed and dotted lines mark the 80% and 20% probabilities, respectively.

In the text

Fig. 6. Probability for a selection of 22 known pulsars within 2 kpc to be found inside their parent environments. Pulsars in this sample have proper motion information vt, which enters as a further constraint on their actual situation (see text). Coloured symbols with black edges indicate the probability for vkick = vt. Boxes indicate the 68% containment interval over possible kick velocities vkick ≥ vt, while whiskers extending above and below the boxes show the 95% containment interval. The horizontal segment inside each box is the median probability.

In the text

	Fig. 7. Probability that Geminga resides in its parent environment for five different assumptions on its age. For each age, a box plot such as those in Fig. 6 is displayed.
In the text