# The Population of TeV Pulsar Wind Nebulae in the H.E.S.S. Galactic Plane Survey

H. Abdalla, P. Aubert, J. Chevalier, A. Fiasson, F. Krayzel, G. Lamanna, J.P.

Lees, G. Maurin, S. Rosier-Lees, D.A. Sanchez, et al.

► To cite this version:

H. Abdalla, P. Aubert, J. Chevalier, A. Fiasson, F. Krayzel, et al. The Population of TeV Pulsar Wind Nebulae in the H.E.S.S. Galactic Plane Survey. Astronomy and Astrophysics - A&A, 2018, 612, pp.A2. 10.1051/0004-6361/201629377. in2p3-01480223

HAL Id: in2p3-01480223

https://hal.in2p3.fr/in2p3-01480223

Submitted on 24 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The population of TeV pulsar wind nebulae in the H.E.S.S. Galactic Plane Survey

H.E.S.S. Collaboration*, H. Abdalla¹, A. Abramowski², F. Aharonian³,⁴,⁵, F. Ait Benkhali³, A. G. Akhperjanian⁶,⁵,†, T. Andersson¹⁰,
E. O. Angüner⁷, M. Arrieta¹⁵, P. Aubert²⁴, M. Backes⁸, A. Balzer⁹, M. Barnard¹, Y. Becherini¹⁰, J. Becker Tjus¹¹, D. Berge¹², S. Bernhard¹³,
K. Bernlöhr³, R. Blackwell¹⁴, M. Böttcher¹, C. Boisson¹⁵, J. Bolmont¹⁶, P. Bordas³, J. Bregeon¹⁷, F. Brun²⁶, P. Brun¹⁸, M. Bryan⁹, T. Bulik¹⁹,
M. Capasso²⁹, J. Carr²⁰, S. Carrigan³,⁴⁵, S. Casanova²¹,³, M. Cerruti¹⁶, N. Chakraborty³, R. Chalme-Calvet¹⁶, R. C. G. Chaves¹⁷,²², A. Chen²³,
J. Chevalier²⁴, M. Chrétien¹⁶, S. Colafrancesco²³, G. Cologna²⁵, B. Condon²⁶, J. Conrad²⁷,²⁸, C. Couturier¹⁶, Y. Cui²⁹, I. D. Davids¹,⁸,
B. Degrange³⁰, C. Deil³, J. Devin¹⁷, P. deWilt¹⁴, L. Dirson², A. Djannati-Ataï³¹, W. Domainko³, A. Donath³, L. O’C. Drury⁴, G. Dubus³²,
K. Dutson³³, J. Dyks³⁴, T. Edwards³, K. Egberts³⁵, P. Eger³, J.-P. Ernenwein²⁰, S. Eschbach³⁶, C. Farnier²⁷,¹⁰, S. Fegan³⁰, M. V. Fernandes²,
A. Fiasson²⁴, G. Fontaine³⁰, A. Förster³, S. Funk³⁶, M. Füßling³⁷, S. Gabici³¹, M. Gajdus⁷, Y. A. Gallant¹⁷,*, T. Garrigoux¹, G. Giavitto³⁷,
B. Giebels³⁰, J. F. Glicenstein¹⁸, D. Gottschall²⁹, A. Goyal³⁸, M.-H. Grondin²⁶, D. Hadasch¹³, J. Hahn³, M. Haupt³⁷, J. Hawkes¹⁴,
G. Heinzelmann², G. Henri³², G. Hermann³, O. Hervet¹⁵,⁴⁴, A. Hillert³, J. A. Hinton³, W. Hofmann³, C. Hoischen³⁵, M. Holler³⁰, D. Horns²,
A. Ivascenko¹, A. Jacholkowska¹⁶, M. Jamrozy³⁸, M. Janiak³⁴, D. Jankowsky³⁶, F. Jankowsky²⁵, M. Jingo²⁳, T. Jogler³⁶, L. Jouvin³¹,
I. Jung-Richardt³⁶, M. A. Kastendieck², K. Katarzyński³⁹, U. Katz³⁶, D. Kerszberg¹⁶, B. Khélifi³¹, M. Kieffer¹⁶, J. King³, S. Klepser³⁷,*,
D. Klochkov²⁹, W. Kluźniak³⁴, D. Kolitzus¹³, Nu. Komin²³, K. Kosack¹⁸, S. Krakau¹¹, M. Kraus³⁶, F. Krayzel²⁴, P. P. Krüger¹, H. Laffon²⁶,
G. Lamanna²⁴, J. Lau¹⁴, J.-P. Lees²⁴, J. Lefaucheur¹⁵, V. Lefranc¹⁸, A. Lemière³¹, M. Lemoine-Goumand²⁶, J.-P. Lenain¹⁶, E. Leser³⁵, T. Lohse⁷,
M. Lorentz¹⁸, R. Liu³, R. López-Coto³, I. Lypova³⁷, V. Marandon³, A. Marcowith¹⁷, C. Mariaud³⁰, R. Marx³, G. Maurin²⁴, N. Maxted¹⁴,
M. Mayer⁷,*, P. J. Meintjes⁴⁰, M. Meyer²⁷, A. M. W. Mitchell³, R. Moderski³⁴, M. Mohamed²⁵, L. Mohrmann³⁶, K. Morå²⁷, E. Moulin¹⁸,
T. Murach⁷, M. de Naurois³⁰, F. Niederwanger¹³, J. Niemiec²¹, L. Oakes⁷, P. O’Brien³³, H. Odaka³, S. Öttl¹³, S. Ohm³⁷, E. de Oña Wilhelmi³,⁴⁶,
M. Ostrowski³⁸, I. Oya³⁷, M. Padovani¹⁷, M. Panter³, R. D. Parsons³, M. Paz Arribas⁷, N. W. Pekeur¹, G. Pelletier³², C. Perennes¹⁶,
P.-O. Petrucci³², B. Peyaud¹⁸, S. Pita³¹, H. Poon³, D. Prokhorov¹⁰, H. Prokoph¹⁰, G. Pühlhofer²⁹, M. Punch³¹,¹⁰, A. Quirrenbach²⁵, S. Raab³⁶,
A. Reimer¹³, O. Reimer¹³, M. Renaud¹⁷, R. de los Reyes³, F. Rieger³,⁴¹, C. Romoli⁴, S. Rosier-Lees²⁴, G. Rowell¹⁴, B. Rudak³⁴, C. B. Rulten¹⁵,
V. Sahakian⁶,⁵, D. Salek⁴², D. A. Sanchez²⁴, A. Santangelo²⁹, M. Sasaki²⁹, R. Schlickeiser¹¹, F. Schüssler¹⁸, A. Schulz³⁷, U. Schwanke⁷,
S. Schwemmer²⁵, M. Settimo¹⁶, A. S. Seyffert¹, N. Shafi²³, I. Shilon³⁶, R. Simoni⁹, H. Sol¹⁵, F. Spanier¹, G. Spengler²⁷, F. Spies², Ł. Stawarz³⁸,
R. Steenkamp⁸, C. Stegmann³⁵,³⁷, F. Stinzing³⁶,†, K. Stycz³⁷, I. Sushch¹, J.-P. Tavernet¹⁶, T. Tavernier³¹, A. M. Taylor⁴, R. Terrier³¹, L. Tibaldo³,
D. Tiziani³⁶, M. Tłuczykont², C. Trichard²⁰, R. Tuffs³, Y. Uchiyama⁴³, K. Valerius³⁶,⁴⁷,*, D. J. van der Walt¹, C. van Eldik³⁶, B. van Soelen⁴⁰,
G. Vasileiadis¹⁷, J. Veh³⁶, C. Venter¹, A. Viana³, P. Vincent¹⁶, J. Vink⁹, F. Voisin¹⁴, H. J. Völk³, T. Vuillaume²⁴, Z. Wadiasingh¹, S. J. Wagner²⁵,
P. Wagner⁷, R. M. Wagner²⁷, R. White³, A. Wierzcholska²¹, P. Willmann³⁶, A. Wörnlein³⁶, D. Wouters¹⁸, R. Yang³, V. Zabalza³³, D. Zaborov³⁰,
M. Zacharias²⁵, A. A. Zdziarski³⁴, A. Zech¹⁵, F. Zefi³⁰, A. Ziegler³⁶, and N. Żywucka³⁸

*(Affiliations can be found after the references)*

Received 22 July 2016 / Accepted 24 February 2017

## ABSTRACT

The nine-year H.E.S.S. Galactic Plane Survey (HGPS) has yielded the most uniform observation scan of the inner Milky Way in the TeV gamma-ray band to date. The sky maps and source catalogue of the HGPS allow for a systematic study of the population of TeV pulsar wind nebulae found throughout the last decade. To investigate the nature and evolution of pulsar wind nebulae, for the first time we also present several upper limits for regions around pulsars without a detected TeV wind nebula. Our data exhibit a correlation of TeV surface brightness with pulsar spin-down power $\dot{E}$. This seems to be caused both by an increase of extension with decreasing $\dot{E}$, and hence with time, compatible with a power law $R_{\text{PWN}}(\dot{E}) \sim \dot{E}^{-0.65 \pm 0.20}$, and by a mild decrease of TeV gamma-ray luminosity with decreasing $\dot{E}$, compatible with $L_{1-10\text{TeV}} \sim \dot{E}^{0.59 \pm 0.21}$. We also find that the offsets of pulsars with respect to the wind nebula centre with ages around 10 kyr are frequently larger than can be plausibly explained by pulsar proper motion and could be due to an asymmetric environment. In the present data, it seems that a large pulsar offset is correlated with a high apparent TeV efficiency $L_{1-10\text{TeV}}/\dot{E}$. In addition to 14 HGPS sources considered firmly identified pulsar wind nebulae and 5 additional pulsar wind nebulae taken from literature, we find 10 HGPS sources that are likely TeV pulsar wind nebula candidates. Using a model that subsumes the present common understanding of the very high-energy radiative evolution of pulsar wind nebulae, we find that the trends and variations of the TeV observables and limits can be reproduced to a good level, drawing a consistent picture of present-day TeV data and theory.

**Key words.** gamma rays: general – catalogs – surveys – ISM: supernova remnants – pulsars: general

## 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. Pulsars can live for up to $10^{5-6}$ kyr, but their magnetic and particle outflow

<footnote>* Corresponding authors: H.E.S.S. Collaboration, e-mail: contact.hess@hess-experiment.eu</footnote>

<footnote>† Deceased.</footnote>
is decreasing steadily. Therefore, most of the observed PWNe
are associated with pulsars that are less than a few 100 kyr old
(Roberts 2004).

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 emit-
ted as pulsed electromagnetic radiation (Abdo et al. 2013). Most
of the pulsar outflow consists of high-energy particles and mag-
netic fields that feed into the growing PWN plasma. This plasma
is dynamically inferior to the ~10⁵<fnref n="1" /> erg carried away by the su-
pernova 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 emis-
sion component, they carry important information that the syn-
chrotron emission, albeit much higher in flux and energy trans-
port, 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 neb-
ulae than those observed in the radio or X-ray regimes. So in
general, the IC image gives a more accurate and complete pic-
ture 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 Galac-
tic sources have meanwhile been associated with TeV pulsar
wind nebulae. Most of these objects are situated in the in-
ner 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 VERI-
TAS (e.g. Aliu et al. 2013), and arrays of air shower detec-
tors such as MILAGRO (Abdo et al. 2009), have also observed
PWNe and contributed very valuable case studies, mostly of sys-
tems 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 re-
lease 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 ob-
jects, and their duty cycle is restricted to dark, cloudless nights.

A systematic search with the *Fermi* Large Area Tele-
scope for GeV pulsar wind nebulae in the vicinity of TeV-
detected sources (Acero et al. 2013) yielded 5 firmly identi-
fied 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 $\ell = 250^\circ$ to $65^\circ$ and latitudes $|b| \le 3.5^\circ$. We
undertake a census of all the firmly identified PWNe detected
with H.E.S.S. and other IACTs, and for the first time comple-
ment this sample with HGPS flux upper limits for all covered
pulsar locations without a corresponding TeV detection. This al-
lows 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<fnref n="1" /> (HGPS; H.E.S.S. Collaboration 2018)
and the ATNF pulsar catalogue<fnref n="2" /> (Manchester et al. 2005, ver-
sion 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<fnref n="3" /> 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 un-
certainty, however, is not very well defined and can be as large
as a factor of 2. For the few cases in which pulsar distance esti-
mations 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 con-
sidered as firmly identified PWNe, we use the identification cri-
teria discussed in the HGPS paper and take as a starting point
the list of all 12 identified PWNe and the 8 identified compos-
ite SNRs (H.E.S.S. Collaboration 2018, Table 3). Most PWNe
in the HGPS are identified by positional and/or morphologi-
cal coincidence with a PWN identified in other wavelengths,
or by their specific (mostly energy dependent) TeV morphol-
ogy. Our selection for this paper also requires that the cor-
responding 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

<footnote n="1">¹ http://www.mpi-hd.mpg.de/HESS/hGPS</footnote>

<footnote n="2">² http://www.atnf.csiro.au/research/pulsar/psrcat</footnote>

<footnote n="3">³ The difference between the two is that the ATNF pulsar catalogue is a full listing of different surveys and targeted observations, including, for instance, Fermi-LAT detected gamma-ray pulsars, whereas the PMPS is a comparably uniform survey of one particular radio instrument and hence it is less prone to observational biases.</footnote>
<table><caption>**Table 1.** HGPS sources considered as firmly identified pulsar wind nebulae in this paper.</caption><thead><tr><th>HGPS name</th><th>ATNF name</th><th>Canonical name</th><th>lg <i>E&#x0307;</i></th><th>&tau;<sub>c</sub><br>(kyr)</th><th>d<br>(kpc)</th><th>PSR offset<br>(pc)</th><th>&Gamma;</th><th><i>R</i><sub>PWN</sub><br>(pc)</th><th><i>L</i><sub>1-10 TeV</sub><br>(10<sup>33</sup> erg s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J1813-178<sup>1</sup></td><td>J1813-1749</td><td></td><td>37.75</td><td>5.60</td><td>4.70</td><td>&lt;2</td><td>2.07 &pm; 0.05</td><td>4.0 &pm; 0.3</td><td>19.0 &pm; 1.5</td></tr><tr><td>J1833-105</td><td>J1833-1034</td><td>G21.5-0.9<sup>2</sup></td><td>37.53</td><td>4.85</td><td>4.10</td><td>&lt;2</td><td>2.42 &pm; 0.19</td><td>&lt;4</td><td>2.6 &pm; 0.5</td></tr><tr><td>J1514-591</td><td>B1509-58</td><td>MSH 15-52<sup>3</sup></td><td>37.23</td><td>1.56</td><td>4.40</td><td>&lt;4</td><td>2.26 &pm; 0.03</td><td>11.1 &pm; 2.0</td><td>52.1 &pm; 1.8</td></tr><tr><td>J1930+188</td><td>J1930+1852</td><td>G54.1+0.3<sup>4</sup></td><td>37.08</td><td>2.89</td><td>7.00</td><td>&lt;10</td><td>2.6 &pm; 0.3</td><td>&lt;9</td><td>5.5 &pm; 1.8</td></tr><tr><td>J1420-607</td><td>J1420-6048</td><td>Kookaburra (K2)<sup>5</sup></td><td>37.00</td><td>13.0</td><td>5.61</td><td>5.1 &pm; 1.2</td><td>2.20 &pm; 0.05</td><td>7.9 &pm; 0.6</td><td>44 &pm; 3</td></tr><tr><td>J1849-000</td><td>J1849-0001</td><td>IGR J18490-0000<sup>6</sup></td><td>36.99</td><td>42.9</td><td>7.00</td><td>&lt;10</td><td>1.97 &pm; 0.09</td><td>11.0 &pm; 1.9</td><td>12 &pm; 2</td></tr><tr><td>J1846-029</td><td>J1846-0258</td><td>Kes 75<sup>2</sup></td><td>36.91</td><td>0.728</td><td>5.80</td><td>&lt;2</td><td>2.41 &pm; 0.09</td><td>&lt;3</td><td>6.0 &pm; 0.7</td></tr><tr><td>J0835-455</td><td>B0833-45</td><td>Vela X<sup>7</sup></td><td>36.84</td><td>11.3</td><td>0.280</td><td>2.37 &pm; 0.18</td><td>1.89 &pm; 0.03</td><td>2.9 &pm; 0.3</td><td>0.83 &pm; 0.11*</td></tr><tr><td>J1837-069<sup>8</sup></td><td>J1838-0655</td><td></td><td>36.74</td><td>22.7</td><td>6.60</td><td>17 &pm; 3</td><td>2.54 &pm; 0.04</td><td>41 &pm; 4</td><td>204 &pm; 8</td></tr><tr><td>J1418-609</td><td>J1418-6058</td><td>Kookaburra (Rabbit)<sup>5</sup></td><td>36.69</td><td>10.3</td><td>5.00</td><td>7.3 &pm; 1.5</td><td>2.26 &pm; 0.05</td><td>9.4 &pm; 0.9</td><td>31 &pm; 3</td></tr><tr><td>J1356-645<sup>9</sup></td><td>J1357-6429</td><td></td><td>36.49</td><td>7.31</td><td>2.50</td><td>5.5 &pm; 1.4</td><td>2.20 &pm; 0.08</td><td>10.1 &pm; 0.9</td><td>14.7 &pm; 1.4</td></tr><tr><td>J1825-137<sup>10</sup></td><td>B1823-13</td><td></td><td>36.45</td><td>21.4</td><td>3.93</td><td>33 &pm; 6</td><td>2.38 &pm; 0.03</td><td>32 &pm; 2</td><td>116 &pm; 4</td></tr><tr><td>J1119-614</td><td>J1119-6127</td><td>G292.2-0.5<sup>11</sup></td><td>36.36</td><td>1.61</td><td>8.40</td><td>&lt;11</td><td>2.64 &pm; 0.12</td><td>14 &pm; 2</td><td>23 &pm; 4</td></tr><tr><td>J1303-631<sup>12</sup></td><td>J1301-6305</td><td></td><td>36.23</td><td>11.0</td><td>6.65</td><td>20.5 &pm; 1.8</td><td>2.33 &pm; 0.02</td><td>20.6 &pm; 1.7</td><td>96 &pm; 5</td></tr></tbody></table>

**Notes.** The sources are sorted by decreasing $\dot{E}$. $\lg \dot{E}$ stands for $\log_{10}(\dot{E}/\text{erg s}^{-1})$, $\tau_c$ is the pulsar characteristic age, $d$ is the pulsar distance, $R_{\text{PWN}}$ is the 1-sigma Gaussian extension and $L_{1-10\text{TeV}}$ is the TeV luminosity. The pulsar distances are printed uniformly here, but their uncertainties might often be larger or not available; see ATNF Catalogue references (http://www.atnf.csiro.au/people/pulsar/psrcat/) for detailed information. The limits are 2-sigma limits (see Sect. 2.3). (*) The luminosity of Vela X is calculated as described in Sect. 2.3.

**References.** Previous publications on these sources: (1) Funk et al. (2007); (2) Djannati-Ataï et al. (2008); (3) Aharonian et al. (2005b); (4) Acciari et al. (2010); (5) Aharonian et al. (2006e); (6) Terrier et al. (2008); (7) Aharonian et al. (2006a); (8) Gotthelf & Halpern (2008); (9) Renaud et al. (2008); (10) Aharonian et al. (2005c); (11) Acero et al. (2013); (12) H.E.S.S. Collaboration (2012a).

<table><caption>**Table 2.** List of ATNF pulsar distance estimates that were modified.</caption><thead><tr><th>PSR</th><th>Distance<br>(kpc)</th><th>Method/adjacent object</th><th>Reference</th></tr></thead><tbody><tr><td>J0205+6449 (3C 58)</td><td>2.0</td><td>H1 absorption</td><td>Kothes (2013)</td></tr><tr><td>J1023-5746 (Westerlund 2)</td><td>8.0</td><td>Westerlund 2 open cluster</td><td>Rauw et al. (2007)</td></tr><tr><td>J1418-6058 (Rabbit)</td><td>5.0</td><td>Fiducial distance to Rabbit PWN</td><td>Ng et al. (2005)</td></tr><tr><td>J1849-000</td><td>7.0</td><td>Scutum arm tangent region</td><td>Gotthelf et al. (2011)</td></tr></tbody></table>

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. Data extracted from the HGPS

The quantities taken from the HGPS catalogue are source position, extension, integral flux >1 TeV, and spectral index $\Gamma$ from the power-law fit of the differential photon flux $\phi_0 \times (E/E_0)^{-\Gamma}$. The extension measure $\sigma$ 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\sigma$ limits following a similar prescription, namely in the cases where the offset was less than $3\sigma$ 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\text{TeV}}$ and index $\Gamma$ is converted to a luminosity between 1 and 10 TeV using

$$L_{1-10\text{TeV}} = 1.92 \times 10^{44} \frac{I_{>1\text{TeV}}}{\text{cm}^{-2}\text{s}^{-1}} \times \frac{\Gamma - 1}{\Gamma - 2} (1 - 10^{2-\Gamma}) \left(\frac{d}{\text{kpc}}\right)^2 \text{erg s}^{-1}, \quad (1)$$

where $d$ is the source distance and the integral flux $I_{>1\text{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 $\sigma_{\Gamma}$ and integral flux errors $\sigma_I$, are

$$\left(\frac{\sigma_L}{L_{1-10\text{TeV}}}\right)^2 = \left[\frac{\sigma_I}{I_{>1\text{TeV}}}\right]^2 + \left[\left(\frac{1}{(\Gamma - 1)(\Gamma - 2)} + \frac{\ln 10}{1 - 10^{\Gamma - 2}}\right)\sigma_{\Gamma}\right]^2. \quad (2)$$
<table><caption>**Table 3.** Pulsar wind nebulae outside the HGPS catalogue.</caption><thead><tr><th>Canonical name</th><th>ATNF name</th><th>lg <i>E</i></th><th>τ<sub>c</sub><br>(kyr)</th><th>d<br>(kpc)</th><th>PSR offset<br>(pc)</th><th>Γ</th><th><i>R</i><sub>PWN</sub><br>(pc)</th><th><i>L</i><sub>1–10 TeV</sub><br>(10<sup>33</sup> erg s<sup>−1</sup>)</th></tr></thead><tbody><tr><td>N157B<sup>1</sup></td><td>J0537−6910</td><td>38.69</td><td>4.93</td><td>53.7</td><td>&lt;22</td><td>2.80 ± 0.10</td><td>&lt;94</td><td>760 ± 80</td></tr><tr><td>Crab Nebula<sup>2</sup></td><td>B0531+21</td><td>38.65</td><td>1.26</td><td>2.00</td><td>&lt;0.8</td><td>2.63 ± 0.02</td><td>&lt;3</td><td>32.1 ± 0.7</td></tr><tr><td>G0.9+0.1<sup>3</sup></td><td>J1747−2809</td><td>37.63</td><td>5.31</td><td>13.3</td><td>&lt;3</td><td>2.40 ± 0.11</td><td>&lt;7</td><td>46 ± 7</td></tr><tr><td>3C 58<sup>4</sup></td><td>J0205+6449</td><td>37.43</td><td>5.37</td><td>2.00</td><td>&lt;2</td><td>2.4 ± 0.2</td><td>&lt;5</td><td>0.23 ± 0.06</td></tr><tr><td>CTA 1<sup>5</sup></td><td>J0007+7303</td><td>35.65</td><td>13.9</td><td>1.40</td><td>&lt;4</td><td>2.2 ± 0.2</td><td>6.6 ± 0.5</td><td>0.71 ± 0.10</td></tr></tbody></table>

**Notes.** See Table 1 for the explanation of the columns. G0.9+0.1 is listed in the catalogue, but not treated in the HGPS analysis pipeline, so we treat it as an HGPS-external result. Offset limits were calculated as for the HGPS (see Sect. 2.3). In the case of N157B and 3C 58, $2\sigma_{PSF}$ was used as conservative extension limit since no value is given in the respective papers.

**References.** (1) H.E.S.S. Collaboration (2012c); Abramowski et al. (2015); (2) extension limit: Aharonian et al. (2004), flux: Aharonian et al. (2006b); (3) Aharonian et al. (2005a); (4) Aleksić et al. (2014); (5) Aliu et al. (2013).

The errors on flux and index are assumed to be independent because the reference energy of 1 TeV is typically very close to the 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<fnref n="4" />.

We also extract flux upper limits from the sky maps of the HGPS data release. The 95% confidence level limit $I_{>1\text{TeV}}$ on the flux is converted as above, assuming a spectral index of $\Gamma = 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⁻¹ for the offset speed (see Sect. 5.2.2), a required flux limit radius $R_{\text{lim}} = R_{\text{PWN}} + R_{\text{off}}$ is derived and a corresponding angular extent $\theta_{\text{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° < \theta_{\text{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°)^2$ to account for the uncontained part of the PWN. If $\theta_{\text{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 un bias 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 $\dot{E}$, which can be observationally determined from its period $P$ and period derivative $\dot{P}$, assuming a neutron star moment of inertia of $I = 10^{45}$ 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 $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 ($O(<10 \text{ 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),

<footnote n="4">⁴ Vela X is the only source where this prescription leads to a significant deviation from previously published dedicated analyses, both because of its energy cut-off and its extended emission component up to 1.2° away from its centre (Abramowski et al. 2012). Therefore, we convert its $I_{>1\text{TeV}}$ to an energy flux using its cut-off spectrum ($\Gamma = 1.35 \pm 0.08$; $\lambda = 0.0815 \pm 0.0115$ for a flux function $F(E) \sim E^{-\lambda}$ exp($-\lambda E$)), which leads to a 17% higher energy flux than when only using the power-law approximation. Furthermore, the extended and faint “ring” emission component noted in Abramowski et al. (2012) is taken into account by applying a correction factor of $1.31 \pm 0.16$. This emission component is derived from the ratio of “Inner” and “Total” integral fluxes presented in Abramowski et al. (2012), Table 3.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Histograms of spatial separation between PMPS pulsars and TeV source components from the HGPSC list. In the high-$\dot{E}$ pulsar sample (left), a clear correlation is seen as a peak at small squared angular distances, whereas the low-$\dot{E}$ associations, if present, are not significant beyond the expected rate of chance coincidences (right). The angular separation cut of $\theta < 0.5^\circ$ applied in the preselection of PWN candidates (Sect. 3.2) is indicated by a dashed vertical line in the left panel.</figcaption></figure>

where the whole HGPS sky map of that time was used along with the PMPS pulsar catalogue to evaluate a detection fraction $N_{\text{detected}}/N_{\text{pulsars}}$ for pulsars in different bands in $\dot{E}/d^2$.

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 $\dot{E}/d^2$ and all “Gaussian components” listed in the unbiased HGPSC component list<fnref n="5" />. The shaded band shows the expectation derived from simulated pulsar samples. It is derived for the same band of $\dot{E}/d^2$, 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 $\dot{E}/d^2 > 10^{34}$ erg s⁻¹ kpc⁻² and is absent for less energetic pulsars. An estimate for the number of chance coincidences for a cut of $0.5^\circ$ 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.

## 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 $\dot{E}/d^2 = 10^{34}$ erg s⁻¹ kpc⁻² and have an angular separation $\theta$ from an HGPS source of less than

$0.5^\circ$. We also require a characteristic age $\tau_c < 10^7$ yr to prevent millisecond pulsars, which are different concerning their nature and physics of emission, from entering the PWN candidate sample<fnref n="6" />. 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 $\dot{E} > 10^{35}$ erg s⁻¹, independent of 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 $\dot{E}/d^2$, 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 pulsars

<footnote n="5">⁵ We use $\dot{E}/d^2$ as an estimator for detectability for consistency with previous works. This is optimal under the assumptions that (a) the TeV luminosity scales linearly with $\dot{E}$, and (b) the sources appear small compared to the correlation radius. Both of these assumptions are questionable, given the large extension of some objects and the weak correlation between $\dot{E}$ and TeV luminosity. For this reason, we cross-checked the study with just $\dot{E}$ as the estimator, and we find very similar results. Presumably, the fact that $d$ only varies by a factor of 10 throughout the population makes the distance correction a subdominant effect against intrinsic luminosity variations.</footnote>

<footnote n="6">⁶ There is only one case of such a coincidence, PSR J1832–0836, which correlates with HESS J1832–085 along with the much more likely ordinary PSR B1830–08.</footnote>
<table><caption>**Table 4.** Candidate pulsar wind nebulae from the pre-selection.</caption><thead><tr><th>HGPS name</th><th>ATNF name</th><th>lg E&#x0307;</th><th>&tau;<sub>c</sub><br>(kyr)</th><th>d<br>(kpc)</th><th>PSR offset<br>(pc)</th><th>&Gamma;</th><th>R<sub>PWN</sub><br>(pc)</th><th>L<sub>1-10 TeV</sub><br>(10<sup>33</sup> erg s<sup>-1</sup>)</th><th>Rating<br>1 2 3 4</th></tr></thead><tbody><tr><td>J1616-508 (1)</td><td>J1617-5055</td><td>37.20</td><td>8.13</td><td>6.82</td><td>&lt;26</td><td>2.34 &pm; 0.06</td><td>28 &pm; 4</td><td>162 &pm; 9</td><td>★★★★</td></tr><tr><td>J1023-575</td><td>J1023-5746</td><td>37.04</td><td>4.60</td><td>8.00</td><td>&lt;9</td><td>2.36 &pm; 0.05</td><td>23.2 &pm; 1.2</td><td>67 &pm; 5</td><td>★★★★</td></tr><tr><td>J1809-193 (1)</td><td>J1811-1925</td><td>36.81</td><td>23.3</td><td>5.00</td><td>29 &pm; 7</td><td>2.38 &pm; 0.07</td><td>35 &pm; 4</td><td>53 &pm; 3</td><td>★★★☆</td></tr><tr><td>J1857+026</td><td>J1856+0245</td><td>36.66</td><td>20.6</td><td>9.01</td><td>21 &pm; 6</td><td>2.57 &pm; 0.06</td><td>41 &pm; 9</td><td>118 &pm; 13</td><td>★★★★</td></tr><tr><td>J1640-465</td><td>J1640-4631 (1)</td><td>36.64</td><td>3.35</td><td>12.8</td><td>&lt;20</td><td>2.55 &pm; 0.04</td><td>25 &pm; 8</td><td>210 &pm; 12</td><td>★★★★</td></tr><tr><td>J1641-462</td><td>J1640-4631 (2)</td><td>36.64</td><td>3.35</td><td>12.8</td><td>50 &pm; 5</td><td>2.50 &pm; 0.11</td><td>&lt;14</td><td>17 &pm; 4</td><td>☆ ★★★</td></tr><tr><td>J1708-443</td><td>B1706-44</td><td>36.53</td><td>17.5</td><td>2.60</td><td>17 &pm; 3</td><td>2.17 &pm; 0.08</td><td>12.7 &pm; 1.4</td><td>6.6 &pm; 0.9</td><td>★★★★</td></tr><tr><td>J1908+063</td><td>J1907+0602</td><td>36.45</td><td>19.5</td><td>3.21</td><td>21 &pm; 3</td><td>2.26 &pm; 0.06</td><td>27.2 &pm; 1.5</td><td>28 &pm; 2</td><td>★★★★</td></tr><tr><td>J1018-589A</td><td>J1016-5857 (1)</td><td>36.41</td><td>21.0</td><td>8.00</td><td>47.5 &pm; 1.6</td><td>2.24 &pm; 0.13</td><td>&lt;4</td><td>8.1 &pm; 1.4</td><td>☆ ★★★</td></tr><tr><td>J1018-589B</td><td>J1016-5857 (2)</td><td>36.41</td><td>21.0</td><td>8.00</td><td>25 &pm; 7</td><td>2.20 &pm; 0.09</td><td>21 &pm; 4</td><td>23 &pm; 5</td><td>★★★★</td></tr><tr><td>J1804-216</td><td>B1800-21</td><td>36.34</td><td>15.8</td><td>4.40</td><td>18 &pm; 5</td><td>2.69 &pm; 0.04</td><td>19 &pm; 3</td><td>42.5 &pm; 2.0</td><td>★★★★</td></tr><tr><td>J1809-193 (2)</td><td>J1809-1917</td><td>36.26</td><td>51.3</td><td>3.55</td><td>&lt;17</td><td>2.38 &pm; 0.07</td><td>25 &pm; 3</td><td>26.9 &pm; 1.5</td><td>★★★★</td></tr><tr><td>J1616-508 (2)</td><td>B1610-50</td><td>36.20</td><td>7.42</td><td>7.94</td><td>60 &pm; 7</td><td>2.34 &pm; 0.06</td><td>32 &pm; 5</td><td>220 &pm; 12</td><td>☆ ★★★</td></tr><tr><td>J1718-385</td><td>J1718-3825</td><td>36.11</td><td>89.5</td><td>3.60</td><td>5.4 &pm; 1.6</td><td>1.77 &pm; 0.06</td><td>7.2 &pm; 0.9</td><td>4.6 &pm; 0.8</td><td>★★★★</td></tr><tr><td>J1026-582</td><td>J1028-5819</td><td>35.92</td><td>90.0</td><td>2.33</td><td>9 &pm; 2</td><td>1.81 &pm; 0.10</td><td>5.3 &pm; 1.6</td><td>1.7 &pm; 0.5</td><td>☆ ★★★</td></tr><tr><td>J1832-085</td><td>B1830-08 (1)</td><td>35.76</td><td>147</td><td>4.50</td><td>23.3 &pm; 1.5</td><td>2.38 &pm; 0.14</td><td>&lt;4</td><td>1.7 &pm; 0.4</td><td>☆ ☆ ★★</td></tr><tr><td>J1834-087</td><td>B1830-08 (2)</td><td>35.76</td><td>147</td><td>4.50</td><td>32.3 &pm; 1.9</td><td>2.61 &pm; 0.07</td><td>17 &pm; 3</td><td>25.8 &pm; 2.0</td><td>☆ ★★☆</td></tr><tr><td>J1858+020</td><td>J1857+0143</td><td>35.65</td><td>71.0</td><td>5.75</td><td>38 &pm; 3</td><td>2.39 &pm; 0.12</td><td>7.9 &pm; 1.6</td><td>7.1 &pm; 1.5</td><td>☆ ★★☆</td></tr><tr><td>J1745-303</td><td>B1742-30 (1)</td><td>33.93</td><td>546</td><td>0.200</td><td>1.42 &pm; 0.15</td><td>2.57 &pm; 0.06</td><td>0.62 &pm; 0.07</td><td>0.014 &pm; 0.003</td><td>☆ ☆ ☆</td></tr><tr><td>J1746-308</td><td>B1742-30 (2)</td><td>33.93</td><td>546</td><td>0.200</td><td>&lt;1.1</td><td>3.3 &pm; 0.2</td><td>0.56 &pm; 0.12</td><td>0.009 &pm; 0.003</td><td>★ ★☆ ☆</td></tr></tbody></table>

**Notes.** See Table 1 for the explanation of the columns. In the rating columns (1: PSR containment, 2: extension, 3: luminosity, 4: surface brightness, see Sect. 6), a big star ★ denotes a quantity that fulfills its requirement, a small star ★ denotes a compatible limit, a lightning symbol ☆ denotes a limit or measurement in conflict with the requirement (see Sect. 6). Numbers in brackets indicate double associations.

(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 the average 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-E pulsars known today.

## 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 E determines the energy injection history of a PWN. It is decaying continually at a rate determined by the so-called spin-down timescale $\tau_c$, following

an evolution similar to that expected from a dipole (see also Appendix B)

$$ \dot{E}(t) = \dot{E}_0 \left(1 + \frac{t}{\tau_0}\right)^{\frac{n+1}{n-1}}, \qquad (3) $$

where $\tau_0$ is the initial spin-down timescale, $\dot{E}_0$ 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 $\tau_0 \sim 10^{2.5-3.5}$ yr, $\dot{E}_0 \sim 10^{37.5-40}$ erg s$^{-1}$, and $n \sim 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 ($E_{\text{rot}} = \dot{E}_0\tau_0(n-1)/2 = I\Omega_0^2/2$, typically $\le 10^{50}$ 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 $\dot{P}$ (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

$$ \tau_c \equiv \frac{P}{2\dot{P}} = (\tau_0 + t) \frac{n-1}{2}. \qquad (4) $$

If $t \gg \tau_0$ and $n=3$, then $\tau_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 $\dot{E}$ and $\tau_c$, i.e.

$$ \dot{E} = \dot{E}_0 \left[ \frac{2}{n-1} \cdot \frac{\tau_c}{\tau_0} \right]^{\frac{n+1}{n-1}}, \qquad (5) $$

or, equivalently, between $\dot{P}$ and $P$ (see Eq. (B.12) in Appendix B), i.e.

$$ \dot{P}(P) = \frac{P_0}{\tau_0} \frac{1}{n-1} \left( \frac{P}{P_0} \right)^{2-n}. \qquad (6) $$
<table><caption>**Table 5.** Flux and luminosity upper limits (95% CL) for regions around pulsars without detected PWN.</caption><thead><tr><th>ATNF name</th><th>lg <i>E</i></th><th>τ<sub>c</sub><br>(kyr)</th><th>d<br>(kpc)</th><th>θ<sub>pred</sub><br>(deg)</th><th>θ<sub>int</sub><br>(deg)</th><th>Significance<br>(σ)</th><th><i>F</i><sub>&gt;1 TeV</sub><br>(10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th><i>L</i><sub>1-10 TeV</sub><br>(10<sup>33</sup> erg s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J1400-6325</td><td>37.71</td><td>12.7</td><td>7.00</td><td>0.150</td><td>0.2</td><td>1.4</td><td>&lt;0.41</td><td>&lt;8.3</td></tr><tr><td>J1124-5916</td><td>37.08</td><td>2.85</td><td>5.00</td><td>0.137</td><td>0.2</td><td>1.0</td><td>&lt;0.27</td><td>&lt;2.8</td></tr><tr><td>J1410-6132</td><td>37.00</td><td>24.8</td><td>15.6</td><td>0.127</td><td>0.2</td><td>2.8</td><td>&lt;0.53</td><td>&lt;54</td></tr><tr><td>J1935+2025</td><td>36.67</td><td>20.9</td><td>6.21</td><td>0.29</td><td>0.4</td><td>1.9</td><td>&lt;0.88</td><td>&lt;14</td></tr><tr><td>J1112-6103</td><td>36.65</td><td>32.7</td><td>12.2</td><td>0.21</td><td>0.4</td><td>3.7</td><td>&lt;1.0</td><td>&lt;62</td></tr><tr><td>J1801-2451</td><td>36.41</td><td>15.5</td><td>5.22</td><td>0.30</td><td>0.4</td><td>1.1</td><td>&lt;0.56</td><td>&lt;6.3</td></tr><tr><td>J1837-0604</td><td>36.30</td><td>33.8</td><td>6.41</td><td>0.42</td><td>0.4</td><td>9.5</td><td>&lt;2.1</td><td>&lt;36</td></tr><tr><td>J1341-6220</td><td>36.15</td><td>12.1</td><td>11.1</td><td>0.129</td><td>0.2</td><td>2.6</td><td>&lt;0.46</td><td>&lt;24</td></tr><tr><td>J1055-6028</td><td>36.08</td><td>53.5</td><td>15.5</td><td>0.25</td><td>0.4</td><td>1.1</td><td>&lt;0.70</td><td>&lt;70</td></tr><tr><td>J1934+2352</td><td>35.96</td><td>21.6</td><td>11.6</td><td>0.175</td><td>0.2</td><td>1.6</td><td>&lt;1.1</td><td>&lt;64</td></tr><tr><td>J1932+2220</td><td>35.88</td><td>39.8</td><td>10.9</td><td>0.29</td><td>0.4</td><td>-0.9</td><td>&lt;0.55</td><td>&lt;27</td></tr><tr><td>J1702-4310</td><td>35.80</td><td>17.0</td><td>5.14</td><td>0.35</td><td>0.4</td><td>0.9</td><td>&lt;0.59</td><td>&lt;6.5</td></tr><tr><td>J1413-6141</td><td>35.75</td><td>13.6</td><td>10.1</td><td>0.161</td><td>0.2</td><td>2.8</td><td>&lt;0.54</td><td>&lt;23</td></tr><tr><td>J1909+0749</td><td>35.65</td><td>24.7</td><td>9.48</td><td>0.24</td><td>0.4</td><td>0.9</td><td>&lt;0.41</td><td>&lt;15</td></tr><tr><td>J1815-1738</td><td>35.59</td><td>40.4</td><td>8.78</td><td>0.37</td><td>0.4</td><td>8.9</td><td>&lt;2.1</td><td>&lt;68</td></tr><tr><td>J1646-4346</td><td>35.56</td><td>32.5</td><td>5.79</td><td>0.48</td><td>0.4</td><td>-2.0</td><td>&lt;0.27</td><td>&lt;3.8</td></tr><tr><td>J1850-0026</td><td>35.52</td><td>67.5</td><td>11.1</td><td>0.44</td><td>0.4</td><td>3.7</td><td>&lt;0.91</td><td>&lt;46</td></tr><tr><td>J1907+0918</td><td>35.51</td><td>38.0</td><td>7.79</td><td>0.40</td><td>0.4</td><td>2.7</td><td>&lt;0.61</td><td>&lt;15</td></tr><tr><td>J1406-6121</td><td>35.34</td><td>61.7</td><td>8.15</td><td>0.56</td><td>0.4</td><td>4.4</td><td>&lt;3.3</td><td>&lt;91</td></tr><tr><td>J1412-6145</td><td>35.08</td><td>50.6</td><td>7.82</td><td>0.51</td><td>0.4</td><td>5.0</td><td>&lt;3.0</td><td>&lt;75</td></tr><tr><td>J1550-5418</td><td>35.00</td><td>1.41</td><td>4.00</td><td>0.29</td><td>0.4</td><td>1.0</td><td>&lt;0.47</td><td>&lt;3.1</td></tr><tr><td>J1841-0524</td><td>35.00</td><td>30.2</td><td>5.34</td><td>0.53</td><td>0.4</td><td>20.9</td><td>&lt;7.3</td><td>&lt;86</td></tr></tbody></table>

**Notes.** In addition to the table variables explained in Table 1, $\theta_{pred}$ is the predicted PWN extension (including offset), $\theta_{int}$ is the correlation radius of the map where the limit is taken from, and $F_{>1\text{TeV}}$ is the actual flux limit (see Sect. 2.3 for details). In the cases of high significance, the pulsar coincides with a TeV source that is not considered the PWN.

Consequently, the power indices of the above relations are only determined by the braking index $n$. Figures 2 show how real 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 $\dot{E}$ and $\tau_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 \sim 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 \sim t^{0.73}$ for $t < \tau_0$ in van der Swaluw et al. (2001) and $R \sim t^{0.3}$ for $t > \tau_0$ 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 collision of 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 plasma bubble while the old nebula from its younger period still remains, typically as an IC-dominated PWN due to its much lower magnetisation.
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Schematic of the Milky Way and its spiral arms, along with firmly identified PWNe, candidates, and energetic pulsars ($\dot{E} > 10^{35}$ 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).</figcaption></figure>

### 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 simplified, 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$^7$, 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 $\tau_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 expecta-

tion (*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 $\tau_c$ and $\dot{E}_0$, 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 $\dot{E}_0$ and $\tau_0$, strongly restricted by Fig. 2 (left), may be larger if the two quantities were anti-correlated such that high-$\dot{E}_0$ pulsars always tend to have a lower $\tau_0$; this is physically plausible because the two quantities are related through the pulsar birth period and magnetic field. On the mathematical side, $\dot{E}_0$, $\eta$, the energy injection range and the background photon density are all parameters with which the TeV luminosity scales in an almost

<footnote>$^7$ A spectral break at lower injection energies is generally necessary to model low-energy data, but since this does not impact the TeV regime, and we therefore cannot constrain it with the data presented in this paper, we focus on the VHE part with a single power law.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 4. Top: PWN extension occurrences over distance from Earth, in comparison to the band of extensions that can be expected to be identified in the HGPS analysis chain. Bottom: distribution of known distances of energetic pulsars ($\dot{E} > 10^{35}$ erg s$^{-1}$).</figcaption></figure>

linear way. In our varied model, we deal with this redundancy by only varying $\dot{E}_0$, 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 $\tau_c$ or $\dot{E}$ 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

$$ \lg Y_{est} = p_0 + p_1 \lg X. \qquad (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 $\sigma_{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 $\chi^2$ 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).

### 5.2.1. Extension

Figure 5 (left) shows the evolution of PWN extension as a function of characteristic age $\tau_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 $\dot{E}$.

A caveat is that there is a selection bias from the fact that extension estimates or limits are only available for sources that are
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Left: PWN extension evolution with time, in comparison to the modelling considered in this work. Right: PWN extension evolution with $\dot{E}$, as fitted in the $R_{PWN}(\dot{E})$ column of Table 6 for pulsar wind nebulae with $\dot{E} > 10^{36}$ erg s$^{-1}$ (see Sect. 5.2.1). The shaded range shows the fit range and standard deviation $\sigma_{lg R}$. 1 dex refers to an order of magnitude and is the unit of the logspace defined $\sigma_{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.</figcaption></figure>

<table><caption>**Table 6.** ASURV fit results.</caption><thead><tr><th></th><th><i>R</i><sub>PWN</sub>(<i>È</i>)</th><th><i>R</i><sub>PWN</sub>(<i>τ</i><sub>c</sub>)</th><th><i>L</i><sub>1–10 TeV</sub>(<i>È</i>)</th><th><i>L</i><sub>1–10 TeV</sub>(<i>τ</i><sub>c</sub>)</th><th><i>S</i>(<i>È</i>)</th><th><i>d</i><sub>P–P</sub>(<i>ε</i><sub>TeV</sub>)</th><th><i>d</i><sub>P–P</sub>(<i>È</i>)</th><th><i>d</i><sub>P–P</sub>(<i>τ</i><sub>c</sub>)</th></tr></thead><tbody><tr><td><i>p</i>-value</td><td>0.012</td><td>0.047</td><td>0.010</td><td>0.13</td><td>0.0013</td><td>0.0004</td><td>0.035</td><td>0.0086</td></tr><tr><td>σ<sub>lg Y</sub></td><td>0.32</td><td>0.39</td><td>0.83</td><td>0.91</td><td>0.28</td><td>0.18</td><td>0.49</td><td>0.42</td></tr><tr><td><i>p</i><sub>0</sub></td><td>1.48 ± 0.20</td><td>0.38 ± 0.22</td><td>33.22 ± 0.27</td><td>34.1 ± 0.4</td><td>30.62 ± 0.13</td><td>1.97 ± 0.16</td><td>1.07 ± 0.25</td><td>-0.9 ± 0.5</td></tr><tr><td><i>p</i><sub>1</sub></td><td>-0.65 ± 0.20</td><td>0.55 ± 0.23</td><td>0.59 ± 0.21</td><td>-0.46 ± 0.36</td><td>0.81 ± 0.14</td><td>0.52 ± 0.07</td><td>-0.75 ± 0.29</td><td>1.4 ± 0.5</td></tr></tbody></table>

**Notes.** $p_0$ and $p_1$ relate to Eq. (7). The *p*-value is calculated after the Cox proportional hazards model. The fit used (within ASURV) is the “EM algorithm”. *P* is given in 0.1 s, $\dot{P}$ in $10^{-13}$ s s$^{-1}$, $\dot{E}$ in $10^{36}$ erg s$^{-1}$, and $\tau_c$ in kyr. $R_{PWN}$ is given in pc, $L_{1-10\,TeV}$ in erg s$^{-1}$, $S$ in erg s$^{-1}$ pc$^{-2}$. The 2D pulsar-PWN offset $d_{P-P}$ is given in parsecs, and $\epsilon_{TeV} = L_{1-10\,TeV}/\dot{E}$ is the apparent TeV efficiency.

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 un bias the sample in the fitting procedure below, we apply a cut of $\dot{E} > 10^{36}$ erg s$^{-1}$, beyond which the likelihood 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 $\tau_c$ or $\dot{E}$ 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 $\dot{E}$ and increasing $\tau_c$. A more general 2D fit of $R_{PWN}(P, \dot{P})$ does not lead to a significant improvement of the fit, nor a lower *p*-value. The parametrisation of $R_{PWN}(\dot{E})$ is shown in Fig. 5 (right) to show that it is indeed suitable for predicting the extensions of the detected young PWNe ($\dot{E} > 10^{36}$ erg s$^{-1}$) reasonably well. The only PWN below $10^{36}$ 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 \sim \tau_c^{0.55 \pm 0.23}$ can be compared to the baseline model in Fig. 5 (left), which assumes the canonical $R \sim 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 $\tau_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 con-
sistent 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 gener-
ally grow with time until an age of few tens of kiloyears is clear
and supported by the ASURV fits. There are, however, few pul-
sar 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, asymmet-
ric crushing of the PWN by the surrounding SNR and, possi-
bly, by asymmetric pulsar outflow. Hobbs et al. (2005) deter-
mined the mean 2D speed for non-millissecond pulsars to be
307 ± 47 km s⁻¹, 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 distribu-
tion and high-velocity outliers, but it is clear that the bulk of
the pulsars have 2D velocities of less than 500 km s⁻¹. 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 increas-
ing offset with rising τc and falling E 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 expan-
sion 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 av-
erage, it is unlikely that such a high fraction of the high-E pul-
sars with TeV-detected nebulae are so fast (1000 km s⁻¹ 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-aged wind 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 contain-
ment ratio as the PWN offset divided by the PWN extension ra-

dius. Given the offset and extension evolution discussed above,
a pulsar is not expected to leave its (then relic) wind nebula be-
fore 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 exten-
sion 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 as-
sociation of a pulsar with a TeV object; for old systems, how-
ever, these MWL data are very difficult to acquire since the syn-
chrotron 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 character-
istic 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^{-12}$ erg s$^{-1}$ cm$^{-2}$ (H.E.S.S. Collaboration 2018), which is equivalent to a luminosity threshold of $3 \times 10^{33}$ 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 $\dot{E} > 10^{35}$ erg s$^{-1}$. As explained in Sects. 2.3 and 3.2, about one-third of these high-$\dot{E}$ 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 $\sim 10^{36}$ 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 cor-
relation of luminosity with pulsar spin-down<fnref n="8" />. The ASURV fit
suggests a relation of $L_{1-10\text{TeV}} \sim \dot{E}^{0.59 \pm 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

<footnote n="8">⁸ The *p*-value without N157B is still 0.06, so the correlation does not
only depend on this one source.</footnote>
<figure><img src="image_6.png" /><figcaption>**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⁻¹). *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).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** *Left:* relation of TeV luminosity and pulsar *E*. *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).</figcaption></figure>

spin-down power falls below ~$10^{36}$ erg s⁻¹. 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<fnref n="9" />. Figure 7 (right) shows the result of the fitted parametrisation $L_{1-10\,\text{TeV}}(\dot{E})$ derived from our data.

In contrast, the $L_{1-10\,\text{TeV}}$ over $\tau_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, *E* is the better variable to characterise the evolutionary state of the PWN luminosity.

### 5.3.2. Apparent TeV efficiency

The TeV efficiency, conventionally defined as $\epsilon_{\text{TeV}} = L_{1-10\,\text{TeV}}/\dot{E}$, is not the real present efficiency of a PWN because $L_{1-10\,\text{TeV}}$ is a result of the whole injection history, whereas $\dot{E}$ 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

<footnote n="9">⁹ The *p*-value for the fit of $L_{1-10\,\text{TeV}}(\dot{E})$ without the limits is 0.31.</footnote>
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

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.

### 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{L_{1-10 \text{ TeV}}}{4\pi R_{\text{PWN}}^2} \approx \frac{F_{1-10 \text{ TeV}}}{\sigma^2}, \quad (8)$$

where $R_{\text{PWN}}$ is the physical PWN radius (in pc), $\sigma$ is its angular extent as seen from Earth, and $F_{1-10\text{TeV}}$ is the integral energy flux between 1 TeV and 10 TeV measured at Earth.

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

As seen in the ASURV fit values in Table 6, a comparatively strong correlation is found, confirming the above findings of a decreasing luminosity and increasing extension of aging pulsar wind nebulae. The measured power-law relation of $S \sim \dot{E}^{0.81 \pm 0.14}$ matches what the model suggests (~0.9 for the part where $\dot{E} < 2 \times 10^{37}$ 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 ($\sigma_{1gS} \sim 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.

### 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 $\dot{E}$.
<figure><img src="image_9.png" /><figcaption>**Fig. 11.** Common logarithmic residuals of rating criteria 2–4, using the standard deviations $\sigma_{lgY}$ 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.</figcaption></figure>

**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 observations at 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 standard deviations, using the measured $\sigma_{lgR} = 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 standard deviations, using the measured $\sigma_{lgL} = 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 $\sigma_{lgS} = 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 as-
sociated with the very nearby old PSR B1742−30, are a spe-
cial case. The pulsar is a factor of 10–100 older than most other
PWNe discussed here, so the extrapolation performed for the rat-
ing 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).

7. Conclusions

In this work we subsume and examine the population of TeV
pulsar wind nebulae found to date. The census presents 14 ob-
jects reanalysed in the HGPS catalogue pipeline, which we con-
sider 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⁻¹ 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 $\dot{E} \ge 10^{37}$ 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\text{TeV}} < 5.7 \times 10^{34}$ erg s$^{-1}$, which is compatible with the predicted $3.3 \times 10^{34}$ erg s$^{-1}$ that can be taken from $L_{1-10\text{TeV}}(\dot{E})$ in Table 6.

In summary, only 5 of the 17 highest-$\dot{E}$ 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 parame-
ters. The main conclusion of this work is that for several observ-
ables, 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 popula-
tion of detected and undetected PWNe suggests.

Using the flux limits for undetected PWNe, we find evi-
dence 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⁻¹ (roughly corre-
sponding to several tens of kiloyears). This was implicitly known
before from the mere non-detection of old pulsar wind nebulae,
but for the first time could be put into a quantitative perspec-
tive 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\text{TeV}} \sim \dot{E}^{0.58 \pm 0.21}$, in consistency with the model that sug-
gests a power index of around 0.5.

Another feature that was discussed on some individual ob-
jects 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 de-
veloping asymmetrically, for instance due to an inhomogeneous
surrounding medium (ISM), this crushing may result in consid-
erable distortion and displacement of the wind nebula. Put to a
population-scoped graph (Fig. 6, left), it becomes clear that pul-
sar proper motions are insufficient to explain the large offsets
observed, which may instead be due to reverse shock interac-
tion being a dominant and frequent cause of pulsar-PWN off-
set in middle-aged systems (see also 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 be-
tween 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 characteris-
tic age and falling pulsar spin-down) could also be shown to
be evident in the data. The fitted relation $R \sim \tau_c^{0.55 \pm 0.23}$ sug-
gests 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 con-
sistent 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 re-
ported 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 sur-
face brightness and pulsar $\dot{E}$ (Fig. 10, left). What stands out is
not only the correlation itself, but also its relatively low scat-
ter. 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 particu-
larly clear about it for the young to middle-aged PWNe we in-
vestigated. A more sensitive data set – as expected from CTA
(Acharya et al. 2013) – will reduce the uncertainties on spec-
tral 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^{36}$ erg s⁻¹) 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 posi-
tion 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 syn-
thesis study, one could use a specific Milky Way model to "cali-
brate" 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 popula-
tion 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<fnref n="10" />. 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<fnref n="11" />, 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 multididi-
mensional 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 mod-
elling code is needed that includes the difficult reverse shock in-
teraction 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 to reliably 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 con-
tribute 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 capa-
bility 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.

Acknowledgements. 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 benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation.

References

Abdo, A. A., Allen, B. T., Aune, T., et al. 2009, ApJ, 700, L127

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17

Abysekara, A. U., Alfaro, R., Alvarez, C., et al. 2015, ArXiv e-print
[arXiv:1508.03327]

Abramowski, A., Acero, F., Aharonian, F., et al. 2012, A&A, 548, A38

Abramowski, A., Aharonian, F., Benkhali, F. A., et al. 2014, MNRAS, 439, 2828

Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2015, Science, 347, 406

Acciari, V. A., Aliu, E., Arlen, T., et al. 2010, ApJ, 719, L69

Acero, F., Ackermann, M., Ajello, M., et al. 2013, ApJ, 773, 77

Acharya, B. S., Actis, M., Aghajani, T., et al. 2013, Astropart. Phys., 43, 3

Aharonian, F. A., & Atoyan, A. M. 1995, International Cosmic Ray Conference, 2, 358

Aharonian, F. A., Atoyan, A. M., & Kifune, T. 1997, MNRAS, 291, 162

Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2004, ApJ, 614, 897

Aharonian, F., Akhperjanian, A. G., Aye, K.-M., et al. 2005a, A&A, 432, L25

Aharonian, F., Akhperjanian, A. G., Aye, K.-M., et al. 2005b, A&A, 435, L17

Aharonian, F. A., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2005c, A&A, 442, L25

<footnote n="10">¹⁰ All plotted parameters were derived from the four parameters *P*, *Ṗ*, *L*₁⁻₁₀TeV, and the PWN extension; the TeV offset was not dealt with in the modelling.</footnote>

<footnote n="11">¹¹ Even considering only the four papers mentioned in Sect. 4.1, Ė₀ and τ₀ vary there by factors of 250 and 6, respectively, compared to 10 and 1.4 in our work (see Table A.1).</footnote>
Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006a, A&A, **448**, L43

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006b, A&A, **457**, 899

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006c, ApJ, **636**, 777

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006d, A&A, **460**, 365

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006e, A&A, **456**, 245

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, A&A, **472**, 489

Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, **499**, 723

Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014, A&A, **567**, L8

Aliu, E., Archambault, S., Arlen, T., et al. 2013, ApJ, **764**, 38

Aliu, E., Archambault, S., Aune, T., et al. 2014, ApJ, **787**, 166

Arzoumanian, Z., Chernoff, D. F., & Cordes, J. M. 2002, ApJ, **568**, 289

Bamba, A., Anada, T., Dotani, T., et al. 2010, ApJ, **719**, L116

Baughman, B. M., Wood, J., & for the HAWC Collaboration. 2015, ArXiv e-prints [arXiv:1508.03497]

Blondin, J. M., Chevalier, R. A., & Frierson, D. M. 2001, ApJ, **563**, 806

Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., **42**, 237

Bucciantini, N. 2011, Astrophys. Space Sci. Proc., **21**, 473

Carrigan, S. 2007, Ph.D. Thesis, Ruperto-Carola University of Heidelberg

Carrigan, S., Hinton, J. A., Hofmann, W., et al. 2008, in International Cosmic Ray Conference, **2**, 659

Chevalier, R. A. 1977, in Astrophys. Space Sci. Lib. 66, Supernovae, ed. D. N. Schramm, **53**

Cordes, J. M., & Lazio, T. J. W. 2002, ArXiv e-prints [arXiv:astro-ph/0207156]

de Jager, O. C., & Djannati-Atai, A. 2009, in Astrophys. Space Sci. Lib., **357**, ed. W. Becker, **451**

de Jager, O. C., & Harding, A. K. 1992, ApJ, **396**, 161

De Jager, O. C., Harding, A. K., Baring, M. G., & Mastichiadis, A. 1995, International Cosmic Ray Conference, **2**, 528

Djannati-Atai, A., de Jager, O. C., Terrier, R., Gallant, Y. A., & Hoppe, S. 2008, in International Cosmic Ray Conference, **2**, 823

Dubois, P. F., Hinsen, K., & Hugunin, J. 1996, Comput. Phys., **10**

Du Plessis, I., de Jager, O. C., Buchner, S., et al. 1995, ApJ, **453**, 746

Faucher-Giguère, C.-A., & Kaspi, V. M. 2006, ApJ, **643**, 332

Funk, S., Hinton, J. A., Moriguchi, Y., et al. 2007, A&A, **470**, 249

Gaensler, B. M., & Slane, P. O. 2006, ARA&A, **44**, 17

Gelfand, J. D., Slane, P. O., & Zhang, W. 2009, ApJ, **703**, 2051

Gonthier, P. L., Van Guilder, R., & Harding, A. K. 2004, ApJ, **604**, 775

Gotthelf, E. V., & Halpern, J. P. 2008, ApJ, **681**, 515

Gotthelf, E. V., Halpern, J. P., Terrier, R., & Mattana, F. 2011, ApJ, **729**, L16

Gould, R. J. 1965, Phys. Rev. Lett., **15**, 577

Gunn, J. E., & Ostriker, J. P. 1969, Nature, **221**, 454

H.E.S.S. Collaboration (Abramowski, A., Acero, F., et al.) 2011a, A&A, **528**, A143

H.E.S.S. Collaboration (Abramowski, A., Acero, F., et al.) 2011b, A&A, **525**, A46

H.E.S.S. Collaboration (Abramowski, A., Acero, F., et al.) 2012a, A&A, **548**, A46

H.E.S.S. Collaboration (Abramowski, A., Acero, F., et al.) 2012b, A&A, **541**, A5

H.E.S.S. Collaboration (Abramowski, A., Acero, F., et al.) 2012c, A&A, **545**, L2

H.E.S.S. Collaboration (Abdalla, H., Abramowski, A., et al.) 2018, A&A, **612**, A1 (H.E.S.S. SI)

Hessels, J. W. T., Nice, D. J., Gaensler, B. M., et al. 2008, ApJ, **682**, L41

Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, **360**, 974

Horan, D., & Wakely, S. 2016, TeVCat online catalog for TeV Astronomy (version 3.400)

Hunter, J. D. 2007, Comp. Sci. Eng., **9**, 90

Isobe, T., Feigelson, E. D., & Nelson, P. I. 1986, ApJ, **306**, 490

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python

Kargaltsev, O., & Pavlov, G. G. 2010, X-ray Astronomy 2009; Present Status, Multi-Wavelength Approach and Future Perspectives, **1248**, 25

Kargaltsev, O., Rangelov, B., & Pavlov, G. G. 2013, ArXiv e-prints [arXiv:1305.2552]

Kennel, C. F., & Coroniti, F. V. 1984a, ApJ, **283**, 694

Kennel, C. F., & Coroniti, F. V. 1984b, ApJ, **283**, 710

Klepser, S. 2012, Astropart. Phys., **36**, 64

Knödléseder, J., Mayer, M., Deil, C., et al. 2013, in International Cosmic Ray Conference [arXiv:1307:5560]

Kothes, R. 2013, A&A, **560**, A18

Lavalley, M., Isobe, T., & Feigelson, E. 1992, in Astronomical Data Analysis Software and Systems I, eds. D. M. Worrall, C. Biemesderfer, & J. Barnes, ASP Conf. Ser., **25**, 245

Lorimer, D. R., Faulkner, A. J., Lyne, A. G., et al. 2006, MNRAS, **372**, 777

Magalhaes, N. S., Miranda, T. A., & Frajuca, C. 2012, ApJ, **755**, 54

Manchester, R. N., Lyne, A. G., Camilo, F., et al. 2001, MNRAS, **328**, 17

Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, **129**, 1993

Marandon, V. 2010, Ph.D. Thesis, Observatoire de Paris, Université Paris 7

Martín, J., Torres, D. F., & Rea, N. 2012, MNRAS, **427**, 415

Mattana, F., Falanga, M., Götz, D., et al. 2009, ApJ, **694**, 12

Mayer, M. 2010a, Diploma thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg

Mayer, M. 2010b, Ph.D. Thesis, U. Erlangen-Nuremberg (main)

Mayer, M., Brucker, J., Holler, M., et al. 2012, ArXiv e-print [arXiv:1202.1455]

Ng, C.-Y., Roberts, M. S. E., & Romani, R. W. 2005, ApJ, **627**, 904

Pacini, F., & Salvati, M. 1973, ApJ, **186**, 249

Porter, T. A., & Strong, A. W. 2005, in International Cosmic Ray Conf., **4**, 77

Possenti, A., Cerutti, R., Colpi, M., & Mereghetti, S. 2002, A&A, **387**, 993

Rauw, G., Manfroid, J., Gosset, E., et al. 2007, A&A, **463**, 981

Rees, M. J., & Gunn, J. E. 1974, MNRAS, **167**, 1

Renaud, M., Hoppe, S., Komin, N., et al. 2008, in AIP Conf. Ser. 1085, eds. F. A. Aharonian, W. Hofmann, & F. Rieger, 285

Reynolds, S. P., & Chevalier, R. A. 1984, ApJ, **278**, 630

Roberts, M. S. E. 2004, The Pulsar Wind Nebula Catalog (March 2005 version)

Spengler, G. 2015, Astropart. Phys., **67**, 70

Taylor, A. R. 2012, in Proc. IAU Symp., **291**

Temim, T., Slane, P., Kolb, C., et al. 2015, ApJ, **808**, 100

Terrier, R., Mattana, F., Djannati-Atai, A., et al. 2008, in AIP Conf. Ser. 1085, eds. F. A. Aharonian, W. Hofmann, & F. Rieger, 312

Torres, D. F., Cillis, A., Martín, J., & de Oña Wilhelmi, E. 2014, J. High Energy Astrophys., **1**, 31

Vallée, J. P. 2008, AJ, **135**, 1301

van der Swaluw, E., Achterberg, A., Gallant, Y. A., & Tóth, G. 2001, A&A, **380**, 309

van der Swaluw, E., Downes, T. P., & Keegan, R. 2004, A&A, **420**, 937

Vernetto, S., Chen, S., Di Sciascio, G., et al. 2013, in International Cosmic Ray Conference

Vink, J., Bamba, A., & Yamazaki, R. 2011, ApJ, **727**, 131

Vorster, M. J., Tibolla, O., Ferreira, S. E. S., & Kaufmann, S. 2013, ApJ, **773**, 139

Weekes, T. C., Cawley, M. F., Fegan, D. J., et al. 1989, ApJ, **342**, 379

Yin, P.-F., Yu, Z.-H., Yuan, Q., & Bi, X.-J. 2013, Phys. Rev. D, **88**, 023001

Yue, Y. L., Xu, R. X., & Zhu, W. W. 2007, Adv. Space Res., **40**, 1491

Zhang, L., Chen, S. B., & Fang, J. 2008, ApJ, **676**, 1210

<footnote n="1">¹ Centre for Space Research, North-West University, 2520 Potchefstroom, South Africa</footnote>

<footnote n="2">² Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany</footnote>

<footnote n="3">³ Max-Planck-Institut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany</footnote>

<footnote n="4">⁴ Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland</footnote>

<footnote n="5">⁵ National Academy of Sciences of the Republic of Armenia, Marshall Baghramian Avenue, 24, 0019 Yerevan, Republic of Armenia</footnote>

<footnote n="6">⁶ Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia</footnote>

<footnote n="7">⁷ Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany</footnote>

<footnote n="8">⁸ University of Namibia, Department of Physics, Private Bag 13301 Windhoek, Namibia</footnote>

<footnote n="9">⁹ GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands</footnote>

<footnote n="10">¹⁰ Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden</footnote>

<footnote n="11">¹¹ Institut für Theoretische Physik, Lehrstuhl IV: Weltraum und Astrophysik, Ruhr-Universität Bochum, 44780 Bochum, Germany</footnote>

<footnote n="12">¹² GRAPPA, Anton Pannekoek Institute for Astronomy and Institute of High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands</footnote>

<footnote n="13">¹³ Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität Innsbruck, 6020 Innsbruck, Austria</footnote>
14 School of Physical Sciences, University of Adelaide, 5005 Adelaide, Australia

15 LUTH, Observatoire de Paris, PSL Research University, CNRS, Université Paris Diderot, 5 place Jules Janssen, 92190 Meudon, France

16 Sorbonne Universités, UPMC Université Paris 06, Université Paris Diderot, Sorbonne Paris Cité, CNRS, Laboratoire de Physique Nucléaire et de Hautes Énergies (LPNHE), 4 place Jussieu, 75252 Paris Cedex 5, France

17 Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France

18 DSM/Irfu, CEA Saclay, 91191 Gif-Sur-Yvette Cedex, France

19 Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland

20 Aix-Marseille Université, CNRS/IN2P3, CPPM UMR 7346, 13288 Marseille, France

21 Instytut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

22 Funded by EU FP7 Marie Curie, grant agreement No. PIEF-GA-2012-332350

23 School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, 2050 Johannesburg, South Africa

24 Laboratoire d'Annecy-le-Vieux de Physique des Particules, Université Savoie Mont-Blanc, CNRS/IN2P3, 74941 Annecy-le-Vieux, France

25 Landessternwarte, Universität Heidelberg, Königstuhl, 69117 Heidelberg, Germany

26 Université Bordeaux, CNRS/IN2P3, Centre d'Études Nucléaires de Bordeaux Gradignan, 33175 Gradignan, France

27 Oskar Klein Centre, Department of Physics, Stockholm University, Albanova University Center, 10691 Stockholm, Sweden

28 Wallenberg Academy Fellow

29 Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany

30 Laboratoire Leprince-Ringuet, École Polytechnique, CNRS/IN2P3, 91128 Palaiseau, France

31 APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Cité, 10, rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France

32 Univ. Grenoble Alpes, IPAG, 38000 Grenoble, France
CNRS, IPAG, 38000 Grenoble, France

33 Department of Physics and Astronomy, The University of Leicester, University Road, Leicester, LE1 7RH, UK

34 Nicolaus Copernicus Astronomical Center, ul. Bartycka 18, 00-716 Warsaw, Poland

35 Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam, Germany

36 Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, 91058 Erlangen, Germany

37 DESY, 15738 Zeuthen, Germany

38 Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland

39 Centre for Astronomy, Faculty of Physics, Astronomy and Informationics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

40 Department of Physics, University of the Free State, PO Box 339, 9300 Bloemfontein, South Africa

41 Heisenberg Fellow (DFG), ITA Universität Heidelberg, Germany

42 GRAPPA, Institute of High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

43 Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-ku, 171-8501 Tokyo, Japan

44 Now at Santa Cruz Institute for Particle Physics and Department of Physics, University of California at Santa Cruz, Santa Cruz, CA 95064, USA

45 Now at: Technische Universität Kaiserslautern, Fachgebiet Physik/Energetische Gebäudeoptimierung, Paul-Ehrlich-Str., Gebäude 29, Raum 214, 67663 Kaiserslautern, Germany

46 Now at: Institute for Space Sciences (CSIC-IEEC), Campus UAB, Carrer de Can Magrans s/n, 08193 Barcelona, Spain

47 Now at: Karlsruhe Institute of Technology, Institute for Nuclear Physics, PO Box 3640, 76021 Karlsruhe, Germany
# 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 $\delta t$. In each step, the amount of spin-down energy converted into relativistic electrons and positrons is given as

$$ \Delta E_p(t) = \eta \int_t^{t+\delta t} \dot{E}(t') dt', \qquad (A.1) $$

where the spin-down evolution $\dot{E}(t)$ of the pulsar (Eq. (3), p. 7) is characterised by the braking index $n$, the initial spin-down timescale $\tau_0$, and the initial spin-down $\dot{E}_0$. The lepton conversion efficiency $\eta$ 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 $\tau_0$ to one on $P_0$ using

$$ \tau_0 = \frac{2\tau_c}{n-1} \left(\frac{P_0}{P}\right)^{n-1}, \qquad (A.2) $$

in this work we take $\tau_0$ 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:

$$ \frac{dN_{\text{inj}}}{dE}(E, t) = \Phi_0(t) \left(\frac{E}{1 \text{ TeV}}\right)^{-\beta} \qquad (A.3) $$

with a power-law index $\beta$. $\Phi_0(t)$ can be calculated imposing

$$ \Delta E_p(t) \stackrel{!}{=} \int_{E_{\min}}^{E_{\max}} \frac{dN_{\text{inj}}}{dE}(E, t) dE. \qquad (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{dN_{\text{cooled}}}{dE}(E, t) = \frac{dN}{dE}(E, t - \delta t) \exp\left(-\frac{\delta t}{\tau_{\text{eff}}(E, t)}\right), \qquad (A.5) $$

with an effective cooling timescale

$$ \tau_{\text{eff}}^{-1} = \tau_{\text{syn}}^{-1} + \tau_{\text{esc}}^{-1} + \tau_{\text{ad}}^{-1}, \qquad (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)

$$ \tau_{\text{syn}}(E, t) = 12.5 \left[ \frac{B(t)}{10 \mu \text{G}} \right]^{-2} \left[ \frac{E}{10 \text{ TeV}} \right]^{-1} \text{kyr} \qquad (A.7) $$

$$ \tau_{\text{esc}}(E, t) = 34 \left[ \frac{B(t)}{10 \mu \text{G}} \right] \left[ \frac{E}{10 \text{ TeV}} \right]^{-1} \times \left[ \frac{R(t)}{1 \text{ pc}} \right]^2 \text{kyr}. \qquad (A.8) $$

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}{E_p}$, is governed by the expansion of the nebula and can be calculated (following de Jager & Harding 1992) from

$$ \frac{dE_{\text{ad}}}{dt} = -\frac{E}{3} \nabla v_{\perp}(R) = \dot{E}_p, \qquad (A.9) $$

with $v_{\perp}(R)$ being the radial component of the particle velocity. In general, its divergence can be calculated to

$$ \nabla v_{\perp}(R) = \frac{1}{R^2} \frac{\partial(R^2 v_{\perp})}{\partial R} \qquad (A.10) $$

$$ = \frac{1}{R(t)^2} \frac{\partial(R(t)^2 v_{\perp}(t))}{\partial t} \frac{\partial t}{\partial R} \qquad (A.11) $$

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<fnref n="12" /> distinct phases, which define the expansion behaviour of the PWN according to the age of the system in terms of the spin-down timescale $\tau_0$ and the reverse-shock interaction time $t_{\text{rs}}$. Usually, the reverse-shock passage and the subsequent reverberations are expected to occur at a time $t_{\text{rs}} > \tau_0$. For this case, the following relations have been derived in the aforementioned works:

$$ R(t) \propto \begin{cases} t^{6/5} & \text{for } t \le \tau_0 \\ t & \text{for } \tau_0 < t \le t_{\text{rs}} \\ t^{3/10} & \text{for } t > t_{\text{rs}}. \end{cases} \qquad (A.12) $$

<footnote n="12">¹² The original version of the model presented in Mayer et al. (2012) does not incorporate a free expansion phase and uses only two evolutionary stages.</footnote>
<table><caption>**Table A.1.** Overview of parameters used for the modelling and the calculation of *varied* model ranges.</caption><thead><tr><th rowspan="2">Parameter description</th><th rowspan="2"></th><th colspan="2">Parameter values</th></tr><tr><th>baseline model</th><th>varied model</th></tr></thead><tbody><tr><td>Braking index</td><td>n</td><td>3.0</td><td>2.5 ... 3.5</td></tr><tr><td>Initial spin-down power</td><td>˙E<sub>0</sub></td><td>(10<sup>39</sup> erg s<sup>-1</sup>)</td><td>1.0 ... 4.0</td></tr><tr><td>Initial spin-down timescale</td><td>τ<sub>0</sub></td><td>(kyr)</td><td>0.32 ... 0.77</td></tr><tr><td>Initial magn. field strength</td><td>B<sub>0</sub></td><td>(μG)</td><td>200 ... 110 ... 270</td></tr><tr><td>Reverse shock interaction timescale</td><td>t<sub>rs</sub></td><td>(kyr)</td><td>4.0 ... 8.0</td></tr><tr><td>PWN radius at t = 3 kyr</td><td>R<sub>3</sub></td><td>(pc)</td><td>6.0 ... 3.0 ... 12.0</td></tr><tr><td>Adopted const. ISM magn. field strength</td><td>B<sub>ISM</sub></td><td>(μG)</td><td>3.0</td></tr><tr><td>Lepton conversion efficiency</td><td>η</td><td></td><td>1.0</td></tr><tr><td>Index of magn. field evolution</td><td>α</td><td></td><td>0.6</td></tr><tr><td>Index of lepton injection spectrum</td><td>β</td><td></td><td>2.0</td></tr><tr><td>Lower bound of lepton energy distribution</td><td>E<sub>min</sub></td><td>(TeV)</td><td>0.03</td></tr><tr><td>Upper bound of lepton energy distribution</td><td>E<sub>max</sub></td><td>(TeV)</td><td>300</td></tr></tbody></table>

In the (supposably much less common) opposite case, $t_{rs} < \tau_0$, the time evolution of $R_{PWN}$ is modified to

$$R(t) \propto \begin{cases} t^{6/5} & \text{for } t \le t_{rs} \\ t^{11/15} & \text{for } t_{rs} < t \le \tau_0 \\ t^{3/10} & \text{for } t > \tau_0. \end{cases} \quad (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_{ISM}, \quad (A.14)$$

assuming a constant and homogeneous ISM contribution of 3 μG, and adopting an index of $\alpha = 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 + \delta t$ is determined by the balance of freshly injected leptons and those cooled out of the respective energy interval,

$$\frac{dN}{dE}(E, t + \delta t) = \frac{dN_{\text{cooled}}}{dE}(E, t) + \frac{dN_{\text{inj}}}{dE}(E, t + \delta t). \quad (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 match the 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.

- $\tau_0$, $\dot{E}_0$: these parameters define the starting point of the pulsar trajectory on Fig. 2 and the total energy budget of the pulsar
<table><caption>**Table A.2.** Evolution of a PWN in our *baseline* model.</caption><thead><tr><th rowspan="2">t<br>(kyr)</th><th colspan="3">Pulsar</th><th colspan="3">PWN</th></tr><tr><th>τ<sub>c</sub><br>(kyr)</th><th>È<br>(10<sup>38</sup> erg s<sup>-1</sup>)</th><th>B<sub>PWN</sub><br>(μG)</th><th>R<sub>PWN</sub><br>(pc)</th><th>L<sub>1-10 TeV</sub><br>(10<sup>33</sup> erg s<sup>-1</sup>)</th><th>Γ</th></tr></thead><tbody><tr><td>0.10</td><td>0.60</td><td>1.39 × 10<sup>39</sup></td><td>148</td><td>0.142</td><td>1.27 × 10<sup>35</sup></td><td>2.08</td></tr><tr><td>0.14</td><td>0.63</td><td>1.23 × 10<sup>39</sup></td><td>140</td><td>0.207</td><td>1.36 × 10<sup>35</sup></td><td>2.11</td></tr><tr><td>0.19</td><td>0.69</td><td>1.05 × 10<sup>39</sup></td><td>131</td><td>0.316</td><td>1.41 × 10<sup>35</sup></td><td>2.14</td></tr><tr><td>0.26</td><td>0.76</td><td>8.63 × 10<sup>38</sup></td><td>122</td><td>0.458</td><td>1.42 × 10<sup>35</sup></td><td>2.17</td></tr><tr><td>0.36</td><td>0.85</td><td>6.78 × 10<sup>38</sup></td><td>113</td><td>0.665</td><td>1.37 × 10<sup>35</sup></td><td>2.19</td></tr><tr><td>0.49</td><td>0.99</td><td>5.07 × 10<sup>38</sup></td><td>103</td><td>0.971</td><td>1.28 × 10<sup>35</sup></td><td>2.22</td></tr><tr><td>0.67</td><td>1.17</td><td>3.61 × 10<sup>38</sup></td><td>94.0</td><td>1.34</td><td>1.16 × 10<sup>35</sup></td><td>2.25</td></tr><tr><td>0.92</td><td>1.42</td><td>2.44 × 10<sup>38</sup></td><td>84.6</td><td>1.84</td><td>1.01 × 10<sup>35</sup></td><td>2.28</td></tr><tr><td>1.27</td><td>1.77</td><td>1.58 × 10<sup>38</sup></td><td>75.6</td><td>2.54</td><td>8.36 × 10<sup>34</sup></td><td>2.30</td></tr><tr><td>1.74</td><td>2.24</td><td>9.82 × 10<sup>37</sup></td><td>67.0</td><td>3.49</td><td>6.72 × 10<sup>34</sup></td><td>2.31</td></tr><tr><td>2.40</td><td>2.89</td><td>5.89 × 10<sup>37</sup></td><td>59.0</td><td>4.79</td><td>5.27 × 10<sup>34</sup></td><td>2.32</td></tr><tr><td>3.29</td><td>3.79</td><td>3.44 × 10<sup>37</sup></td><td>51.7</td><td>6.58</td><td>4.04 × 10<sup>34</sup></td><td>2.33</td></tr><tr><td>4.52</td><td>5.02</td><td>1.96 × 10<sup>37</sup></td><td>45.0</td><td>8.30</td><td>3.19 × 10<sup>34</sup></td><td>2.35</td></tr><tr><td>6.21</td><td>6.71</td><td>1.10 × 10<sup>37</sup></td><td>39.0</td><td>9.13</td><td>2.46 × 10<sup>34</sup></td><td>2.38</td></tr><tr><td>8.53</td><td>9.03</td><td>6.05 × 10<sup>36</sup></td><td>33.7</td><td>10.0</td><td>1.84 × 10<sup>34</sup></td><td>2.39</td></tr><tr><td>11.7</td><td>12.2</td><td>3.30 × 10<sup>36</sup></td><td>29.1</td><td>11.0</td><td>1.35 × 10<sup>34</sup></td><td>2.39</td></tr><tr><td>16.1</td><td>16.6</td><td>1.79 × 10<sup>36</sup></td><td>25.1</td><td>12.2</td><td>9.71 × 10<sup>33</sup></td><td>2.38</td></tr><tr><td>22.1</td><td>22.6</td><td>9.63 × 10<sup>35</sup></td><td>21.6</td><td>13.4</td><td>6.92 × 10<sup>33</sup></td><td>2.35</td></tr><tr><td>30.4</td><td>30.9</td><td>5.16 × 10<sup>35</sup></td><td>18.6</td><td>14.7</td><td>4.87 × 10<sup>33</sup></td><td>2.32</td></tr><tr><td>41.8</td><td>42.2</td><td>2.76 × 10<sup>35</sup></td><td>16.1</td><td>16.2</td><td>3.39 × 10<sup>33</sup></td><td>2.29</td></tr><tr><td>57.4</td><td>57.9</td><td>1.47 × 10<sup>35</sup></td><td>13.9</td><td>17.8</td><td>2.32 × 10<sup>33</sup></td><td>2.25</td></tr><tr><td>78.8</td><td>79.3</td><td>7.84 × 10<sup>34</sup></td><td>12.1</td><td>19.6</td><td>1.58 × 10<sup>33</sup></td><td>2.21</td></tr><tr><td>108</td><td>109</td><td>4.17 × 10<sup>34</sup></td><td>10.6</td><td>21.5</td><td>1.06 × 10<sup>33</sup></td><td>2.17</td></tr><tr><td>149</td><td>149</td><td>2.21 × 10<sup>34</sup></td><td>9.33</td><td>23.7</td><td>6.96 × 10<sup>32</sup></td><td>2.14</td></tr><tr><td>204</td><td>205</td><td>1.17 × 10<sup>34</sup></td><td>8.26</td><td>26.0</td><td>4.53 × 10<sup>32</sup></td><td>2.11</td></tr><tr><td>281</td><td>281</td><td>6.23 × 10<sup>33</sup></td><td>7.37</td><td>28.6</td><td>2.91 × 10<sup>32</sup></td><td>2.08</td></tr><tr><td>386</td><td>386</td><td>3.30 × 10<sup>33</sup></td><td>6.62</td><td>31.5</td><td>1.83 × 10<sup>32</sup></td><td>2.05</td></tr><tr><td>530</td><td>530</td><td>1.75 × 10<sup>33</sup></td><td>6.00</td><td>34.6</td><td>1.14 × 10<sup>32</sup></td><td>2.03</td></tr><tr><td>728</td><td>728</td><td>9.29 × 10<sup>32</sup></td><td>5.49</td><td>38.1</td><td>7.03 × 10<sup>31</sup></td><td>2.00</td></tr><tr><td>1000</td><td>1000</td><td>4.92 × 10<sup>32</sup></td><td>5.06</td><td>41.9</td><td>4.26 × 10<sup>31</sup></td><td>1.98</td></tr></tbody></table>

**Notes.** *t* is the true age of the pulsar, τ<sub>c</sub> its characteristic age, È its spin-down luminosity. B<sub>PWN</sub> is the magnetic field in the PWN, R<sub>PWN</sub> the PWN radius, L<sub>1-10 TeV</sub> the TeV luminosity, and Γ the gamma-ray index between 1 and 10 TeV.

(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 × 10<sup>49</sup> erg.

- B₀: the initial B-field is chosen such that, using Eq. (A.14) for its decay, Crab-like young PWNe have (present) fields B<sub>PWN</sub> ~ 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<sub>RS</sub>, 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.

- B<sub>ISM</sub>: 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<sub>min</sub>, E<sub>max</sub>: 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.

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
<figure><img src="image_10.png" /><figcaption>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 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.</figcaption></figure>

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 $\dot{E}$ 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 $\dot{E}$ and $\tau_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 $\Omega$ as

$$ \dot{E} = -k' \Omega^4. \qquad (B.1) $$

Since the angular momentum loss rate is

$$ \dot{J} = \frac{\dot{E}}{\Omega} = k' \Omega^3 \qquad (B.2) $$

it follows that the velocity loss rate is

$$ \dot{\Omega} = \frac{\dot{J}}{I} = -k \Omega^3, \qquad (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, \qquad (B.4) $$

which turns Eq. (B.1) into

$$ \dot{E} = -k' \Omega^{n+1}. \qquad (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}}. \qquad (B.6) $$

Using Eq. (B.5), and $P = 2\pi/\Omega$, and differentiating $P$ one obtains

$$ \dot{E}(t) = \dot{E}_0 \left(1 + \frac{t}{\tau_0}\right)^{-\frac{n+1}{n-1}} \qquad (B.7) $$

$$ P(t) = P_0 \left(1 + \frac{t}{\tau_0}\right)^{\frac{1}{n-1}} \qquad (B.8) $$

$$ \dot{P}(t) = \frac{P_0}{\tau_0(n-1)} \left(1 + \frac{t}{\tau_0}\right)^{-\frac{n-2}{n-1}}. \qquad (B.9) $$

The canonical formulae to calculate $\dot{E}$ and $\tau_c$ from $P$ and $\dot{P}$ then yield

$$ \dot{E}(t) = 4\pi^2 I \frac{P(t)}{P(t)^3} = \frac{4\pi^2 I}{\tau_0 P_0^2 (n-1)} \left(1 + \frac{t}{\tau_0}\right)^{-\frac{n+1}{n-1}} \qquad (B.10) $$

$$ \tau_c(t) \equiv \frac{P(t)}{2\dot{P}(t)} = \frac{n-1}{2}(t + \tau_0) \qquad (B.11) $$
(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$, $\tau_c$ is $\tau_0 (n - 1)/2$ and increases steadily.

For the relation of $\dot{P}$ and $P$, Eqs. (B.8) and (B.9) furthermore imply

$$ \dot{P}(P) = \frac{P_0}{\tau_0} \frac{1}{n-1} \left(\frac{P}{P_0}\right)^{2-n}, \qquad (B.12) $$

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

In order to see what happens if $\dot{E}$ is plotted against $\tau_c$, one has to resolve the dependency on $t$ to arrive at

$$ \dot{E}(\tau_c) = \dot{E}_0 \left( \frac{2}{n-1} \frac{\tau_c}{\tau_0} \right)^{-\frac{n+1}{n-1}}. \qquad (B.13) $$

Clearly, the evolution curve of a pulsar on the $\dot{E}-\tau_c$ diagram starts at a point $[\tau_0 (n - 1)/2, \dot{E}_0]$, which depends on $\tau_0$, $n$, and $\dot{E}_0$, 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 $\tau_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:

$$ E_{\text{dep}}(t) = \int_{0}^{t} \dot{E}(t')dt' \qquad (B.14) $$

$$ = \frac{n-1}{2} \dot{E}_0 \tau_0 - \dot{E}(t) \tau_c(t). \qquad (B.15) $$

For $t \to \infty$, $\dot{E}(t) \tau_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 $I \Omega_0^2/2$, the total rotational energy of the pulsar. Unfortunately, $n$, $\dot{E}_0$, and $\tau_0$ are three unknown initial properties of the pulsar, so it cannot be measured. Unlike that, the second term $\dot{E} \tau_c$, which represents the present budget of rotational energy, can be calculated from the measured $P$ and $\dot{P}$. The ordinary (low-aged) pulsar with the maximum present budget of energy is PSR J0537-6910 in N 157B, with $7.6 \times 10^{49}$ erg, which is a lower limit to the maximal initial rotational energies that can be reached.