© ESO 2018

1 Introduction

Pulsar wind nebulae (PWNe) are clouds of magnetised electron-positron plasma that can span many parsecs and are observed via their synchrotron or inverse Compton (IC) radiation (see Gaensler & Slane 2006, for a comprehensive review on the subject). They are created inside supernova remnants (SNRs) by the energetic outflow (“wind”) of a pulsar, which is a swiftly rotating neutron star that is the compact leftover of the supernova explosion. The pulsar wind runs into the supernova ejecta and develops a standing shock wave beyond which the PWN builds up as an expanding bubble of diffuse plasma. outflow is decreasing steadily. Therefore, most of the observed PWNe are associated

It is instructive to consider the energetics of a typical PWN system. A pulsar releases a total amount of energy of up to 10⁴⁹ –10⁵⁰ erg over its lifetime, but only ≲10 % of this energy is emitted as pulsed electromagnetic radiation (Abdo et al. 2013). Most of the pulsar outflow consists of high-energy particles and magnetic fields that feed into the growing PWN plasma. This plasma is dynamically inferior to the ~ 10⁵¹ erg carried away by the supernova blast wave around it. A good portion of the PWN energy is radiated off, predominantly through synchrotron emission in the first few thousand years, which can be observed in the X-ray and radio bands. Besides that, a few percent of the PWN energy are converted to IC radiation in the TeV regime. In Gould (1965), but also in later works (De Jager et al. 1995; Du Plessis et al. 1995; Aharonian & Atoyan 1995), it was already suggested that this could allow for the detection of TeV emission. And even though the IC photons are an energetically subdominant emission component, they carry important information that the synchrotron emission, albeit much higher in flux and energy transport, does not give access to; they emerge predominantly from homogeneous, time-constant CMB and IR photon seed fields, and therefore trace the electron plasma independent of the time- and space-varying magnetic fields. In Aharonian et al. (1997), it was suggested that the TeV nebulae could be much larger nebulae than those observed in the radio or X-ray regimes. So in general, the IC image gives a more accurate and complete picture of the electron population than the synchrotron photons.

Indeed, since the TeV detection of the Crab PWN in 1989 with the Whipple telescope (Weekes et al. 1989), tens of Galactic sources havemeanwhile been associated with TeV pulsar wind nebulae. Most of these objects are situated in the inner Galaxy; many were therefore discovered and extensively investigated from the southern hemisphere using the H.E.S.S. Imaging Atmospheric Cherenkov Telescope (IACT) array (e.g. Aharonian et al. 2006d), which can observe the inner Milky Way at low zenith angles and high sensitivity. The northern IACT systems MAGIC (e.g. Aleksić et al. 2014) and VERITAS (e.g. Aliu et al. 2013), and arrays of air shower detectors such as MILAGRO (Abdo et al. 2009), have also observed PWNe and contributed very valuable case studies, mostly of systems evolving in the less dense outer Milky Way regions. Also HAWC shows promising potential to contribute new data soon (Abeysekara et al. 2015) but has not provided a major data release yet. In the 1–10 TeV regime, IACTs generally have a better angular resolution and sensitivity than air shower arrays, even though their fields of view (FOV) are limited to one or few objects, and their duty cycle is restricted to dark, cloudless nights.

A systematic search with the Fermi Large Area Telescope for GeV pulsar wind nebulae in the vicinity of TeV-detected sources (Acero et al. 2013) yielded 5 firmly identified high-energy gamma-ray PWNe and 11 further candidates. The PWN detections were also often complemented by multi-wavelength observations in the X-ray or radio bands (see e.g. Kargaltsev et al. 2013).

In this paper, we proceed along the lines of previous work that aimed at a uniform analysis of the whole population of TeV pulsar wind nebulae, such as Carrigan (2007), Carrigan et al. (2008), Marandon (2010), and Mayer (2010a). To do so, we take advantage of the newly released TeV source catalogue (H.E.S.S. Collaboration 2018), which is based on the nine-year H.E.S.S. Galactic Plane Survey (HGPS). It provides a uniform analysis of source sizes, positions, and spectra based on data taken during nearly 3000 h of observations. It covers the Galactic plane at longitudes ℓ = 250° to 65° and latitudes |b|≲3.5°. We undertake a census of all the firmly identified PWNe detected with H.E.S.S. and other IACTs, and for the first time complement this sample with HGPS flux upper limits for all covered pulsar locations without a corresponding TeV detection. This allows for a less biased judgement of the whole population. We compare the common properties and trends of this population to those found in the numerous efforts to theoretically describe the nature of pulsar wind nebulae.

2 Observational data

2.1 HGPS and ATNF catalogues as data sources

We use two different sets of astronomical tables: the H.E.S.S. Galactic Plane Survey¹ (HGPS; H.E.S.S. Collaboration 2018) and the ATNF pulsar catalogue² (Manchester et al. 2005, version 1.54). For most purposes in this paper, the HGPS source catalogue and the full ATNF listing are used. Only the TeV-PSR spatial correlation study in Sect. 3.1 makes use of less biased listings, namely the HGPS components list (HGPSC) and Parkes Multibeam Pulsar Survey (PMPS; Manchester et al. 2001; Lorimer et al. 2006, and references therein), which is a subset³ of the ATNF pulsar catalogue. The HGPSC components list is an unbiased representation of the TeV objects in terms of Gaussian components, which does not invoke a priori knowledge of source associations or other prejudiced assumptions.

For the pulsar distances, we choose the distance estimates of Cordes & Lazio (2002) provided by the ATNF team. Their uncertainty, however, is not very well defined and can be as large as a factor of 2. For the few cases in which pulsar distance estimations were added or replaced from references other than the ATNF pulsar catalogue, these values are listed in Table 2.

2.2 Firmly identified TeV pulsar wind nebulae

To determine which of the known TeV sources should be considered as firmly identified PWNe, we use the identification criteria discussed in the HGPS paper and take as a starting point the list of all 12 identified PWNe and the 8 identified composite SNRs (H.E.S.S. Collaboration 2018, Table 3). Most PWNe in the HGPS are identified by positional and/or morphological coincidence with a PWN identified in other wavelengths, or by their specific (mostly energy dependent) TeV morphology. Our selection for this paper also requires that the corresponding pulsar has been detected and timed; if this is not the case, the properties of the source cannot be put into the physics context of this study, despite its identified PWN nature. This excludes the PWNe in SNRs G327.1−1.1 and G15.4+0.1, and the identified composite SNRs CTB 37A and W41 (see H.E.S.S. Collaboration 2018, Table 3 and references therein). In composite SNRs, the PWN component is mostly believed to outshine the potential contribution from the SNR shell in TeV gamma-rays, and we assume here that this is the case for TeV sources identified as composite SNRs with the exception of HESS J1640−465. For this object, detailed observations with H.E.S.S. suggest that a significant part of the TeV emission may originate from the SNR shell (Abramowski et al. 2014). Therefore, we exclude HESS J1640−465 from firm identification and consider it a PWN candidate. The sample we arrive at is listed in Table 1.

In addition to the firmly identified objects found in the HGPS, we include five HGPS-external PWNe, among them G0.9+0.1, which is inside the plane scan, but was not re-analysed with the HGPS pipeline. These PWNe are displayed using distinct symbols in the figures throughout this work. This latter group, listed in Table 3, is based both on dedicated H.E.S.S. observations outside of the scope of the HGPS and on data from other IACTs.

We do not include detections that are only reported from (direct) air shower detectors, such as MILAGRO, HAWC, or ARGO-YBJ, because their angular and spectral uncertainties are much higher, making the source resolution and pulsar association more difficult and the spectral statements more uncertain.

2.3 Dataextracted from the HGPS

The quantities taken from the HGPS catalogue are source position, extension, integral flux >1 TeV, and spectral index Γ from the power-law fit of the differential photon flux ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$. The extension measure σ is given as the standard deviation of a circular Gaussian function. Extension upper limits were used as provided in the catalogue, namely in cases where the extension is not more than two standard deviations larger than the systematic minimum extension of 0.03°.

Offsets between pulsar and PWN centroid position were calculated and, where necessary, converted to 2σ limits following a similar prescription, namely in the cases where the offset was less than 3σ above a systematic minimum of 0.0056°, which is a typical value for the systematic positional uncertainty of H.E.S.S.

The integral photon flux I_{>1 TeV} and index Γ is converted to a luminosity between 1 and 10 TeV using ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$(1)

where d is the source distance and the integral flux I_{>1 TeV} is taken from the Flux_Map column of the catalogue, which is recommended there as the most reliable estimate of the integral flux. The errors, propagated from the index errors σ_Γ and integral flux errors σ_I, are $${\left(\frac{{\sigma }_L}{\text{lumi}}\right)}^2\text{ \&}={\left[\frac{{\sigma }_I}{I_{>1\text{tev}}}\right]}^2\text{nn\&}+{\left[\left(\frac{1}{(\Gamma -1)(\Gamma -2)}+\frac{\ln 10}{1-{10}^{\Gamma -2}}\right){\sigma }_{\Gamma }\right]}^2.$$(2)

The errors on flux and index are assumed to be independent because the reference energy of 1 TeV is typically very close tothe mean pivot energy of the fits. The errors on the distance estimation of pulsars are not available consistently and are likely not Gaussian in most cases, so they are not treated here and remain a systematic uncertainty. For uniformity, the power-law integration is also used in the few cases where a high-energy cut-off is found to be significant, as the cut-off has very little influence on the integral⁴ .

We also extract flux upper limits from the sky maps of the HGPS data release. The 95 % confidence level limit I>1 TeV on the flux is converted as above, assuming a spectral index of Γ = 2.3, which is the typical TeV index also used in several pipeline analysis steps of the HGPS analysis (H.E.S.S. Collaboration 2018). The flux limits are available for integration radii of 0.1° , 0.2° , and 0.4° ; the latter of which is only available internally and will not be part of the public HGPS data release. For pulsars that qualify for an upper limit, we use the baseline model (Appendix A) to estimate the PWN extension. Assuming 1000 km s^-1 for the offset speed (see Sect. 5.2.2), a required flux limit radius Rlim = R_PWN + Roff is derived and a corresponding angular extent θ_pred as seen from Earth is calculated. If this extension is below 0.4°, the value is rounded up to the next available correlation radius and a flux limit is looked up in the respective limit map. In the case of 0.4° < θ_pred < 0.6°, we assume that the source could have been detected, and calculate a limit from the 0.4° map, scaling it up by ${({\theta }_{\text{pred}}/0.4\text{dg})}^2$ to account for the uncontained part of the PWN. If θ_pred > 0.6°, no limit is calculated since one cannot exclude that a potential weak and undetected PWN emission would have been confused with background in the background subtraction of the HGPS pipeline.

2.4 Caveats of the HGPS

The HGPS data contain unbiased observations, a priori targeted observations, and re-observations of hotspots. It is therefore impossible to raise truly objective and statistically robust statements on chance coincidence detections of TeV objects near energetic pulsars. A way to unbias the data would be to remove all deep and targeted observations from the catalogue construction pipeline, which would obviously discard very interesting parts of the data set and lead to a different catalogue content. We refrain from this exercise here, trying to make use of the richness that is present in the full data set and catalogue.

A uniform source analysis, as provided in the HGPS and exploited here, has many advantages with regard to a population study. The fluxes and extensions are determined with one software version, data quality cut, analysis algorithm, and event selection cut set, leading to values that are comparable and consistently defined among all sources. The disadvantage of uniformity is that it comes with a lack of adjustment. Customised data quality cuts can allow for the detection of weaker sources or for lower systematic uncertainties for very strong sources. This is deliberately not done here.

Besides this, the energy threshold and sensitivity of Cherenkov telescopes vary with the zenith angle of observation, and therefore with the declination of a given sky region. The IACT data thus are intrinsically not completely uniform across different sky regions.

3 Correlation of TeV sources and pulsars

The total energy output of a pulsar at a given time is characterised by its spin-down power Ė , which can be observationally determined from its period P and period derivative Ṗ, assuming a neutron star moment of inertia of I = 10⁴⁵ g cm² (see also Appendix B for the basic formulae of pulsar evolution). Pulsars deploy most of their spin-down energy within few tens of kiloyears. The pulsar wind nebulae thereby created are loss-dominated ever after that period, when the electrons are diffused and lose their energy through radiative or adiabatic cooling with cooling times of $\mathcal{O}(10\text{kyr})$ (see Sect. A.3). Therefore, the natural expectation for a bright PWN is that it has to have an accordingly young ($\mathcal{O}($<10 kyr)) and still energetic pulsar nearby.

Observationally, this is supported by the fact that most TeV pulsar wind nebulae (and sources in general) are found at Galactic latitudes <0.5°; if pulsars were to grow TeV nebulae in their late stage of evolution, then TeV sources should also be more numerous at higher latitudes, where many old pulsars drift off to.

3.1 Spatial correlation

A way to find general support for the association of energetic pulsars and TeV sources was explored by Carrigan et al. (2008), where the whole HGPS sky map of that time was used along with the PMPS pulsar catalogue to evaluate a detection fraction Ndetected∕Npulsars for pulsars in different bands in Ė∕d2.

To investigate whether this spatial correlation is still manifest in the data, Fig. 1 shows the distribution of angular distances between all pulsars of a given range in Ė∕d2 and all “Gaussian components” listed in the unbiased HGPSC component list⁵ . The shaded band shows the expectation derived from simulated pulsar samples. It is derived for the same band of Ė∕d2, calculating 30 000 randomisations of the PMPS pulsar sample. The observed Galactic latitude and longitude distributions of the pulsars are preserved in the reshuffling. A significant correlation beyond chance coincidences is found for pulsars with Ė∕d2 > 10³⁴ erg s-1 kpc-2 and is absent for less energetic pulsars. An estimate for the number of chance coincidences for a cut of 0.5° yields a value of 9.7, while 35 HGPSC components are actually found. Using the full ATNF catalogue instead of PMPS and the HGPS source catalogue instead of the components list, the study is more similar to the source selection we do in the following, but involves statistically less unbiased samples. The estimated number of chance coincidences derived in this case is 11.5.

Fig. 1 Histograms of spatial separation between PMPS pulsars and TeV source components from the HGPSC list. In the high-Ė pulsar sample (left), a clear correlation is seen as a peak at small squared angular distances, whereas the low-Ė associations,if present, are not significant beyond the expected rate of chance coincidences (right). The angular separation cut of θ < 0.5° applied in the preselection of PWN candidates (Sect. 3.2) is indicated by a dashed vertical line in the left panel.

3.2 Pulsar wind nebulae preselection candidates and flux limits

The strategy employed to select and evaluate unconfirmed PWN candidates in this paper is a two-step procedure: First, a loose preselection of candidates has been carried out. Secondly, these candidates are distinctly marked in the various observables correlation plots of Sect. 5, leading to a subsequent judgement on their physical plausibility to be a PWN in the post-selection of Sect. 6.

The criteria we impose for the preselection are that a pulsar should be more energetic than Ė ∕d2 = 10³⁴ erg s-1 kpc-2 and have an angular separation θ from an HGPS source of less than 0.5°. We also require a characteristic age τc < 10⁷ yr to prevent millisecond pulsars, which are different concerning their nature and physics of emission, from entering the PWN candidate sample⁶ . While these criteria are arbitrary to some extent, we note that, as a preselection, they were chosen to be relatively loose and amply include all firmly identified PWNe.

Energetic pulsars that do not have an HGPS source nearby or that coincide with an HGPS source that is already firmly associated to another pulsar are selected for the calculation of a flux upper limit. In the latter case, the flux of the established source is not subtracted, since one cannot isolate one from the other and the conservative flux limit is therefore on top of the emission of the main source. In the limit calculation step, we include all pulsars with Ė > 10³⁵ erg s-1, independentof their distance. For very old and extended objects, a large distance can even be favourable because only then can their full supposed extent be covered within the H.E.S.S. FOV, leading to a meaningful flux limit (see also Sect. 2.3).

For the same reason as in the selection of firmly identified PWNe, we deliberately choose not to treat pulsar systems in which the pulsar is not clearly identified in terms of period, derivative (presumably because the pulsar beam does not intersect Earth), and distance. We require a known pulsar distance so as to be able to quantify TeV properties, such as luminosity and extension, and compare them with the firmly identified population. But we should note that this implies that we cannot consider among PWN candidates the TeV sources coincident with PSR J1459−6053, PSR J1813−1246 and PSR J1826−1256 (see H.E.S.S. Collaboration 2018), which are pulsars that are detected in high-energy gamma-rays but not in the radio domain.

As a caveat of our cut in Ė∕d2, we note that potential ancient nebulae from very old pulsars cannot make it into our selection and are not be considered in this work (except for being included in terms of a flux limit). Figure 1 (right) shows that the TeV detection of such ancient nebulae has to be treated as hypothetical, judging from the global catalogue point of view we adopt in this paper.

The result of the preselection is that besides the 14 firmly identified PWNe we consider here, 18 additional PWN candidates pass the criteria; two of these additional candidates have two pulsars they could be associated with and four pulsarshave two possible TeV counterparts. The 5 HGPS-external PWNe also match the criteria. We exclude the γ-ray binary PSR B1259−63 here. While the TeV source is believed to contain the wind nebula of its pulsar, the TeV emission is clearly impacted by the binary nature of the object and therefore out of the scope of this paper. Also, the obvious TeV shells that were omitted from the standard HGPS pipeline are excluded here, although coincident pulsars are allowed to be included in the limits listing if they qualify.

Among the pulsars without a matching detected TeV source, 65 with Ė > 10³⁵ erg s-1 are selected for the limit calculation; however the assumed PWN extension and offset are small enough to calculate a flux limit with the HGPS maps for only 22 of those. Of these limits, 3 appear to be on top of significant emission for various reasons: PSR J1837−0604 coincides with the PWN HESS J1837−069. The limit of PSR J1815−1738 is integrated over 0.4° and therefore contains parts of the emission of HESS J1813−178. PSR J1841−0524 is situated within the very large HESS J1841−055, possibly consisting of multiple sources; the Ė∕d2 of this object is too low for it to qualify as a candidate.

The pulsars selected as firm PWNe from the HGPS catalogue, as external PWNe, candidate PWNe, and for flux limits are listed in Tables 1, 3–5, respectively. They are shown in the Ė – τc and Ṗ –P planes in Fig. 2. The plots also show ATNF pulsars without detected TeV wind nebula for comparison and highlight some prominent or special objects with labels. These are labeled throughout the paper for orientation.

As expected, the preselection candidates are young, but on average somewhat older than the already established PWNe. This is likely because only young wind nebulae have a detectable extended X-ray counterpart, which allows for a firm identification. Most of the candidates have previously been hypothesised to be a PWN or to have a PWN component. The only substantially older pulsar is PSR B1742−30, which is selected thanks to its very low distance despite its low Ė. We cannot display this pulsar in all plots of this paper, but we discuss it as a special case in Sect. 6.

Fig. 2 Left: spin-down power Ė and characteristic age τc of pulsars with a firmly identified PWN, candidate PWN, and without TeV counterpart (grey dots). The black line and shaded band show the injection evolution of the modelling used in this paper. The dashed lines indicate lines of constant total remaining energy Ė τ; see Appendix B. Hence a model curve that starts at Ė0τ0 = 1049 erg represents a pulsar with total initial rotational energy of 1049 erg. Since both Ė and τc depend on P and Ṗ , the axes in this plot do not represent independent quantities. Right: same data, shown in the commonly used view, using the independently measured P and Ṗ .

3.3 Location in the Galaxy

In order to assess the reach of the population study presented in this work it is instructive to display the positions of Galactic PWNe together with the sensitivity (or depth) of the H.E.S.S. Galactic Plane Survey. The map in Fig. 3 visualises the 2D projection of the Galactic distribution of very energetic pulsars (Ė > 10³⁵ erg s-1). The symbols distinguish between pulsars with firmly identified wind nebulae, candidate PWNe, and pulsars at >10³⁵ erg s-1 for which no TeV wind nebula has been detected so far. For reference, the map comprises a schematic representation of the spiral arms of the Milky Way according to the parametrisation of Vallée (2008). The overlaid blue and yellow curves define the accessible range of the HGPS for point-like sources at an integrated luminosity (1–10 TeV) of 1% and 10% of the Crab luminosity, respectively (for details see H.E.S.S. Collaboration 2018).

For sources of 10 % Crab luminosity, the HGPS covers approximately one quarter of our Galaxy, and generally does not reach much farther from Earth than the distance to the Galactic centre. For extended objects, the horizon can be expected to be closer, and for close-by extended sources, the H.E.S.S. FOV can limit the capability of isolating them from the background.

Most of the detected PWNe are located close to one of the nearby dense spiral arm structures, where pulsars are expected to be born. In particular, the Crux Scutum arm hosts half of all HGPS pulsar wind nebulae. Several high-Ė pulsars are on closer spiral arms but are not detected.

A way to look at the sensitivity to extended PWNe is shown in the upper part of Fig. 4, where the extension is plotted against distance from Earth. To guide the eye, two lines indicate the range of detected extensions between the systematic minimum of about 0.03° and the maximum extension in HGPS of ~0.6° (Vela X, see Sect. 5.2.1). As can be seen in the lower panel of Fig. 4, most PWNe are detected around 5.1 kpc, which is theaverage distance of PWNe in Table 1. This allows for the determination of radii between 3 and at least 60 pc.

We conclude that both the H.E.S.S. FOV (5°) and angular resolution (0.03°) are adequate to study the wind nebulae of most of the high-Ė pulsars known today.

Fig. 3 Schematic of the Milky Way and its spiral arms, along with firmly identified PWNe, candidates, and energetic pulsars (Ė > 1035 erg s-1) without detected TeV wind nebula. The yellow and blue curves outline the sensitivity horizon of the HGPS for point-like sources with an integrated gamma-ray luminosity (1–10 TeV) of 1% and 10% of the Crab luminosity, respectively (see H.E.S.S. Collaboration 2018, for details).

Fig. 4 Top: PWN extension occurrences over distance from Earth, in comparison to the band of extensions that can be expectedto be identified in the HGPS analysis chain. Bottom: distribution of known distances of energetic pulsars (Ė > 1035 erg s-1).

4 Theoretical notion of pulsar wind nebulae

Before discussing the properties of the PWNe and PWN candidates we found, this section recapitulates some concepts of the theoretical understanding of pulsar wind nebulae.

A PWN is usually considered to be a calorimetrical, dynamical object around a pulsar. It stores and displays the radiative output of the pulsar during tens of kiloyears while at the same time undergoing a substantial dynamical evolution inside the host SNR. Expressed in terms of a diffusion equation, this means that it is energised by the magnetic and particle flux from the pulsar, and cooled by radiative (synchrotron emission and IC scattering), adiabatic, and escape losses (e.g. Martín et al. 2012; Zhang et al. 2008, and references therein). In the context of this work, acceleration and injection mechanisms are not considered in detail. Pulsars are regarded as particle-dominated, diffuse injectors of electrons. Here and in the following, the term “electrons” always refers to the full electron and positron outflow.

4.1 Injection evolution

The energy outflow of the pulsar Ė determines the energy injection history of a PWN. It is decaying continually at a rate determined by the so-called spin-down timescale τ, following anevolution similar to that expected from a dipole (see also Appendix B) $$\text{edot}(t)={\text{ edot}}_0{\left(1+\frac{t}{{\tau }_0}\right)}^{-\frac{n+1}{n-1}},$$(3)

where τ₀ is the initial spin-down timescale, Ė₀ is the initial spin-down luminosity, n is the so-called “braking index” (e.g. Pacini & Salvati 1973), and t is the time since the birth of the pulsar. Values typically considered are τ₀ ~ 10^2.5−3.5 yr, Ė₀ ~ 10^37.5−40 erg s-1, and n ~3 (Martín et al. 2012; Zhang et al. 2008; Vorster et al. 2013; Gelfand et al. 2009). This indicates that most of the pulsar rotational energy budget (${\varphi }_0\times {(E/E_0)}^{-\Gamma }$, typically ≲10⁵⁰ erg; see Appendix B) is spent in the first few thousand years.

The present spin-down luminosity can be calculated from the period P and its time derivative Ṗ (Gaensler & Slane 2006, Eq. (1)). Another parameter that can be derived from the pulsar ephemeris is the so-called characteristic age, which is defined as ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$(4)

If t ≫ τ₀ and n = 3, then τc is an estimator for the true age t of a pulsar. Independent of this condition, though, Eqs. (3) and (4) imply a straight power-law correlation between Ė and τc , i.e. $\text{edot}={\text{ edot}}_0{\left[\frac{2}{n-1}\cdot \frac{\text{age}}{{\tau }_0}\right]}^{-\frac{n+1}{n-1}},$(5)

or, equivalently, between Ṗ and P (see Eq. (B.12) in Appendix B), i.e. $$\text{pdot}(P)=\frac{P_0}{{\tau }_0}\frac{1}{n-1}{\left(\frac{P}{P_0}\right)}^{2-n}\cdot$$(6)

Consequently, the power indices of the above relations are only determined by the braking index n. Figures 2 show howreal pulsars populate these diagrams. They are born on the upper left of the plots and move towards the lower right as their spin-down decays. Pulsar population synthesis studies have shown that this distribution can be reproduced assuming magnetic dipole spindown (n = 3; e.g. Faucher-Giguère & Kaspi 2006, and references therein). Some such studies found evidence for pulsar magnetic field decay, but on timescales of several Myr (e.g. Gonthier et al. 2004). As this is much longer than the PWN evolution timescales we consider, in the baseline model of this paper we assume that the injection evolution is dictated by an average braking index n = 3, which is a compromise between theoretical expectation, observed pulsar Ė and τc, and the measured braking indices (see Appendix A for more details).

4.2 Dynamical evolution

The dynamical evolution of PWNe can generally be divided into three distinct stages (Gaensler & Slane 2006; Gelfand et al. 2009; van der Swaluw et al. 2001, 2004, and others): the free expansion (<2–6 kyr), reverse shock interaction (until some tens of kyr), and relic stage. In the free expansion phase, the plasma bubble grows inside the unshocked ejecta of the SNR, whose forward and reverse shocks do not interact with the PWN. This phase is comparably well understood because of numerous analytical (Rees & Gunn 1974; Kennel & Coroniti 1984a,b) and numerical (Martín et al. 2012; Bucciantini 2011, and references therein) works on the subject mostly focussed on the Crab nebula case, but applicable to other young PWNe. The PWN is growing fast (Chevalier 1977, R ~ t^1.2), attenuating the magnetic field strength and synchrotron radiation, while IC emission from the accumulating electrons quickly increases in the beginning and then decreases very slowly (Torres et al. 2014). This early stage is the only phase where the IC scattering on synchrotron photons (synchrotron self-Compton emission) can also play a dominant role.

The second phase begins after a few thousand years, when the PWN has grown to a size of the order of ~ 10 pc and encounters the reverse shock of the SNR, which may be moving spatially inwards (Blondin et al. 2001). Since the total dynamic energy in the SNR exceeds that of the PWN by one or two orders of magnitude, the PWN may be compressed again by up to a factor of 10 (Gelfand et al. 2009) and experiences a series of contractions and expansions until a steady balance is reached. After that, the wind nebula continues to grow at a much slower pace, like R ~ t^0.73 for t < τ₀ in van der Swaluw et al. (2001) and R ~ t^0.3 for t > τ₀ in Reynolds & Chevalier (1984). In the work of Gelfand et al. (2009), where a spherically symmetric case was simulated, the oscillations were found to lead to dramatic changes in the synchrotron and IC luminosities, making the TeV emission disappear completely for several thousand years. In reality, where the SNR develops asymmetrically and the pulsar has a proper motion, these drastic changes are presumably washed out to some degree, leading to a more continuous behaviour. Still, the collisionof PWN bubble and reverse shock heavily depends on the evolution of the whole system and its interaction with the surroundings, making such evolved PWNe very diverse, non-uniform objects (see also de Jager & Djannati-Ataï 2009).

This non-uniformity becomes even more pronounced if the pulsar, owing to its proper motion or a tilted crushing of the nebula, spatially leaves the main PWN bubble or even the SNR. In that case, which is called the relic stage, the pulsar can form a local plasmabubble while the old nebula from its younger period still remains, typically as an IC-dominated PWN due to its much lower magnetisation.

4.3 Modelling

The interpretation of the data we present and of the log-linear trends we fit to the evolution plots require a comparison to what can be expected in theory with the basic concepts outlined above. To do so, we built a simpified, time-dependent model for the evolution of the VHE electron population and TeV emission of PWNe. We deliberately opted for a simple model because we do not need it to contain detailed parameters that our TeV data does not allow us to investigate.

The model we describe in Appendix A assumes a time-dependent injection of electrons with a fixed power-law spectrum⁷ , but decreasing total power according to Eq. (3). Following analytical formulae for the expansion, the cooling from synchrotron, adiabatic, inverse Compton, and escape losses is applied to the electron population as a function of time. The respective characteristic age τc is always tracked as well to compare the model correctly to data. The photon emission is calculated for each time step from the electron population, including the full Klein-Nishina formula.

The strategy for the comparison of PWN data and theory is to define the parameters of the model such that it reflects both the average trend of PWN evolution (baseline model) and the scatter of individual wind nebulae around that average expectation (varied model). This means that, unlike other works, we do not model individual objects in their particular multi-wavelength context. Instead, we try to find out what the typical evolution is and what the typical variations need to be in order to produce the picture we obtain for the whole population. The band of the varied model can therefore be interpreted as the area where a synthesised population would be found (in the absence of detection selection effects).

As it turns out in the following, we succeeded in finding such a model describing the evolution that a typical PWN in a typical, dense spiral-arm surrounding undergoes. Since this one model implies an evolution curve for every observable we consider, both along τc and Ė₀ , a good leverage on its absolute parameters is given. Starting from the baseline model, the parameters are varied with the aim to realistically reproduce the scatter of measured PWN observables. This way, the scatter of observables itself is exploited as another observable, with the large number of curves leading again to a good handle on the scatter.

It should be noted though that intrinsic (physical) and analytical (mathematical) correlations between parameters are neglected in the varied model. For instance, the scatter ranges of Ė₀ and τ₀ , strongly restricted by Fig. 2 (left), may be larger if the two quantities were anti-correlated such that high-Ė₀ pulsars always tend to havea lower τ₀; this is physically plausible because the two quantities are related through the pulsar birth period and magnetic field. On the mathematical side, Ė₀, η, the energy injection range and the background photon density are all parameters with which the TeV luminosity scales in an almost linear way. In our varied model, we deal with this redundancy by only varying Ė₀ , but similar results can be achieved if one of the other factors is varied instead. See also Appendix A.7 for this and other caveats of the model.

5 Properties of TeV pulsar wind nebulae

In this section we present and discuss the distributions and correlations of TeV wind nebulae and their respective pulsars. For each topic we describe what we present, discuss potential biases, and then interpret what we find, using the modelling described in Appendix A where needed and appropriate. The presented plots serve to evaluate the plausibility of our current candidate sample (Sect. 6) and may prove useful in investigating future PWN candidates.

5.1 Fitting and statistical treatment of uncertainties

The properties of PWN are intrinsically scattered (see Sect. 4.2) and all observables are calculated using a distance estimation based on the dispersion measure of the pulsar and a model of the Galactic free electron distribution, whose uncertainty is not statistically well described. Consequently, the probability density functions (p.d.f.) of our observables (size, luminosity, and offset) for a given τc or Ė are dominated by the scatter of intrinsic properties and errors in the distance estimation and not by our statistical uncertainty.

As a consequence, in the cases where we pursue a fit of observables with the aim of testing the significance of a correlation or extracting an estimator function, we follow the approach put forward by Vink et al. (2011) and Possenti et al. (2002). They performed a least-squares fit of the respective observable with residuals calculated in common logarithmic space. The fit function is a (log-)linear function, expressed generally as $$\text{lgt}{Y}_{\text{est}}=p_0+p_1\text{ lgt}X.$$(7)

In order not to be restricted to detected objects but also to include the valuable limits from pulsars without VHE emission, we use the asurv code (Lavalley et al. 1992) for the minimisation. It allows us to apply statistical methods to test for the existence of a correlation, such as the Cox proportional hazards model, or to perform a multivariate regression including limits (see Isobe et al. 1986, for an overview on the statistics inside asurv). Besides the parameters p_i of our function, asurv also determines the variation σ_{lg Y} that the data are scattered with.

Owing to the existing selection biases and the uncertain p.d.f. shapes involved, the derived estimator function might not always approximate a virtual true evolution function, but rather evaluate the unweighted average trend of the examined data points. Table 6 summarises the fit results that are referred to in the following paragraphs. The p-values are taken from the Cox proportional hazards model, which is a regression method for data with upper limits. This model was originally developed for biostatistical applications, where it is extensively used. As described in Isobe et al. (1986), Section III, the model provides an equivalent χ² for the null hypothesis (no correlation), which can be transformed to a p-value. For the linear regressions and parameter determinations, the expectation maximisation (EM) algorithm is used, which is an iterative least-squares method that allows for the inclusion of limits (Isobe et al. 1986, Sect. IV).

5.2 Morphological properties

The morphological parameters provided by the HGPS catalogue are source position and extension. As a pulsar and its PWN evolve, the PWN is thought to become increasingly extended and offset from the pulsar position (see Sect. 4). This basic evolutionary behaviour can be found unmistakably in Figs. 5 and 6 (left).

Fig. 5 Left: PWN extension evolution with time, in comparison to the modelling considered in this work. Right: PWN extension evolution with Ė, as fitted in the RPWN(Ė) column of Table 6 for pulsar wind nebulae with Ė > 1036 erg s-1 (see Sect. 5.2.1). The shaded range shows the fit range and standard deviation σlg R . 1 dex refers to an order of magnitude and is the unit of the logspace defined σlg Y . For clarity, this plot excludes PWN candidates and divides the sample into nearby and far pulsar wind nebulae to illustrate the potential selection or reconstruction bias (see text). The dot-dashed and dotted lines indicate the systematic minimum of 0.03° and the maximum measured extension in the HGPS of 0.6°, respectively, which are both projected to the average PWN distance of 5.1 kpc.

Fig. 6 Left: spatial offset of pulsar and TeV wind nebula as a function of pulsar characteristic age. The dashed and dotted lines indicate the systematic minimum of about 0.015° and the association criterion of 0.5°, which are both projected to the average PWN distance of 5.1 kpc. The solid line shows the offset one would expect from pulsar motion only (assuming a large pulsar velocity of 500 km s−1 ). Right: time evolution of the containment ratio. Since the pulsar motion can be assumed to be constant and the expansion decelerates, one expects the containment fraction to increase and eventually pass unity after some tens of kiloyears. The dotted horizontal line shows the rating criterion (offset/extension <1.5) applied in the post-selection of candidates (Sect. 6).

5.2.1 Extension

Figure 5 (left) shows the evolution of PWN extension as a function of characteristic age τc . We can determine extensions beyond a systematic minimum of around 0.03° and at least up to the observed extension of Vela X, at around 0.6° (see Sect. 3.3). As shown in Fig. 4, most known pulsars lie at distances that therefore allow for the measurement of PWN extensions between 3 to 60 pc. In Fig. 5 (right), where the extensions are plotted against pulsar spin-down, far and close-by systems are distinguished. This elucidates our ability to resolve far and near systems and shows the plain correlation of size with Ė .

A caveat is that there is a selection bias from the fact that extension estimates or limits are only available for sources that are detected. Systems that are too faint or too large to be detected with our sensitivity and FOV are missing in the PWN sample. Since we cover a wide range of different distances, sources that are large and bright, or faint and small, can still be represented to some level in the sample. However, if there is a state in which PWNe are faint and large at the same time, it might be that they cannot be detected at any distance. From the current understanding of PWN theory, this can be the case for PWNe of ages beyond few tens or hundreds of kiloyears, so the study presented here has to be taken with some caution in that regime. To unbias the sample in the fitting procedure below, we apply a cut of Ė > 10³⁶ erg s-1, beyond which the likeliness of detection is reasonably high and the detected objects can be considered representative for their stage of evolution.

A measurement bias we may have is that the limited FOV might truncate the tails of the source for very extended sources. This effect was suggested by Vernetto et al. (2013) as an explanation for the differences between some IACT spectra and the results of the air shower detector ARGO. We cannot entirely verify or falsify this claim here, but since only few sources approach the critical regime beyond 1°, it is presumably a minor effect in this study.

A possible physics bias that might enhance the effect seen in Fig. 5 (right) is that close-by objects are on average located farther away from the Galactic centre and therefore in less dense surroundings than far objects. This might influence the average dynamical evolution they experience.

Fitting the data to check for correlations with τc or Ė yields the results shown in Table 6. The low p-values and non-zero p_1,2 confirm, on 2–3 standard deviation confidence levels, that the extension increases along the evolution of a PWN, i.e. with falling Ė and increasing τc . A more general 2D fit of R_PWN(P, Ṗ) does not lead to a significant improvement of the fit, nor a lower p-value. The parametrisation of R_PWN(Ė) is shown in in Fig. 5 (right) to show that it is indeed suitable for predicting the extensions of the detected young PWNe (Ė > 10³⁶ erg s-1) reasonably well. The only PWN below 10³⁶ erg s-1, CTA 1, does not follow the extrapolation of that trend and appears to be dynamically different from the rest of the population.

The relation R ~ τc^0.55±0.23 can be compared to the baseline model in Fig. 5 (left), which assumes the canonical R ~ t^1.2 and t^0.3 , at early and late times, respectively, and thus encloses the measured value well. The conversion between true age and τc according to Eq. (4) is taken into account in the displayed model curves.

Comparing the data with the model, the initial and fast free expansion can accommodate the non-detections of extensions of very young pulsar wind nebulae, while the slope of t^0.3 for evolved PWNe (Reynolds & Chevalier 1984) is roughly consistent with the comparably small extensions of the few older PWNe in the sample. One has to keep in mind that the curve at high ages is more an upper limit than a prediction because a potential crushing (as a sudden decrease in size after the free expansion) is not included in the model. The absolute scale of the curves is a free parameter of the model, but later turns out to be constrained by the surface brightness values we measure (Sect. 5.3.3).

In conclusion, the fact that TeV pulsar wind nebulae generally grow with time until an age of few tens of kiloyears is clear and supported by the asurv fits. There are, however, few pulsar systems older than that to place stringent constraints on the model at later evolution stages.

5.2.2 PSR-TeV offset

An offset between a pulsar and its TeV wind nebula can be caused by a combination of pulsar proper motion, asymmetric crushing of the PWN by the surrounding SNR and, possibly, by asymmetric pulsar outflow. Hobbs et al. (2005) determined the mean 2D speed for non-millisecond pulsars to be 307 ± 47 km s-1, and the velocity distributions were found to be compatible with a Maxwellian distribution. Other works, such as Arzoumanian et al. (2002), suggest a more complex distribution and high-velocity outliers, but it is clear that the bulk of the pulsars have 2D velocities of less than 500 km s-1. In Fig. 6 (left), the offset against characteristic age is compared to a 2D velocity of 500 km s⁻¹. The true age of young pulsars can be less than τc, in which case those points, shown in true age, may even move to the left, and thus enlarge the distance to the 500 km s⁻¹ line. To give an idea of which offsets can be detected, lines for the maximum offset implied by our angular selection criterion and systematic minimum resolution are also shown for the mean PWN distance of 5.1 kpc.

The asurv fits (Table 6) suggest that the trend of increasing offset with rising τc and falling Ė is statistically manifest in the data. What is interesting beyond this general increase is that 5 of 9 pulsars with ages beyond 7 kyr in Fig. 6 (left) are more offset from their PWN than expected from mere pulsar motion. At these ages, PWNe presumably are beyond their free expansion phase and have started interaction with the SNR reverse shock or surrounding medium. While the velocity distribution of pulsars can have outliers that are significantly faster than average, it is unlikely that such a high fraction of the high-Ė pulsars with TeV-detected nebulae are so fast (1000 km s-1 or more would be required). This suggests that the asymmetric evolution of the PWN, caused by interaction with the reverse shock and/or asymmetric surrounding medium (Blondin et al. 2001), is in fact the dominant offset mechanism for middle-agedwind nebulae. Further support for this conclusion comes from the very few measured pulsar transverse velocity vectors that are currently available in the ATNF catalogue for our PWN sample (e.g. for Vela X and HESS J1825 − 137). These vectors do not consistently point away from the PWN, as one would expect from a pulsar motion dominated offset.

5.2.3 Containment

Containment of a pulsar in its TeV wind nebula, although not strictly binding in the relic stage, is often taken as an argument to claim the PWN nature of an object. We define the containment ratio as the PWN offset divided by the PWN extension radius. Given the offset and extension evolution discussed above, a pulsar is not expected to leave its (then relic) wind nebula before some tens of kiloyears; yet the ratio should increase and approach unity at some point, unless the relative movement is in the direction of the line of sight. Figure 6 (right) shows the evolution of the containment ratio with characteristic age.

An additional caveat to mention for this quantity is that no upper or lower limit can be calculated if both offset and extension are already limits, which is the case for 7 of the 19 firmly identified objects in our sample. Another selection bias concerns the identification itself. Good reasons, such as observations at other wavelengths, are required to argue for a non-contained association of a pulsar with a TeV object; for old systems, however, these MWL data are very difficult to acquire since the synchrotron component has become very faint. This bias can be regarded as intrinsic to the decomposition of old, “dissolving” PWNe, whose remains become inevitably difficult to associate with the pulsar as time passes.

In Figure 6 (right), most young pulsars are well contained in their nebulae, but there are a few older pulsar wind nebulae that were firmly associated to a pulsar close to or slightly beyond their (1σ Gaussian) extension radius.

5.3 Luminosity, limits, and derived parameters

5.3.1 Luminosity

From Sect. 4 and in our model, the TeV luminosity of pulsar wind nebulae is expected to rise quickly within the first few hundred years and decay slowly over many thousands of years. Figure 7 (left) shows the evolution of luminosity with pulsar spin-down power and Fig. 8 (left) the evolution with characteristic age.

Figure 8 (right) indicates the distribution of luminosities. The average detection threshold of energy flux between 1 and 10 TeV is at around 10⁻¹² erg s-1 cm^-2 (H.E.S.S. Collaboration 2018), which is equivalent to a luminosity threshold of 3 × 1033 erg s-1 at the mean PWN distance of 5.1 kpc. In this work, we reduce the selection bias present in previous studies by involving flux upper limits for all eligible pulsars with Ė > 10³⁵ erg s-1. As explained in Sects. 2.3 and 3.2, about one-third of these high-Ė pulsars in the ATNF catalogue can be expected to have an extension small enough from which it is possible to extract a meaningful limit. So again, PWNe that are very large, presumably with ages beyond a few tens of kiloyears (below ~ 10³⁶ erg s-1), might be truncated from our data set. Figure 7 (right) is the equivalent of Fig. 5 (right), showing the luminosities and limits in two bands of distance. The expected extensions and derived limits are listed in Table 5. In the fit below we add further flux limits calculated for the pulsars associated with the candidate PWNe, applying the same calculation method as for the limits in Table 5. This adds 11 further valid limits, which are also included in Fig. 7 (right). This flux limit can actually be below the flux of the candidate, for instance if the candidate is more extended than predicted by the model (e.g. in the case of HESS J1023−575).

The primary feature of the data is a mild but stable correlation of luminosity with pulsar spin-down⁸ . The asurv fit suggests a relation of L_{1−10 TeV} ~Ė^0.59±0.21 (see Table 6). The model supports this, indicating a power index of around 0.5. The slow but steady decay, combined with the growing extension, is what hampers a TeV detection once the pulsar spin-down power falls below ~10³⁶ erg s-1. This decay could not be observed in other works before (e.g. Mattana et al. 2009; Kargaltsev & Pavlov 2010) owing to the missing upper limits⁹. Figure 7 (right) shows the result of the fitted parametrisation L_{1−10 TeV}(Ė) derived from our data.

In contrast, the L_{1−10 TeV} over τc (Fig. 8, left) is scattered widely and a correlation is statistically not clear (see Table 6). This, however, matches the broad scatter suggested by the varied model (shaded area). Apparently, Ė is the better variable to characterise the evolutionary state of the PWN luminosity.

Fig. 7 Left: relation of TeV luminosity and pulsar Ė. Right: same as Fig. 5, but for luminosity and including the upper limits from the left panel and 11 upper limits calculated for pulsars of candidate PWNe (see text).

Fig. 8 Left: evolution of TeV luminosity with characteristic age. Right: distribution of 1–10 TeV luminosity for PWNe and PWN candidates treated in this work.

5.3.2 Apparent TeV efficiency

The TeV efficiency, conventionally defined as εTeV = L_{1−10 TeV}∕Ė, is not the real present efficiency of a PWN because L_{1−10 TeV} is a result of the whole injection history, whereas Ė characterises the present outflow of the pulsar. Therefore, TeV pulsar wind nebulae can in principle have TeV efficiencies exceeding unity.

Figure 9 (left) shows the evolution of the efficiency with the pulsar characteristic age. Interestingly, the efficiency seems to be scattered more than suggested by the varied model, unlike in Fig. 8 (left). To shed light on the cause of this it is illustrative to plot TeV efficiency versus the PSR-PWN offset for different groups of characteristic age (Fig. 9, right). With the sample of detected PWNe, a relatively clear correlation can be confirmed also in the asurv fit (Table 6). Apparently, all low-efficiency PWNe are found at low offsets from their pulsar and all high-efficiency wind nebulae have larger offsets. To some level, this correlation is trivial because both efficiency and offset increase with time. After subdividing the sample into different age groups, however, it becomes clear that the plot does not only sort by age; instead, even for PWNe with similar ages, efficiency and offset are correlated and the age groups overlap each other. In the plot, a bias might occur because low-efficiency, high-offset systems may be difficult to identify, but the absence of high-efficiency, low-offset systems must be genuine. A second systematic effect may be that both efficiency and offset depend on the PSR distance estimation d (or its square), so if there were a strong bias in d, it would also appear as a trend in the plot (though hardly at scales of more than a factor of 10).

This is a mild indication that a pronounced offset, as may be induced by SNR reverse shock crushing, comes with a high TeV efficiency, while systems that interact less with their surroundings remain fainter overall. One cannot disentangle at this point whether the crushing itself heats up the plasma bubble or whether the correlation is indirect because denser environments might provide both crushing material and a higher IC target photon density. More case studies in the future might clarify the situation here.

Fig. 9 Left: evolution of TeV efficiency (L1−10 TeV∕Ė) with pulsar characteristic age. Right: TeV efficiency as a function of pulsar offset, plotted for pulsars of different age groups. High-offset systems tend to be more TeV-efficient than low-offset systems.

5.3.3 Surface brightness

All of the TeV quantities discussed so far rely on the knowledge of the distance to a given pulsar system, which in many cases, however, is not very well constrained observationally. A quantity that is independent of the distance is the TeV surface brightness, defined as $$S=\frac{\text{lumi}}{4\pi R\text{ind}PW{N}^2}\approx \frac{F_{1-10\text{TeV}}}{{\sigma }^2},$$(8)

where R_PWN is the physical PWN radius (in pc), σ is its angular extent as seen from Earth, and F_1−10TeV is the integral energy flux between 1 TeV and 10 TeV measured atEarth.

Figure 10 (left) shows the dependence of surface brightness on the pulsar’s Ė . Like the extension, S can only be calculated for detected systems and therefore suffers a selection bias expected to become more important with decreasing Ė . Below a spin-down power of 10³⁶ erg s-1, the data sample is truncated at low surface brightness values.

As seen in the asurv fit values in Table 6, a comparably strong correlation is found, confirming the above findings of a decreasing luminosity and increasing extension of ageing pulsar wind nebulae. The measured power-law relation of S ~Ė^0.81±0.14 matches what the model suggests (~0.9 for the part where Ė < 2 × 1037 erg s-1). We find that the surface brightness gives a strong handle on the self-consistency of the model because it links the dynamical evolution (i.e. the extension) to the spectral evolution (i.e. the flux). That is, the scales of the extension and luminosity evolutions cannot be adjusted independently; they must lead to a consistent surface brightness scale.

An interesting feature to note is that the scatter suggested by the varied model seems to be much larger than what is found in the data (σlg S ~ 0.3). This might indicate that flux and extension are not as independent as implied by a free variation of the respective model parameters. Another effect might be the missing systematic scatter of S from the distance measurements. If the scatter of luminosity and extension measurements were dominated by the errors on the distance, the varied model shown here would implicitly include that scatter, and therefore overestimate the actual source-intrinsic scatter. This in turn would lead to a spread of the predicted surface brightness evolution that is too large.

Fig. 10 Left: dependence of surface brightness S (see text) on pulsar spin-down Ė. Right: TeV photon index over pulsar spin-down Ė for all detected PWNe and candidates.

5.3.4 Photon index

The average photon index in our sample (firm identifications) is ~ 2.3, and about half of the PWN indices deviate significantly from that. Figure 10 (right) shows the relation of photon index and pulsar Ė.

A selection bias can be expected because non-detections do not appear in the plot and very soft spectrum sources are more difficult to detect than hard spectrum sources.

The general range of measured indices (1.9…2.8) is in accordance with the model; most of the firm identifications lie in the predicted range of the varied model or have error bars that are compatible with this model. The precise index is a product of the lepton spectral energy distribution, in particular of elderly cooled electrons (see Fig. A.1, right) and the IC target photon fields, the combination of which on average seems to be appropriate in our model. The two exceptions are the Crab nebula and N157B, for which TeV emission is likely dominated either by IC scattering off their own synchrotron radiation (SSC, for Crab), or dominated by a very high surrounding photon field (N157B). These special features are not incorporated in our generic model.

The peak of the IC emission does not have a clear tendency in our model, although a mild trend for an increasing peak position seems to be manifest in Fig. A.1 (left) beyond ages of few kiloyears. Also, such a trend is not generally agreed on between different modelling codes. The MILAGRO and HAWC observations of the ancient Geminga PWN indicate a multi-TeV nebula (Abdo et al. 2009; Baughman et al. 2015), presumably with a high-peaking spectrum, despite its age of ~ 300 kyr. The data discussed in this paper do not allow for a clear statement here, but show that the trend, if present, is weak.

6 Review of pulsar wind nebula candidates

In the previous section, the preselection candidate PWNe defined in Sect. 3.2 were shown along with the firmly associated PWNe and with average model expectations. Some of the candidates very consistently lie among other PWNe and close to the model prediction, while others do not. In order to compare the candidates among each other, in this section we apply uniform post-selection (“rating”) criteria to all of them.

It is important to note that such a rating only evaluates the plausibility of a given candidate in the context of firmly identified PWNe or of our model (which is adjusted to the PWNe). Therefore, a badly rated candidate may either be an atypical PWN, or an object that contains a PWN alongside a second source (such as a stellar cluster or SNR), or no PWN at all. Arguments from observationsat other wavelengths are ignored in this uniform approach here since they are not available for all candidates. Consequently, our rating evaluates the plausibility of a PWN candidate by how normal the TeV properties of the PWN candidates are.

We evaluate four criteria: three are comparisons to the model evolution and one concerns the containment of the pulsar inside the PWN. Specifically, we apply the following criteria:

• 1.

  Containment ratio (Fig. 6, right): the pulsar offset should be <1.5 extension radii.

• 2.

  TeV extension versus age (Fig. 11, left): Log-residual from model (Fig. 5, left) should be within 2 standarddeviations, using the measured σlg R = 0.39 (Table 6).

• 3.

  TeV luminosity versus pulsar spin-down (Fig. 11, middle): Log-residual from model (Fig. 7, left) should be within 2 standarddeviations, using the measured σlg L = 0.83 (Table 6).

• 4.

  Surface brightness versus pulsar spin-down (Fig. 11, right): Log-residual from model (Fig. 10, left) should be within 2 standard deviations, using the measured σlg S = 0.30 (Table 6).

Table 4 shows the ratings of the considered candidates. There are 20 PSR-TeV pairs, in which there are two TeV double associations (one HGPS source qualifying for two pulsars) and four PSR double associations (one PSR qualifying for two HGPS sources).

Ten of the candidate PSR-TeV pairs fulfill all criteria and seem to be plausible TeV pulsar wind nebula associations. All of these candidate pairs have already been discussed as possible TeV PWNe, namely HESS J1616−508 and HESS J1804−216 (both in Aharonian et al. 2006c), HESS J1809−193 and HESS J1718−385 (both in Aharonian et al. 2007), HESS J1857+026 (Hessels et al. 2008), HESS J1908+063 (aka MGRO J1908+06; e.g. Aharonian et al. 2009; Aliu et al. 2014), HESS J1640−465 (Abramowski et al. 2014, PWN hypothesis disfavoured, though), HESS J1708−443 (H.E.S.S. Collaboration 2011a), HESS J1023−575 (coinciding with massive stellar cluster Westerlund 2, H.E.S.S. Collaboration 2011b), and HESS J1018−589B (the extended additional component close to the binary HESS J1018−589A, H.E.S.S. Collaboration 2012b).

Of the ten disfavoured candidates, one is an alternative association for the above strong candidate HESS J1616−508, disfavoured due to its offset. Similarly, PSR J1811−1925 is a second pulsar in the area of HESS J1809−193, already argued in Aharonian et al. (2007) to be the less likely counterpart of the two pulsars that can be considered. HESS J1026−582 was previously hypothesised to be a PWN, but receives an unfavourable rating due to its pulsar offset, although the HGPS analysis may not be optimal to reveal the morphology of this hard-spectrum source (H.E.S.S. Collaboration 2011b).

The two sources HESS J1745−303 and J1746−308, both associated with the very nearby old PSR B1742−30, are a special case. The pulsar is a factor of 10–100 older than most other PWNe discussed here, so the extrapolation performed for the rating cannot be considered to be very robust. In fact, these two objects could not be represented in most of the figures because they are too far off the axis ranges. They obtain a bad rating mostly because they are both too underluminous and too small for their age. It could well be, though, that HESS J1746−308 is a late-phase PWN, created locally near the pulsar after the main relic PWN bubble has become very faint and/or has dissolved. The predicted size of the PWN according to our model would be 32 pc or 9° in the sky, which is impossible to detect with state-of-the-art IACT analysis methods.

In conclusion, about half of the PWN candidates evaluated in this work are viable PWNe, judging by their TeV and pulsar properties in relation to the population as such. The number of disfavoured candidates (10) matches well with the expectation of ~ 10 chance coincidences evaluated in Sect. 3.1. Hence, it seems plausible that most of the ten high-rated candidates are indeed genuine pulsar wind nebulae. If this were the case, a total of 25 in 78 HGPS sources would be pulsar wind nebulae (including G0.9+0.1 here).

Fig. 11 Common logarithmic residuals of rating criteria 2–4, using the standard deviations σlg Y explained in Sect. 6. Left: extension with respect to the model shown in Fig. 5 (left). Middle: same for luminosity (Fig. 7, left). Right: same for surface brightness (Fig. 10, left). In all cases, limits are shown separately as outlined histograms.

7 Conclusions

In this workwe subsume and examine the population of TeV pulsar wind nebulae found to date. The census presents 14 objects reanalysed in the HGPS catalogue pipeline, which we consider to be firmly identified PWNe, and five more objects found outside that catalogue range or pipeline. In addition to those, we conclude that there are ten strong further candidates in the HGPS data. Most of the PWNe are located in the bright and dense Crux Scutum arm of the inner Milky Way. A spatial correlation study confirmed the picture drawn in earlier studies, namely that only young, energetic pulsars grow TeV pulsar wind nebulae that are bright enough for detection with presently available Cherenkov telescopes. For the first time, flux upper limits for undetected PWNe were given around 22 pulsars with a spin-down power beyond 10³⁵ erg s-1 and with expected apparent extensions (plus offsets) below 0.6° in the sky.

Of the 17 most energetic ATNF pulsars, with a spin-down power Ė ≥ 10³⁷ erg s-1, 11 have either an identified TeV wind nebula (9) or candidate (2) featured in the present study. Of the remaining 6,

• 3 are included in the flux upper limits in Table 5;

• 3 are out of the range of the HGPS:

  • PSR J2022+3842: SNR G076.9+01.0, contains an X-ray PWN; not reported in TeV.

  • PSR J2229+6114: boomerang, contains an X-ray PWN; detected by MILAGRO and VERITAS, but of unclear nature in TeV.

  • J0540−6919: in the Large Magellanic Cloud; a limit is given in Abramowski et al. (2015). Converting the limit to luminosity yields L_{1−10 TeV} < 5.7 × 1034 erg s-1, which is compatible with the predicted 3.3 × 1034 erg s-1 that can be taken from L_{1−10 TeV}(Ė) in Table 6.

In summary, only 5 of the 17 highest-Ė pulsars remain without a detected potential counterpart in the TeV band.

Figures 5 to 10 showed a variety of trends between pulsar and TeV wind nebula parameters, and consistently compared them to a simple one-zone time-dependent emission model of the TeV emission with a varied range of model input parameters. The main conclusion of this work is that for several observables, a trend was found in the data and the trends suggested by our model are consistent with these findings. With only a moderate variation of the model input parameters, we can mimic the spreads of the observables, although the precise value of the parameter ranges is subject to the model caveats discussed in Sect. 4.3. Our first-order understanding of the evolution of TeV pulsar wind nebulae with ages up to several tens of kiloyears therefore seems to be compatible with what the whole population of detected and undetected PWNe suggests.

Using the flux limits for undetected PWNe, we find evidence that the TeV luminosity of PWNe decays with time while they expand in (angular) size, preventing the detection of those whose pulsar has dropped below ~ 10³⁶ erg s-1 (roughly corresponding to several tens of kiloyears). This was implicitly known before from the mere non-detection of old pulsar wind nebulae, butfor the first time could be put into a quantitative perspective here, both by fitting data and limits, and by comparing the data to model predictions. The power-law relation between TeV luminosity and pulsar spin-down power could be estimated as L_{1−10 TeV} ~Ė^0.58±0.21, in consistency with the model that suggests a power index of around 0.5.

Anotherfeature that was discussed on some individual objects before (e.g. Aharonian et al. 2005c; Temim et al. 2015) is the “crushing” of PWNe, which can be exerted by the inward-bound reverse front of the supernova shock wave. For SNRs developing asymmetrically, for instance due to an inhomogeneous surrounding medium (ISM), this crushing may result in considerabledistortion and displacement of the wind nebula. Put to a population-scoped graph (Fig. 6, left), it becomes clear that pulsar proper motions are insufficient to explain the large offsets observed, which may instead be due to reverse shock interaction being a dominant and frequent cause of pulsar-PWN offset in middle-aged systems (seealso de Jager & Djannati-Ataï 2009). Furthermore, the offset appears to relate to high efficiency (Fig. 9, right), suggesting that the PWN either gains energy and brightness through the process that causes the offset or that dense surroundings amplify both the IC luminosity and the offset between pulsar and wind nebula. While the evidence for this at present is not very strong, following up with expanded future studies is certainly worthwhile.

The expansion of PWNe with time (i.e. rising characteristic age and falling pulsar spin-down) could also be shown to be evident in the data. The fitted relation R ~ τc^0.55±0.23 suggests an average expansion coefficient in between those expected theoretically (1.2 and 0.3). The data set is not comprehensive enough to do a fit with two power laws, but appears to be consistent with the model. Notably, and in coherence with what was discussed already in Aharonian et al. (1997), this expansion is not so clear in X-rays, where the synchrotron emission always remains very local because it only traces the young particles in areas of high magnetic field relatively close to the pulsar. Most of the old objects (>30 kyr) in Kargaltsev et al. (2013) are therefore smaller than 1 pc in their bright core emission. On the other hand, in a limited sample of eight PWNe, Bamba et al. (2010) have reported the existence of an additional extended and expanding X-ray emission component, which might be the emission from the particles we see in TeV.

An interesting relation was found between the PWN surface brightness and pulsar Ė (Fig. 10, left). What stands out is not only the correlation itself, but also its relatively low scatter. This might either suggest that luminosity and extension are more correlated than reproduced in our model (such that a high-luminosity outlier is always balanced by an accordingly large extension), or it is an indication that the large scatter in all the other plots is dominated by the distance uncertainty, which is cancelled out in the surface brightness parameter. If this latter were true, it would mean that PWNe in fact evolve even more uniformly than suggested by our varied model.

The evolution trend of the photon index remains an open issue in this study. Neither the data nor the model are particularly clear about it for the young to middle-aged PWNe we investigated. A more sensitive data set – as expected from CTA (Acharya et al. 2013) – will reduce the uncertainties on spectral indices and reduce the selection bias by detecting more soft-spectrum PWNe.

Since both the H.E.S.S. Galactic Plane Survey and the ATNF pulsar database only cover a fraction of the Milky Way, depending on TeV and pulsar brightnesses, this study suffers from several selection biases discussed throughout the text. For TeV-bright, high-Ė, young pulsar systems (>10³⁶ erg s-1) we achieve a relatively good coverage, whereas for systems beyond some tens of kiloyears of age we likely miss many sources. In the plots discussing flux-related quantities, this is partly compensated by the inclusion of flux limits, allowing for statements that consider non-detections. For extension- and position-related quantities, however, we can only rely on the detected cases. It requires a full population synthesis study to judge whether some of the correlations are genuine or include side effects of other correlations or selection biases. Our plots and fits are meant to draw attention to where correlations may lurk and we encourage further work on this matter beyond the scope of this paper.

One presumably very influential parameter ignored in this study is the density of matter and background light at the position of each pulsar. It is likely due to such circumstances that Vela X, 3C 58, and CTA 1 are so faint, and N 157B (in the Large Magellanic Cloud) is so bright. In the scope of a population synthesis study, one could use a specific Milky Way model to “calibrate” the calorimetric objects that TeV pulsar wind nebulae are assumed to be.

On the modelling side, we are able to describe the trends and scatter of the TeV properties of the present PWN population with a relatively simple time-dependent modelling. Its 12 free parameters (7 of which were varied for the varied model) were well below the 4 × 19 observed parameters that the firmly identified PWNe provided¹⁰. It is remarkable that the adaptive parameters need to be varied in a fairly small range, compared to what one may fathom from the modelling literature¹¹ , while still producing sufficient scatter in the predicted observables (even excluding distance uncertainties and target photon densities as additional factors). Whether this indicates that the underlying variations of the individual PWN parameters are indeed small, or whether this is because the parameters are (anti-)correlated (see Sect. 4.3), cannot be clarified in this work. This requires a deeper physical model of the pulsars and possibly a multidimensional likelihood fit to correctly quantify all correlations and identify the true distributions of its parameters.

In the CTA era, most of the PWNe that will be detected in addition to the now assessed population will be middle-aged and old systems that are too faint or too extended to be detected with current instruments. Also, SKA (Taylor 2012) will enlarge the sample of pulsars detected in our Galaxy. To gain new insights from studying these systems, a solid and publicly available modelling code is needed that includes the difficult reverse shock interaction phase of a PWN in a reproducible way. This may help to understand the effect and influence of the amount of crushing and pulsar offset of the PWN, which is likely an influential factor of later PWN evolution.

On the analysis side, it would be beneficial to (i) improve the angular resolution and get to smaller scales of extension, (ii) find ways toreliably disentangle overlapping sources and their spectra, and (iii) aim for detecting objects larger than the IACT camera FOV. The latter is also of interest because pulsar systems in our Galactic neighbourhood, at few hundred parsecs from Earth, are considered plausible candidates to strongly contribute to the cosmic-ray electron and positron fluxes at Earth (e.g. Yin et al. 2013). The CTA cameras will provide us with a larger FOV (Acharya et al. 2013), which improves the capability of mapping out close-by PWNe. Detecting TeV objects even larger than that FOV will require better modelling and/or treatment of the cosmic-ray background event distribution and its systematics (e.g. Spengler 2015; Klepser 2012). In parallel, more generalised analysis packages with wholistic likelihood approaches (Knödlseder et al. 2013) might help us to unriddle sources that occult each other in the densely populated arms of the Galaxy.

Acknowledgments

We would like to thank Rolf Bühler for fruitful discussions on the PWN and Crab subjects. This work made extensive use of the ATNF Pulsar Catalogue (Manchester et al. 2005), the Python packages NumPy/SciPy (Dubois et al. 1996; Jones et al. 2001), and Matplotlib (Hunter 2007). The TeVCat Online Catalog for TeV astronomy (Horan & Wakely 2016) was a helpful resource in our literature research. The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the French Ministry for Research, the CNRS-IN2P3 and the Astroparticle Interdisciplinary Programme of the CNRS, the U.K. Science and Technology Facilities Council (STFC), the IPNP of the Charles University, the Czech Science Foundation, the Polish Ministry of Science and Higher Education, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the Innsbruck University, the Austrian Science Fund (FWF), and the Austrian Federal Ministry for Science, Research and Economy, and by the University of Adelaide and the Australian Research Council. We appreciate the excellent work of the technical support staff in Berlin, Durham, Hamburg, Heidelberg, Palaiseau, Paris, Saclay, and in Namibia in the construction and operation of the equipment. This work benefitted from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation.

Appendix A:

Appendix A Basic modelling of TeV pulsar wind nebulae

In the interpretation of the TeV characteristics of the PWN population described in this paper we have made use of a time-dependent one-zone model. It allows us to trace the evolution of the VHE lepton population, and hence the radiative output of a PWN, based on a few general assumptions. The specific model we adopt here was introduced by Mayer et al. (2012), but extended and improved for this work. Its essential traits are outlined in the following.

A.1 Spin-down evolution and energy conversion into energetic leptons

The model allows us to calculate the evolution of the non-thermal emission of a PWN in discrete time steps with an adaptive step size δt. In each step, the amount of spin-down energy converted into relativistic electrons and positrons is given as $$\Delta E_{\text{p}}(t)=\eta {\displaystyle {\int }_t^{t+\delta t}\text{ edot}}(t^{\prime })\text{d}{t}^{\prime },$$(A.1)

where the spin-down evolution Ė(t) of the pulsar (Eq. (3), p. 7) is characterised by the braking index n, the initial spin-down timescale τ₀, and the initial spin-down Ė₀. The lepton conversion efficiency η can be adjusted to account for additional cooling effects, but in this work is set to 1. This neglects the sub-percent fraction of magnetic energy release that should technically be missing in the particle outflow, but is negligible here.

While it is possible to transfer the dependency on τ₀ to one on P₀ using ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$(A.2)

in this work we take τ₀ as the free parameter.

A.2 Lepton injection spectrum

For the energy spectrum of leptons freshly injected into the nebula we assume the following power-law shape: ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$(A.3)

with a power-law index β. Φ₀ (t) can be calculated imposing ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$(A.4)

The lepton energies needed to deliver the relevant X-ray and gamma-ray energies cover a range of E_min to E_max. Varying the boundary energies essentially changes the number of particles contained in the IC-relevant energy range, and thus the efficiency, but does not fundamentally change the relative evolution of observables. A low-energy break in the injection spectrum is often applied in literature (e.g. Torres et al. 2014), but only impacts the lower ends of the emission spectra. We omit it here because it neither influences, nor is constrained by our data.

A.3 Cooling mechanisms

Cooling is approximated as $$\frac{\text{d}{N}_{\text{cooled}}}{\text{d}E}(E,t)=\frac{\text{d}N}{\text{d}E}(E,t-\delta t)\exp \left(-\frac{\delta t}{{\tau }_{\text{eff}}(E,t)}\right),$$(A.5)

with an effective cooling timescale $${\tau }_{\text{eff}}^{-1}={\tau }_{\text{syn}}^{-1}+{\tau }_{\text{esc}}^{-1}+{\tau }_{\text{ad}}^{-1},$$(A.6)

which comprises synchrotron, escape, and adiabatic losses. This strategy, as well as the expressions for the first two terms, are adopted from Zhang et al. (2008) ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$

Here, R(t) and B(t) describe the time evolution of the PWN radius and the magnetic field strength inside the PWN (cf. Sect. A.4 below). The timescale for adiabatic losses, ${\tau }_{\text{ad}}=-\frac{E}{\stackrel{.}{E_{\text{p}}}}$, is governed by the expansion of the nebula and can be calculated (following de Jager & Harding 1992) from $$\frac{\text{d}{E}_{\text{ad}}}{\text{d}t}=-\frac{E}{3}\nabla {\overrightarrow{v}}_{\perp }(R)={\dot{E}}_{\text{p}},$$(A.9)

with v_⊥(R) being the radial component of the particle velocity. In general, its divergence can be calculated to $$\text{begin}align*\nabla {\overrightarrow{v}}_{\perp }(R)\text{ \&}=\frac{1}{R^2}\frac{\partial (R^2{\overrightarrow{v}}_{\perp })}{\partial R}\text{nn\&}=\frac{1}{R{(t)}^2}\frac{\partial (R{(t)}^2{\overrightarrow{v}}_{\perp }(t))}{\partial t}\frac{\partial t}{\partial R}\text{nnend}align*$$

making use of the radial evolution function R(t) given in the next section. In addition to the above formulation we take into account losses originating from inverse Compton (IC) emission. This is achieved by subtracting the IC emissivity in each time step dependent on the electron energy (see Appendix A.5 for further details on the IC emissivity).

A.4 Dynamical evolution

In order to take into account that the growth rate of a PWN strongly depends on its evolutionary state, the model builds on analytical studies of the development of PWNe inside their SNR environment (e.g. Chevalier 1977; Reynolds & Chevalier 1984). The time evolution implemented in the model comprises three¹² distinct phases, which define the expansion behaviour of the PWN according to the age of the system in terms of the spin-down timescale τ₀ and the reverse-shock interaction time t_rs. Usually, thereverse-shock passage and the subsequent reverberations are expected to occur at a time t_rs > τ₀. For this case, the following relations have been derived in the aforementioned works: $$R(t)\propto \text{nn begin{cases}}{t}^{6/5}\text{ \&{for begininlinestripnssi752}\$\text{t}⩽{\tau }_{\text{0}}$$(A.12)

In the (supposably much less common) opposite case, t_rs < τ₀, the time evolution of R_PWN is modified to $$R(t)\propto \text{nn begin{cases}}{t}^{6/5}\text{ \&{for begininlinestripnssi758}\$\text{t}⩽{\text{t}}_{\text{rs}}$$(A.13)

As a simplification, the crushing of the PWN by the SNR reverse shock is not modelled here. Such crushing presumably reduces the radius between free expansion and reverse shock interaction phase.

The magnetic field evolution is adapted from Zhang et al. (2008) as $$B(t)=\frac{B_0}{1+{\left(\frac{t}{{\tau }_0}\right)}^{\alpha }}+B_{\text{ISM}},$$(A.14)

assuming a constant and homogeneous ISM contribution of 3 μG , and adopting an index of α = 0.6 in order to satisfy the conservation of magnetic flux.

A.5 Time-dependent lepton energy distribution and radiative processes

The framework laid out in the previous sections allows us to calculate the energy distribution of the leptons contained in the PWN at any given time. More specifically, the number of leptons with energy E residing in the nebula at a time t + δt is determined by the balance of freshly injected leptons and those cooled out of the respective energy interval, $$\frac{\text{d}N}{\text{d}E}(E,t+\delta t)=\frac{\text{d}{N}_{\text{cooled}}}{\text{d}E}(E,t)+\frac{\text{d}{N}_{\text{inj}}}{\text{d}E}(E,t+\delta t).$$(A.15)

The iterative evaluation of Eq. (A.15) then yields the lepton energy distribution as a function of time. The time binning is adjusted adaptively to guarantee high precision at a still reasonable computation cost (see Mayer 2010b, Sect. 5.2.2. and Fig. 5.8 for details on this).

From the lepton distribution, the photon population arising from synchrotron emission and inverse Compton scattering as the most important processes can be obtained. The physics of these processes is described in the comprehensive review article by Blumenthal & Gould (1970), which we follow in the implementation of the radiation mechanisms within our model. The target photon fields required as an input for IC scattering are CMB, starlight, and infrared photons. While the uniform CMB component is modelled as a black-body spectrum with an energy density of 0.26 eV cm⁻³ and temperature of 2.7 K, the starlight and infrared components can be adopted from the Galprop code (Porter & Strong 2005). In order to derive a representative radiation field composition for the baseline model, the Galprop fields at the positions of all firmly identified PWNe were averaged, using the mean temperature and energy densities as input for the respective black-body spectra. Following this set-up, the energy densities of the starlight and infrared fields are 1.92 eV cm⁻³ and 1.19 eV cm⁻³, respectively. The temperatures at the spectral peaks are 107 K for the infrared and 7906 K for the starlight field component.

A.6 Results of the time-dependent modelling

In summary, the model takes the parameters listed in Table A.1. The table contains two compilations of parameters: the first one states the values used for the baseline model, which is depicted as a black line throughout the population plots in this paper; the second one gives the ranges of parameters we used to mimic the intrinsic spread of the PWN properties. The PWN evolution implied by our baseline model is listed in Table A.2.

The considerations that went into the choice of the model parameters and ranges are the following:

• We want to mimic a typical PWN in a typical (dense spiral arm) surrounding. For this reason, we do not give objects like Vela X, 3C 58, or CTA 1 too much consideration in the adjustment of the parameters. This can make the model differ from the fit results, which take all objects into account.

• n: the braking index defines the slope of the pulsar trajectory in Fig. 2, which has to be ~3–4 to matchthe measured pulsar population. The theoretical expectation is that n = 3 if the energy loss is dominated by magnetic dipole radiation, whereas a spin-down dominated by gravitational radiation leads to a longer energy release through n = 5 (e.g. Yue et al. 2007). By contrast, the few direct measurements of braking indices presently available lie in the range of 0.9–2.8 (for a compilation see Magalhaes et al. 2012), indicating a much faster spin-down decay. In this study, we set it to the canonical n = 3.

• τ₀, Ė₀: these parameters define the starting point of the pulsar trajectory on Fig. 2 and the total energy budget of the pulsar (see Appendix B and Fig. 2). With the chosen combination and the canonical n = 3, the pulsar energy outflow evolves along the path where ATNF pulsars are actually found and starts with a total energy of Ė₀τ₀ = 3.1 × 1049 erg.

• B₀: the initial B-field is chosen such that, using Eq. (A.14) for its decay, Crab-like young PWNe have (present) fields BPWN ~ 100 μG, while older objects at some point arrive at few tens of μG or less. This is consistent with the ranges found in other modelling works, such as Torres et al. (2014) and Zhang et al. (2008), in which this scale is set with the goal of producing realistic X-ray luminosities.

• t_rs, R₅: these parameters determine the dynamical evolution and are set such that the PWN extension trajectory evolves roughly through the middle of the firmly identified PWNe. They also have strong influence on the surface brightness plot Fig. 10 (left), which interlinks them with the luminosity related parameters.

• BISM: set to the canonical 3 μG.

• η: the lepton efficiency can account for a substantial fraction of energy going into magnetic fields or hadron acceleration, neither of which we assume to be large. Hence, we set η = 1.

• α: set to 0.6 in order to satisfy the conservation of magnetic flux.

• β: an injection index of 2 is a typical value found to lead to good agreement with observed spectral indices here and in other works. The variation we induce produces a realistic variation of gamma-ray photon indices.

• E_min, E_max: these energy bounds mainly determine the ranges of the synchrotron and IC photon spectra, and therefore also the amount of photons produced specifically in the 1–10 TeV band considered here. They are not constrained by our plots beyond this efficiency variation they can provoke.

The parameter set of the baseline model was also used to construct the spectral energy distributions (SEDs) shown in Fig. A.1. These sample SEDs illustrate the time evolution of the radiative output of a generic PWN according to the presented model.

Fig. A.1 Modelled spectral energy distribution (SED) of a generic PWN with parameters according to baseline model given in Table A.1. See Appendix sectionlinking A.7 for caveats of the SEDs. Left: time evolution of the SED, ranging from 1 kyr to 200 kyr. Right: decomposition of the SED of a middle-aged PWN (10 kyr; black dashed curve) into contributions by leptons from various injection epochs (coloured lines). The grey-shaded bands indicate the energy range of 1–10 TeV explored in this paper.

The set of model curves in Fig. A.1 (left) traces the various evolutionary stages of the energy flux at PWN ages ranging between 0.5 kyr and 150 kyr, calculated in equidistant steps on a logarithmic timescale. Even though both the synchrotron and IC contributions obviously undergo significant development with increasing age of the system, the decline of the synchrotron energy flux (due to its strong dependence on the decaying magnetic field strength) is more pronounced than that of the IC component.

Figure A.1 (right) depicts the SED of a generic middle-aged PWN decomposed into numerous contributions from individual epochs. The dominance of the very youngest leptons in producing the synchrotron component (most notably the X-ray part) is manifest in this plot. By contrast, accumulated leptons from various ages contribute to the IC radiation, in particular in the TeV energy band.

A.7 Caveats

As already emphasised in Sect. 4.3 and Appendix A.2, the aim of this model is to serve for the interpretation of the TeV data we have. Spectral breaks, potential reverberation compressions, and other aspects that cannot be judged with the present data are therefore omitted on purpose. The multi-wavelength spectra it predicts, though found in the right order of magnitude, may therefore not be very accurate at energies other than the TeV regime.

Another caveat to note is the correlation of parameters. We vary only 7 of the 12 parameters (the target photon field could additionally be regarded as a 13th parameter), but the variations in the model can of course also be achieved by varying more of the parameters by a smaller magnitude. A variation of Ė is for instance indistinguishable, from the point of view of the TeV properties, from a variation of lepton efficiency. So the variation solution we found leads to a sensible range in predicted observable ranges, but is not unique. Similarly, a correlation of two parameters can mean that larger variations are possible, such as in the example described in Sect. 4.3.

Appendix B Derivation of basic formulae around the relation of $\stackrel{.}{\text{textsfslshape}E}$ and τc

Since the following relations are relatively fundamental to the energy input evolution of PWNe, but still rather hard to find in recent literature, we briefly want to wrap up what Eqs. (3)–(5) and (A.2) are derived from.

As pointed out by Gunn & Ostriker (1969), the energy loss rate of a rotating magnetic dipole depends on the angular velocity Ω as ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$(B.1)

Since the angular momentum loss rate is ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$(B.2)

it follows that the velocity loss rate is ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$(B.3)

where I is the neutron star moment of inertia. To generalise this relation for the non-dipole case, the index 3 is replaced by the braking index n, $$\dot{\Omega }=\frac{\dot{J}}{I}=-k{\Omega }^n,$$(B.4)

which turns Eq. (B.1) into ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$(B.5)

The general solution of the differential equation (Eq. (B.4)) can be written as $$\Omega (t)={\Omega }_0{\left(1+\frac{t}{{\tau }_0}\right)}^{-\frac{1}{n-1}}.$$(B.6)

Using Eq. (B.5), and P = 2π∕Ω, and differentiating P one obtains ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$ ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$ ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$

The canonical formulae to calculate Ė and τc from P and Ṗ then yield ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$

(cf. the notation in Gaensler & Slane 2006, Eqs. (5) and (6)). Note that neither of these expressions relies on the dipole hypothesis of n = 3. At the birth of the pulsar, t = 0, τc is τ₀ (n − 1)∕2 and increases steadily.

For the relation of Ṗ and P, Eqs. (B.8) and (B.9) furthermore imply $$\text{pdot}(P)=\frac{P_0}{{\tau }_0}\frac{1}{n-1}{\left(\frac{P}{P_0}\right)}^{2-n},$$(B.12)

which can be taken to discuss plausible braking indices directly from Fig. 2 (right).

In order to see what happens if Ė is plotted against τc, one has to resolve the dependency on t to arrive at $$\text{edot}(\text{age})={\text{edot}}_0{\left(\frac{2}{n-1}\frac{\text{age}}{{\tau }_0}\right)}^{-\frac{n+1}{n-1}}.$$(B.13)

Clearly, the evolution curve of a pulsar on the Ė -τc diagram starts at a point [τ₀ (n − 1)∕2, Ė₀], which depends on τ₀, n, and Ė₀ , but the slope of the power law is only dictated by the braking index n. This index is not predetermined to be 3 by the way τc is constructed.

Assuming that Eq. (B.7) describes the energy outflow of the pulsar throughout its lifetime, one can calculate the energy deposited up to a certain time as follows: ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$

For t →∞, Ė (t) τc(t) vanishes, so the first term represents the total energy budget that is emitted and, using Eq. (B.7), can be made equivalent to ${\varphi }_0\times {(E/E_0)}^{-\Gamma }$, the total rotational energy of the pulsar. Unfortunately, n, Ė₀ , and τ₀ are three unknown initial properties of the pulsar, so it cannot be measured. Unlike that, the second term Ė τc, which represents the present budget of rotational energy, can be calculated from the measured P and Ṗ . The ordinary (low-aged) pulsar with the maximum present budget of energy is PSR J0537−6910 in N 157B, with 7.6 × 1049 erg, which is a lower limit to the maximal initial rotational energies that can be reached.

References

All Tables

All Figures

Fig. 1 Histograms of spatial separation between PMPS pulsars and TeV source components from the HGPSC list. In the high-Ė pulsar sample (left), a clear correlation is seen as a peak at small squared angular distances, whereas the low-Ė associations,if present, are not significant beyond the expected rate of chance coincidences (right). The angular separation cut of θ < 0.5° applied in the preselection of PWN candidates (Sect. 3.2) is indicated by a dashed vertical line in the left panel.

In the text

Fig. 2 Left: spin-down power Ė and characteristic age τc of pulsars with a firmly identified PWN, candidate PWN, and without TeV counterpart (grey dots). The black line and shaded band show the injection evolution of the modelling used in this paper. The dashed lines indicate lines of constant total remaining energy Ė τ; see Appendix B. Hence a model curve that starts at Ė0τ0 = 1049 erg represents a pulsar with total initial rotational energy of 1049 erg. Since both Ė and τc depend on P and Ṗ , the axes in this plot do not represent independent quantities. Right: same data, shown in the commonly used view, using the independently measured P and Ṗ .

In the text

Fig. 3 Schematic of the Milky Way and its spiral arms, along with firmly identified PWNe, candidates, and energetic pulsars (Ė > 1035 erg s-1) without detected TeV wind nebula. The yellow and blue curves outline the sensitivity horizon of the HGPS for point-like sources with an integrated gamma-ray luminosity (1–10 TeV) of 1% and 10% of the Crab luminosity, respectively (see H.E.S.S. Collaboration 2018, for details).

In the text

	Fig. 4 Top: PWN extension occurrences over distance from Earth, in comparison to the band of extensions that can be expectedto be identified in the HGPS analysis chain. Bottom: distribution of known distances of energetic pulsars (Ė > 1035 erg s-1).
In the text

Fig. 5 Left: PWN extension evolution with time, in comparison to the modelling considered in this work. Right: PWN extension evolution with Ė, as fitted in the RPWN(Ė) column of Table 6 for pulsar wind nebulae with Ė > 1036 erg s-1 (see Sect. 5.2.1). The shaded range shows the fit range and standard deviation σlg R . 1 dex refers to an order of magnitude and is the unit of the logspace defined σlg Y . For clarity, this plot excludes PWN candidates and divides the sample into nearby and far pulsar wind nebulae to illustrate the potential selection or reconstruction bias (see text). The dot-dashed and dotted lines indicate the systematic minimum of 0.03° and the maximum measured extension in the HGPS of 0.6°, respectively, which are both projected to the average PWN distance of 5.1 kpc.

In the text

Fig. 6 Left: spatial offset of pulsar and TeV wind nebula as a function of pulsar characteristic age. The dashed and dotted lines indicate the systematic minimum of about 0.015° and the association criterion of 0.5°, which are both projected to the average PWN distance of 5.1 kpc. The solid line shows the offset one would expect from pulsar motion only (assuming a large pulsar velocity of 500 km s−1 ). Right: time evolution of the containment ratio. Since the pulsar motion can be assumed to be constant and the expansion decelerates, one expects the containment fraction to increase and eventually pass unity after some tens of kiloyears. The dotted horizontal line shows the rating criterion (offset/extension <1.5) applied in the post-selection of candidates (Sect. 6).

In the text

	Fig. 7 Left: relation of TeV luminosity and pulsar Ė. Right: same as Fig. 5, but for luminosity and including the upper limits from the left panel and 11 upper limits calculated for pulsars of candidate PWNe (see text).
In the text

	Fig. 8 Left: evolution of TeV luminosity with characteristic age. Right: distribution of 1–10 TeV luminosity for PWNe and PWN candidates treated in this work.
In the text

	Fig. 9 Left: evolution of TeV efficiency (L1−10 TeV∕Ė) with pulsar characteristic age. Right: TeV efficiency as a function of pulsar offset, plotted for pulsars of different age groups. High-offset systems tend to be more TeV-efficient than low-offset systems.
In the text

	Fig. 10 Left: dependence of surface brightness S (see text) on pulsar spin-down Ė. Right: TeV photon index over pulsar spin-down Ė for all detected PWNe and candidates.
In the text

	Fig. 11 Common logarithmic residuals of rating criteria 2–4, using the standard deviations σlg Y explained in Sect. 6. Left: extension with respect to the model shown in Fig. 5 (left). Middle: same for luminosity (Fig. 7, left). Right: same for surface brightness (Fig. 10, left). In all cases, limits are shown separately as outlined histograms.
In the text

Fig. A.1 Modelled spectral energy distribution (SED) of a generic PWN with parameters according to baseline model given in Table A.1. See Appendix sectionlinking A.7 for caveats of the SEDs. Left: time evolution of the SED, ranging from 1 kyr to 200 kyr. Right: decomposition of the SED of a middle-aged PWN (10 kyr; black dashed curve) into contributions by leptons from various injection epochs (coloured lines). The grey-shaded bands indicate the energy range of 1–10 TeV explored in this paper.

In the text