Halo fraction in TeV-bright pulsar wind nebulae

The discovery of extended TeV emission around the Geminga and PSRB0656+14 pulsars, with properties consistent with free particle propagation in the interstellar medium (ISM), has led to the suggestion of “TeV halos” as a separate source class, which is distinct from pulsar wind nebulae. This has sparked considerable discussion on the possible presence of such halos in other systems. In deﬁning halos as regions where the pulsar no longer dominates the dynamics of the interstellar medium, yet where an over-density of relativistic electrons is present, we make an assessment of the current TeV source population associated with energetic pulsars in terms of size and estimated energy density. Based on two alternative estimators, we conclude that a large majority of the known TeV sources have emission originating in the zone that is energetically and dynamically dominated by the pulsar (i.e. the pulsar wind nebula), rather than from a surrounding halo of escaped particles diffusing into the ISM. Furthermore, whilst the number of established halos will surely increase in the future since there is a known large population of older, less energetic pulsars, we ﬁnd that it is unlikely that such halos contribute signiﬁcantly to the total TeV γ -ray luminosity from electrons accelerated in pulsar wind nebulae due to their lower intrinsic surface brightness.


Introduction
The size of the zone within which a pulsar's influence is dominant, that is, the pulsar wind nebula (PWN), is primarily dictated by the injected energy and the external pressure. Within this zone, relativistic particle propagation may be dominated by advection, rather than diffusion. Particle acceleration to >TeV energies is thought to occur at or near the pulsar wind termination shock, located at ∼0.1 pc from the pulsar; however, relativistic particles are present, and radiating, on potentially much larger scales. The usual formalism applied to study the evolution of a PWN is magneto-hydrodynamics (MHD) and in this description the escape of particles cannot be easily included, but it is generally accepted that for young pulsars, such as the Crab, the confinement within the PWN is very effective. At some point in the evolution of the PWN, the mixing of the low density PWN material with the interstellar medium (ISM) or the reverse shock of an associated supernova remnant (SNR) is likely, and the PWN boundary becomes less clear cut and escape becomes possible. Escaping particles encounter the conditions of the ISM (or perhaps modified conditions associated with the earlier SNR, cf. Evoli et al. 2018) and their transport becomes dominated by diffusion, thus potentially forming a detectable halo of escaped particles and/or reaching the Earth to contribute to the measured local electron and positron fluxes. The recent detection of extended TeV emission around the two very nearby low power pulsars Geminga and PSR B0656+14 (Abeysekara et al. 2017a,b) has led to an extensive discussion of TeV halos and their contribution to the general population of TeV emitting PWN (Linden et al. 2017).
We propose a definition of an electron halo (referred to as "halo" hereafter) to be the presence of an over-density of relativistic electrons around a source or acceleration site in a zone in which the source itself does not dominate the dynamics or composition of the interstellar medium. This implies that diffusion is always expected to dominate in halos unless sources exist in regions where there are very large scale flows, such as a galactic wind or outflow. Alternative definitions have been suggested that are related to the emission region size in different wavebands and the dominant transport mechanism (Linden et al. 2017). The size of the emission region only corresponds to a physical boundary in the case of emission from un-cooled particles, for example, in the radio band; the highest energy electrons cool very rapidly and hence never reach any boundary or physical transition of the system. However, in all wavebands, the extent of the observed emission may also vary due to instrumental sensitivity and biases. Whilst a definition based on transport is appealing, in practice, it is extremely difficult to disentangle the transport processes inside the region that is energetically dominated by the pulsar to define a boundary at which diffusion begins to dominate (see e.g. H.E.S.S. Collaboration 2019). We note that halos may not always be symmetric. Electrons follow magnetic field lines in the interstellar turbulence around the source; additionally, in cases where the bulk of emitting electrons is located at distances from the source that are smaller than the coherence length of the turbulence, the halo is expected to appear asymmetric and filamentary .
The potential existence of γ-ray emission in halo regions around energetic pulsars has long been known; their detectability has been predicted to depend on a low diffusion coefficient (Aharonian 1995). The formation of a halo is possible around any source class from which cosmic-ray electrons may escape before cooling. For GeV-emitting electrons, this has long been assumed to be the case for SNR, for example. For the >10 TeV electrons producing TeV γ-rays, radiative lifetimes are much shorter, ∼10 4 (B/10 µG) −2 (E e /10 TeV) −1 yr, and it is less clear if escape before cooling is possible for SNR or for other galactic sources.
After long being discussed as a potential source of the locally measured cosmic-ray electrons (see e.g. Atoyan et al. 1995), pulsar wind nebulae are the only class of sources for which the escape of TeV-emitting electrons is now firmly established (Abeysekara et al. 2017a). There are also indications of missing (apparently escaped) high energy electrons in the case of the Vela system (Hinton et al. 2011); however, alternative explanations, such as rapid cooling due to an enhanced magnetic field, are possible (Bao & Chen 2019). Apparently escaped relativistic electrons are detected in X-rays around the Guitar Nebula (e.g. Johnson & Wang 2010) and the Lighthouse Nebula (e.g. Pavan et al. 2014). The most obvious way to differentiate between escaped and confined radiating particles is to estimate the PWN boundary using multi-wavelength data. This process is often problematic due to the effects of instrumental sensitivity, non-uniform magnetic fields, and particle cooling. For example, it is very often the case that the zone in which imaging of a PWN is possible in the X-ray domain, with instruments such as Chandra or XMM-Newton, is much smaller than the physical size of the PWN as determined in other wavelengths. Within the X-ray domain, the physical PWN size is also often energy dependent, which is interpreted as a signature of the rapid cooling of the highest energy electrons producing the keV synchrotron emission. Indeed, the typical cooling time of electrons emitting photons with characteristic energy hν c is ∼10 3 yr (B/10 µG) −3/2 (hν c /5 keV) −1/2 . In the radio domain, the cooling effect is unimportant, but surface brightness sensitivity is usually only sufficient for young and compact sources.
Here, we consider various estimates for the expected size of the nebulae around pulsars that have been associated to TeV emission, comparing these estimates to the measured sources sizes. We also assess the fraction of the power that is present in sources with and without halos and hence their contribution to the total γ-ray emission of all pulsars within star-forming systems.

Pulsar wind nebula evolution
According to the definition found above, halos may only exist around PWN whose electrons and positrons have started to escape into the surrounding, unperturbed ISM. It is therefore instructive to briefly recall the main stages of the evolution of a PWN. The environment of pulsars changes dramatically over time, firstly as contained within an evolving supernova remnant (SNR), and finally within the general ISM when the "kick" velocity received by the pulsar at birth moves it beyond the decelerated shell of the host SNR. There is considerable literature associated with PWN evolution, including several reviews, see in particular Gaensler & Slane (2006). In general, however, the existing work focuses on X-ray and radio emission, rather than TeV emission, and/or exclusively on the early to middle stages ( 100 kyr) of PWN evolution. Here we briefly consider the physical properties of the region from which TeV emission originates during the lifetime of a pulsar. Figure 1 illustrates three stages in the evolution of a TeVemitting PWN. In chronological order, we depict the following: first, the system at early times t 10 kyr after the supernova in the upper left panel, then intermediate times t ∼ 10−100 kyr in the upper right panel, and, finally, late times t 100 kyr in the lower panel. Hereafter, we refer to these three stages as stage 1, stage 2, and stage 3, respectively. In all three panels of this sketch, the kick velocity that is initially imparted to the pulsar during the supernova explosion is assumed to point towards the left, and the ISM density gradient in which the SNR evolves points "upwards". The areas shaded in grey correspond to the SNR, and the surrounding solid, dashed, or dotted green lines denote the location of its forward shock. The black dots show the location of the pulsar, the PWN is shaded in blue, and the pulsar wind termination shock is represented with the thin solid blue line inside the PWN. The inset in the lower panel corresponds to an enlargement of the innermost regions of the PWN in stage 3. The high-energy electrons and/or positrons that are responsible for the TeV emission from PWN are usually thought to be accelerated at the termination shock. We note that electrons and positrons are often assumed to be accelerated in equal quantities, but this may not be the case. Giacinti & Kirk (2018) found that an individual pulsar may favour either positrons or electrons at the highest energies. Hereafter we use the word "electrons" indiscriminately when referring to electrons and/or positrons, unless stated otherwise. In each panel, with thin magenta lines, we sketched a few illustrative trajectories for some of these 10 TeV electrons. The dotted red lines delineate the typical extent of the 1 TeV γ-ray emission resulting from these electrons (see Fig. 1).

Stage 1: < 10 kyr
At early times (stage 1), the pulsar is still relatively close to its birthplace. The relativistic electron-positron wind that surrounds it, and which is decelerated to non-relativistic speeds at the termination shock, inflates a nebula inside the parent SNR. The ejecta from the supernova expands at supersonic speeds in the surrounding ISM, creating a forward shock (labelled "FS"). Inside the volume delimited by the contact discontinuity ("CD") between the ejecta and the shocked ISM, a reverse shock ("RS") appears as the expansion slows down. At these early times t 10 kyr, the reverse shock has not reached the PWN yet. The highenergy electrons accelerated at the PWN are thought to remain confined inside at this stage, and the TeV γ-rays emitted by these electrons should therefore come from within the nebula. Therefore, TeV emission from PWNe in stage 1 cannot be associated to halos. For the sake of simplicity, we drew the nebula as a spherically symmetric system, but we note that the reality may be more complex, see for example, the Crab Nebula where prominent filaments from the Rayleigh-Taylor instability are present. Eventually, the reverse shock returns inwards, towards the centre of the explosion. As is depicted in Fig. 1, the SNR expands more slowly in the direction where the ISM density is higher, and the reverse shock may then reach the PWN earlier from some directions than in others. At the end of stage 1 and the beginning of stage 2, the reverse shock crushes and disrupts the PWN; against which the PWN forward shock rebounds with the system experiencing several reverberations between the shock and the PWN (Blondin et al. 2001). Upper left panel: early times, t 10 kyr ("stage 1"), when the PWN is contained inside the SNR and before the reverse shock (RS) interacts with it. The SNR forward shock (FS) and contact discontinuity (CD) are plotted with green lines. The electrons that are responsible for the TeV γ-ray emission of the nebula are thought to be confined within the nebula at this stage. Upper right panel: intermediate times, t ∼ 10−100 kyr ("stage 2"), after the PWN is disrupted by the reverse shock, but before the pulsar escapes its SNR. At this stage, TeV γ-ray emitting electrons start to escape from the PWN into the SNR and possibly into the ISM. Lower panel: system at late times, t 100 kyr ("stage 3"), when the pulsar has escaped from its currently fading parent SNR. At this stage, high-energy electrons escape into the surrounding ISM, and may, only then, form a halo under our definition. See the text in Sect. 2 for more details. The key is in the lower left corner. In all three panels, the ISM density gradient is upwards, and the pulsar "kick" velocity is towards the left.
ISM as well as on the direction and velocity of the pulsar. The nebula is disrupted and the pulsar can be strongly off-centre with respect to the PWN. In our sketch, we drew the typical geometry one would expect for a system evolving in a smooth region of the ISM with a steady density gradient perpendicular to the pulsar velocity (see, for example, Figs. 12 and 13 in Slane et al. 2018 for hydrodynamical simulations of Vela X, as well as Fig. 6 in Temim et al. 2015 andFig. 8 in Kolb et al. 2017 for simulations of G 327.1−1.1). At this stage, high-energy electrons start to escape from the PWN and propagate into the surrounding SNR; further escape into the surrounding ISM then becomes possible. A study of electron escape from Vela X is presented in Hinton et al. (2011). In Fig. 1, we show that the extent of the resulting TeV emission may now be greater than that of the nebula. This still does not constitute a true halo, in that the bulk of the TeV emission originates from electrons that do not propagate in the "unperturbed" ISM. Stage 2 may, however, correspond to TeV halos under the definition of Linden et al. (2017), as particle propagation may be diffusive within the expanding region producing non-thermal emission. We note that the sketch in Fig. 1 of Sudoh et al. (2019) corresponds to a standard PWN-SNR system in its stage 2, albeit in the specific case of a homogeneous ISM density and no pulsar motion.

Stage 3: >100 kyr
At late times, typically t 100 kyr (stage 3, see Fig. 1), the pulsar has finally escaped from its parent SNR due to its kick velocity. At this stage, the SNR expands very slowly and fades away. The pulsar propagates in the ISM and forms a bow-shock PWN, with a tail of shocked pulsar wind trailing behind (blue region in Fig. 1). The inset shows an enlargement of the region around the head of the nebula. Observations of X-ray "filaments" around some bow-shock PWNe demonstrate that high-energy electrons can escape into the surrounding ISM. See, for example, Johnson & Wang (2010) and Hui et al. (2012) for the Guitar Nebula, and Pavan et al. (2014Pavan et al. ( , 2016 for the Lighthouse Nebula. Theoretical studies also suggest that some of the high-energy electrons inside a bow-shock PWN can escape into the ISM, even in regions that are not far from the head of the nebula: see, for example, Bucciantini (2018) and Barkov et al. (2019). These particles should then diffuse in the surrounding turbulent interstellar magnetic fields and emit TeV γ-rays in a volume that is substantially larger than that of the PWN.
Only at this stage, and under the condition that escaped relativistic electrons do not dominate the energy density of the ISM, should the TeV emission be considered as a true halo system. In general, only may old PWNe with ages 100 kyr (or, at least, a few tens of kyr) be surrounded by halos of TeV γ-ray emission. The recent HAWC detection of halos at TeV energies around the two stage 3 PWNe of Geminga pulsar and PSR B0656+14 (Abeysekara et al. 2017a) is consistent with this picture. We consider halos at this stage to constitute a clearly phenomenologically distinct source class from the PWNe of stages 1 and 2, whilst acknowledging that stage 2 systems may present ambiguous properties (see Table 1). A113, page 3 of 10   Table 1 provides a summary of the properties of several wellstudied PWN systems in different evolutionary stages, which roughly correspond to the three stages illustrated in Fig. 1. Where the evolutionary stage is denoted as "1b" or "2b" in Table 1, the system may be considered to be currently between stages. This assignment is roughly based on the pulsar characteristic age; however, this may be a poor estimate of the true system age.

Properties of known PWNe
In particular, the morphology of HESS J1825-137 suggests transitory behaviour with evolutionary models favouring an older age for the system (H.E.S.S. Collaboration 2019; Van Etten & Romani 2011).
To explore the relative sizes of the different components in a PWN-SNR system and their evolution, Table 1 compiles multiwavelength measurements of the emission extent. The size of the SNR, R SNR , is given by the radius of a shell, which is often detected from radio information. The size of the central PWN, R PWN , is obtained from the extent of radio emission located immediately around the pulsar. The X-ray size, R X−ray , provided for comparison is a measure of the central PWN size, which traces to the youngest energetic particles in the system. For the synchrotron emission seen in the radio to X-ray, the observed sizes relate to the magnetic field distribution as well as the particle distribution. In contrast the TeV size, R TeV , only depends on the particle distribution, at least in the usual case of close to uniform radiation density for inverse Compton scattering.
At early evolutionary times (stage 1), both X-ray and TeV emission fill the available region for a similar overall size as the PWN remains confined by a shock front. Towards later times (stages 2-3), the discrepancy between X-ray and TeV size becomes more pronounced as high energy electrons, producing synchrotron X-ray emission, cool close to the source; whereas older, lower energy electrons continue to produce significant TeV emission by IC scattering out to large radii due to particle transport and/or proper motion of the pulsar. However, the measured extension of X-ray emission may also be limited due to absorption along the line of sight and instrument sensitivity. Typically, R TeV is therefore somewhat larger than R X−ray for older systems, with R TeV /R X−ray ≈ 5 for HESS J1825-137 and ≈100 for Geminga, whilst for the younger system MSH 15-52, R TeV /R X−ray ≈ 1. We note, however, that here R TeV is given by the one sigma radius of the emission from Abdalla et al. (2018b), whereas R X−ray is given by the total visible extent of X-ray emission.
Lastly, as a further scale comparison, the pulsar kick velocity v and age of the system t are used to obtain an estimated displacement of the pulsar from its birth place. In Table 1, a value of 300 km s −1 is adopted, corresponding to the average of known pulsar kick velocities (Hansen & Phinney 1997). In comparing the emission extent for specific PWN-SNR systems across multiple wavelengths, differences in projection effects arising from variation in the extent along the line of sight have been neglected and they are assumed to be small with respect to the distance to the system.

Energy density estimates
Using two different methods, we now estimate the energy densities ε e in relativistic electrons and positrons around pulsars with established associated TeV emission; see Table 2. Our selection of TeV sources comprises those identified as PWNe by the HESS collaboration in their study of the PWN population (Abdalla et al. 2018b), which has been expanded to include Geminga and PSR B0656+14 (Abeysekara et al. 2017a). Comparing these energy densities with the typical energy density of the ISM, ε ISM , allows one to determine in a first approximation whether the electrons that are responsible for the VHE γ-ray emission occupy the relatively unperturbed ISM (ε e ε ISM ) or if they are still contained in a region that is energetically and dynamically dominated by the pulsar (ε e ε ISM ).  (Manchester et al. 2005) are also shown and are numbered in cases where the association is not unique. We note thatĖ and τ c are the spin down luminosity and characteristic age of the associated pulsar. The sources are sorted by increasingĖ. The estimated distance to the pulsars, which may have large uncertainties, is taken from the ATNF catalogue (Manchester et al. 2005), unless a different reference is provided. The size is the 68% containment radius of the

Estimator using pulsar properties
Measured properties of the associated pulsar can be used to estimate the total power injected into accelerated electrons. The energy density is estimated simply as ε e = E inj /V, where V = 4πR 3 /3 and R is the 68% containment radius for the TeV emission. For Geminga and PSR B0656+14, we calculated, according to the morphology reported in Abeysekara et al. (2017a), the radius at which 68% of the total γ-ray emission is contained. For the remaining sources, the sizes reported in Tables 1 and  3 of Abdalla et al. (2018b) were scaled to correspond to 68% containment and upper limits are kept as reported. The following two cases are considered: (A) E inj =Ė 0 τ 0 , whereĖ 0 is the initial spin down power and τ 0 is the characteristic spin-down timescale; and (B) E inj =Ėτ c , where τ c is the characteristic age of the pulsar andĖ is the present spin-down power. Case A represents the case of injection of essentially the full rotational energy of the pulsar soon after birth. We adopt values for the initial characteristic spin-down timescale corresponding to an extreme value of the birth period P 0 = 10 ms in order to provide an upper limit to the energy content.This assumption leads to inferred energy densities that are larger than 10 eV cm −3 for all the sources considered here. Whilst for older sources, where cooling losses cannot be ignored, this is clearly an overestimate; it is worth noting that a low birth period can easily result in very high energy densities for systems of age < a few 10s of kyrs.
In Case B, we adopt the current spin down power of the pulsar as a measure of power input, again neglecting energy losses, hence: E inj =Ėτ c . The results with this estimator are shown in the two upper panels of Fig. 2 and in the second column from the right in Table 2. The estimates in this case are typically a factor of 10−100 lower than those of Case A. This case is pessimistic in the sense that it ignores the early evolution of the pulsar where significant energy injection may have taken place at higherĖ; however, it is optimistic in that it ignores losses and energy input in other channels (e.g. expansion of the nebula, magnetic fields). We consider Case B to be a better guide, in general, particularly for evolved sources as losses are unlikely to be negligible for TeV γ-ray emitting particles in the early life of the PWN when magnetic fields are high.
Clearly, more sophisticated treatments are possible, but all rely on models of the evolution of pulsars and their nebulae, which are not currently well constrained by the available data. Instead, we turned to existing measurements at TeV energies to empirically establish the current day electron population in (and/or around) these PWN.

Estimator using TeV γ-ray luminosity
The second method adopted consists of estimating ε e from the measured VHE γ-ray emission of these objects. For Geminga and PSR B0656+14, we used the luminosities and spectral indices measured by the HAWC Collaboration over the energy range E γ = 8−40 TeV (Abeysekara et al. 2017a). For all the other PWNe considered here, we used the measurements of their luminosities and spectral indices over the energy range E γ = 1−10 TeV, as reported in Abdalla et al. (2018b). Since the electrons emitting at such high energies only provide a fraction of the total energy density ε e , we needed to make a phenomenological assumption on the shape of the electron spectrum outside the energy range where it is constrained by these observations. To do so, we note that the electron spectra inferred for those PWNe that are observed at multiple wavelengths are roughly compatible with broken power-law spectra, with a low-energy break at E low ∼ 100 GeV, and a high-energy break at E high ∼ 1−10 TeV.
At energies E low E e E high , the spectrum is roughly compatible with an ∼ E −2 e spectrum, and it softens at E e E high . At E e E low , the spectrum is, on the contrary, harder than E −2 e , implying that the electrons at these low energies hardly contribute to the total energy density in relativistic electrons; see, for example, Meyer et al. (2010) for the Crab Nebula. We set E low = 100 GeV and neglected the contribution of all electrons with E e < E low in our calculations of ε e . In line with the theoretical expectation for particle acceleration at an ultra-relativistic shock with isotropic scattering, we assume that the electron spectrum between E low and E high follows a power-law ∝ E −2.2 e (e.g. Achterberg et al. 2001). In order to test the dependence of our results on the a priori unknown value of E high , we carried out our calculations for three different values: E high = 1, 3, and 10 TeV. For each of these values of E high , we tried two different slopes for the electron spectrum at E e > E high . First, we took a spectrum ∝E −3.2 e , which is consistent with the expectation for a cooled E −2.2 e spectrum 1 . Second, we took a spectrum ∝E −Γ e , where Γ was adjusted individually for each source such that the spectral index of its calculated γ-ray emission is consistent with the spectral index measured in the relevant energy range.
For each of the six electron spectra described above, we calculated, for each source, the γ-ray emission from inverse Compton scattering of these electrons on background radiation fields (CMB, UV, optical, and IR), using the Eq. (2.48) of Blumenthal & Gould (1970), which takes Klein-Nishina effects into account. The spectra of the background UV, optical, and IR photons were calculated using the model of Popescu et al. (2017) for the interstellar radiation fields. These fields were determined individually for each PWN, depending on its location in the Milky Way. Assuming the same source radii R as in the first method, we then inferred the electron densities ε e that are required to fit the observed γ-ray luminosities in the energy ranges E γ = 1−10 TeV for the H.E.S.S. paper sources, and E γ = 8−40 TeV for the two HAWC sources. We note that since we only calculated typical average energy densities ε e , our calculations are formally equivalent to assuming that the electrons are uniformly distributed within the volumes V = 4πR 3 /3 of their sources. When we used the three electron spectra that are ∝E −3.2 at E e > E high , we only fitted the (total) measured γ-ray luminosity. In contrast, for the three electron spectra with adjustable slopes at E e > E high , we simultaneously fitted the measured γ-ray spectral index. For each source, we calculated six values of ε e -one per electron spectrum. Hereafter, we provide the average of these values and indicate the spread with asymmetric errors.
The results are given in the last column of Table 2, and they are shown in the lower row of Fig. 2 as a function of the spindown power (lower left) or characteristic age (lower right) of the associated pulsar. As in the upper row of Fig. 2, the colour of each symbol refers to the source radius R (see the colour bar on the right-hand side), and the area shaded in grey corresponds to the region where ε e < 0.1 eV cm −3 , that is, where ε e ε ISM . Since the values of ε e calculated with this second method are derived from measurements, they are likely a better indicator of the actual electron energy densities at the sources Top panels: energy density calculated as ε e =Ėτ c /V, described in Sect. 3.1. Bottom panels: energy density calculated using the electron spectrum derived from TeV γ-ray measurements, described in Sect. 3.2. The shaded regions correspond to an energy density lower than that of the ISM under any reasonable assumption. Systems for which measured properties are taken from the HGPS are indicated by circles, whereas stars indicate the systems where HAWC data is used (see also Table 2). than those derived with the estimator of Sect. 3.1. Moderately large error-bars are present on the results due to the uncertainties on the shape of the electron spectrum. One can see that almost all of the objects considered here lie above or partly outside the grey band, except for the two HAWC sources (Geminga and PSR B0656+14).

Current halo fraction in TeV-bright PWNe
Figure 2 demonstrates that the VHE γ-ray emission from most TeV-bright PWNe is due to electrons that are contained in a region which is influenced energetically or dynamically by the pulsar (i.e. the PWN itself). Therefore, these sources do not constitute proper halos in which the TeV emission arises from a low energy density zone around the PWN. It is clear from Fig. 2 that the only really unambiguous halo cases in the sample of established TeV-emitting PWN are those powered by Geminga and PSR B0656+14, as established by Abeysekara et al. (2017a). These two HAWC sources stand out as clear halo cases: The electrons responsible for their TeV γ-ray emission propagate in a region where their contribution to the local total energy density is negligible (ε e < 0.1 eV cm −3 ). We note that there is a strong correlation (coefficient 0.98) between the two different energy density estimates shown in Fig. 2. This indicates that the γ-ray properties of PWN are typically more closely related to the current spin-down power, rather than the early lifetime of their pulsars. Whilst for the majority of systems the estimator using γ-ray luminosity results in consistent or lower energy densities than the Case B estimates of Sect. 3.1, for a small number of PWN, the energy density is larger. This is likely due to significant energy injected during the early phase of pulsar evolution.
The threshold in ε e below which a source can be classified as a halo is, of course, neither sharp, nor universal. Several A113, page 7 of 10 of the sources are close to our indicative grey band. Some of these sources might be reaching the stage of their life when they start to transition towards a halo, although we caution that the presence of incorrect pulsar associations in our sample is not excluded. Several objects are ambiguous cases in which more careful study of morphology, energetics, and environment would be needed for a clear classification. The ambiguous cases have characteristic ages that are consistent with stage 2 of Fig. 1, and it seems plausible that some of these objects are hybrids, exhibiting halo as well as PWN components. HESS J1825−137 is a candidate for this mixed situation. The energy-dependence of the morphology strongly suggests advection dominated transport within the PWN, but the fringes of the emission extend to very large distances and may indicate the presence of unconfined TeV-emitting particles (H.E.S.S. Collaboration 2019). The emission profile measured with HESS exhibits a change in slope on the southern side of the nebula, which may hint at a boundary between a high electron energy density core (the PWN) and a surrounding diffusive halo. This effect is more pronounced in linear strips and averages out over the azimuthal angle around the pulsar (see Figs. 4 and 5 of H.E.S.S. Collaboration 2019). Vela X is a well-known example where the TeV IC emission is concentrated within a central region; it is more extensive than the X-ray synchrotron emission, yet well-contained within the surrounding radio SNR such that we consider this to be a system entering stage 2 (Hinton et al. 2011;Bao & Chen 2019).
In our sample of TeV-bright systems with firm pulsar associations listed in Table 2, we note that the Crab and Geminga, both frequently used as canonical examples of the class, have the highest and lowestĖ, respectively. This suggests that both are extreme systems rather than typical of the population.

Total halo fraction in Galactic PWNe
Whilst the fraction of halos in current TeV catalogues is low, these catalogues are clearly biased towards compact, that is, PWN-like, rather than halo-like, objects, given the reduction in sensitivity to lower surface brightness emission at a fixed flux for resolved objects. In addition, the highest power objects are young (Table 2) and would not be expected to form halos yet; only close-by (old, low-power) pulsars are likely to exhibit detectable halos. However, the older systems which do exhibit halos are much more numerous, and it is important to consider the entire pulsar population to assess what fraction of the integrated TeV γ-ray emission is likely to occur within, and outside of, PWN. Figure 3 shows the cumulative fraction of the energetic output of Galactic ATNF pulsars; this excludes millisecond pulsars 2 (Manchester et al. 2005) as a function of their characteristic age. The energy fraction is extracted from the beaming angle corrected (Tauris & Manchester 1998) pulsarĖ, which is shown assuming the current pulsar properties. In order to better assess the contribution of different aged pulsars to the total energy in Galactic electrons, the total energy output of the pulsars was integrated over timescales of 10 4 and 10 5 yr (or the pulsar age if shorter), corresponding to the radiative lifetime in the ISM of around 10 and 1 TeV electrons, respectively. The energy output of the pulsars were assumed to evolve identically as 1 + t τ 0 −α , with a τ 0 set to 10 3 yr and α = (n + 1)/(n − 1) assumed to be equal to 2 as for a pulsar braking index n = 3. Assuming that the power injected in to relativistic electrons evolves in the same manner, a simple evolutionary model using the GAMERA package  (Manchester et al. 2005) created using the instantaneous energy output into electrons E e as well as the energy output integrated up to 10 4 and 10 5 yr electron lifetimes.
For comparison purposes, we show the cumulative γ-ray luminosity at 1 TeV from a generic PWN model as a function of system age, both for a constant magnetic field and a magnetic field evolving in time. The cumulative fraction of pulsars indicates that the energy output into electrons and γ-ray luminosity is dominated by a small fraction of young, energetic pulsars. (Hahn 2016) was constructed to predict the expected evolution of the PWN γ-ray luminosity at 1 TeV, both in a constant magnetic field environment and for a magnetic field evolving as (t/τ 0 ) −0.5 . The cumulative energy fraction in 1 TeV γ-rays from this model is also shown in Fig. 3. As expected, these curves are approximately consistent with the cumulative energy fraction in the known pulsar population integrated over an electron lifetime of 10 4 yr, corresponding to 10 TeV electrons, which are responsible for 1 TeV γ-ray emission. Figure 3 shows that the contribution of 10 5 yr old pulsars to the total TeV γ-ray luminosity of PWNe (black curves) is small, in general, and it depends on the cooling history of electrons injected early in the life of pulsars, when most of the rotational energy is lost. The γ-ray luminosity of older halo systems is low, as the majority of the power is injected into the PWN in the early phase of the pulsars evolution, with particle cooling expected prior to escape into the ISM. By contrast, stage 2 systems may account for up to a third of the total γ-ray luminosity. The overall efficiency for pulsars as cosmic-ray electron sources depends on the interplay of the electron lifetime and effectiveness of particle escape.

Electron diffusion and escape
The electron diffusion coefficient measured by HAWC around Geminga and PSR B0656+14 is close to the Bohm limit and is two orders of magnitude smaller than the one inferred from the boron-to-carbon ratio (Abeysekara et al. 2017a). This naturally raises the question of whether the electrons responsible for these two halos probe the unperturbed interstellar turbulence or the turbulence generated by the cosmic rays that have escaped from these (or nearby) sources. Cosmic-ray self-confinement around hadronic sources has been studied by Ptuskin Evoli et al. (2018) point out that self-confinement could also occur around sources of leptons, suggesting that this may explain the small diffusion coefficient found by HAWC around Geminga and PSR B0656+14. Our finding in Sect. 3.2 that ε e ε ISM around these sources casts doubts on the ability of their electrons to modify the properties of the interstellar turbulence substantially. It cannot be excluded, though, that larger electron currents during an earlier phase, or cosmic rays from a nearby hadronic source, such as the parent SNRs of these pulsars, could have modified the turbulence in the halo regions around these two PWNe. However, it is worth noting that there is currently no need to invoke cosmic-ray-driven instabilities to explain the data. In particular, the value of the cosmic-ray diffusion coefficient measured by HAWC is compatible with theoretical expectations for isotropic turbulence with strengths and coherence lengths that are in the relevant ranges for interstellar turbulence (López-Coto & Giacinti 2018). If cosmic rays diffuse substantially faster in the Galactic halo than in the local ISM, tensions with the boron-tocarbon ratio might be avoided. Indeed, the boron-to-carbon ratio provides an estimate of the amount of matter that cosmic rays traverse before reaching Earth. Assuming homogeneous, isotropic diffusion in the entire Galaxy, one finds a typical diffusion coefficient value of ∼10 28 cm 2 s −1 for GeV cosmic rays. However, this does not rule out the existence of regions with substantially slower diffusion; for instance, an example is the Galactic disc, so long as the grammage accumulated by cosmic rays in these regions does not overshoot that from the boron-to-carbon ratio. Cosmic rays would diffuse faster in the Galactic halo if it contains an out-of-plane, ordered magnetic field and has lower turbulence levels than the disc (e.g. . H.E.S.S. observations of electrons at Earth with >10 TeV energies place stronger constraints on the possibility of a low diffusion coefficient in our local environment than the boron-tocarbon ratio does. As demonstrated in Hooper & Linden (2018), cosmic rays must diffuse substantially faster through the bulk of the Galactic disc than around Geminga for local known pulsars to explain the detection of >10 TeV electrons at Earth. Alternatively, the scenario of a low diffusion coefficient in our local ISM would remain viable if a nearby undetected pulsar contributes to the high-energy end of the all-electron spectrum measured at Earth ). In the future, observations of halos around pulsars in different γ-ray energy ranges should help to constrain the spectrum and nature of the turbulence probed by the emitting electrons. Extended GeV emission has recently been measured with the Fermi satellite around Geminga (Di Mauro et al. 2019).
Here, we discuss the halo fraction in TeV-bright PWNe; however, the formation of halo-like emission at TeV energies due to electron escape is a general phenomenon expected from cosmic accelerators. One possible example is the TeV emission that is detected beyond the shell in the SNR RX J1713.7-3946, which is hypothesised to be either due to a shock precursor or the escape of energetic particles from the shock region (H.E.S.S. Collaboration 2018). We note, however, that in this system and other candidate cases for leptonic halos, there is often ambiguity with a hadronic emission scenario.

Future prospects
Future detections of diffusive halos would considerably help to constrain this picture; the prospects for this were already discussed in Linden et al. (2017). New halos may have already been detected: the HAWC Collaboration recently announced the detection of two more candidate halos -HAWC J0543+233 around PSR B0540+23 (Ė = 4.1 ×10 34 erg s −1 , distance = 1.56 kpc, τ c = 253 kyr, found in a search for extended sources with disc radius 0.5 • (Riviere et al. 2017)) and HAWC J0635+070 around PSR J0633+0632 (Ė = 1.2 × 10 35 erg s −1 , distance = 1.35 kpc, τ c = 59 kyr, with a Gaussian 1σ extent of 0.65 • ± 0.18 • (Brisbois et al. 2018)). The pulsars that may be powering these sources are in an advanced state of their evolution, and both fulfil the condition of being older than a few tens of kyr. A rough estimate of their energy densities can be made using the method of Sect. 3.1, which yields values of 0.6 eV cm −3 (using the 0.5 • disc radius in lieu of the 68% containment radius) and 0.09 eV cm −3 (for a 23 pc 68% containment radius) for HAWC J0543+233 and HAWC J0635+070, respectively. Whilst these values are certainly consistent with those of Geminga and PSR B0656+14, using the estimator based on pulsar properties, it is difficult to provide a more accurate placement in Fig. 2 using the method of Sect. 3.2 due to the lack of more detailed information about their size and spectral properties. Because of the spin-down power and characteristic age of their central source, they should be located inside the grey band, corresponding to ε e ε ISM , and remain clear candidates to be TeV halos; although, HAWC J0635+070 may be an ambiguous case. Despite their low luminosity, the total number of known halo sources may significantly increase with improvements in sensitivity of current and future γ-ray experiments, as a large population of older, less powerful pulsars is known. Linden et al. (2017) propose that TeV-bright sources coincident with PWNe are a distinguishable halo phenomenon; most of the sources included in that study are consistent with the classical picture of a middle-age PWN (stage 2 of Fig. 2) with effective confinement, rather than an escaped population of electrons propagating in the ISM. Examples where the X-ray emission region boundary does not correspond to a physical discontinuity include the Crab nebula, MSH 15-52, and G327.1-1.1, as listed in Table 1, as well as G359.23-0.82 ("the Mouse" nebula), as shown in Fig. 6 of Gaensler & Slane (2006). The tail of the Mouse nebula is clearly longer in radio than in X-rays. From the observational point of view, it must be noted that the X-ray and TeV γ-ray size of these sources are expected to differ even in fully confined systems, for reasons discussed in Sect. 2. Observational γ-ray data may be used to determine whether individual sources are better classified as PWNe, halos, or hybrid systems. Both the energy density and particle transport properties within these systems can be derived by model fitting to the spectral energy distributions and angular emission profiles, respectively, thus enabling classification according to the evolutionary stage. We note, however, that neither property can be determined in a truly model-independent manner.

Conclusions
From this study, we conclude that the fraction of TeV-bright PWNe exhibiting halo properties is low; however, the fraction of PWNe that produce electron halos is presumably large and this is a generic feature of all pulsars with ages 10 5 yr. In the future, CTA (Cherenkov Telescope Array Consortium et al. 2019), SGSO/SWGO (Albert et al. 2019), and LHAASO (Bai et al. 2019) will detect new candidates and the number of halos will increase further. Nevertheless, the total halo fraction of TeV-bright γ-ray PWN is likely to remain small due to their comparatively low γ-ray luminosity (Fig. 3) and large angular size, limiting future detections towards local objects. A study of the prospects for future detections of halos by SGSO/SWGO is presented in the Sect. 4.1 of Albert et al. (2019).