TeV Halos are Everywhere: Prospects for New Discoveries

TL;DR

This paper discusses the widespread presence of TeV halos around pulsars, predicts future detections by HAWC and CTA, and emphasizes their significance for understanding high-energy astrophysical processes.

Contribution

It provides the first comprehensive prediction of the number of TeV halos detectable by upcoming gamma-ray observatories, highlighting their potential ubiquity and scientific importance.

Findings

Future HAWC and CTA will detect approximately 50-240 TeV halos.

The existing HESS catalog may contain 10-50 unidentified TeV halos.

TeV halos are crucial for studying pulsar evolution and cosmic ray sources.

Abstract

Milagro and HAWC have detected extended TeV gamma-ray emission around nearby pulsar wind nebulae (PWNe). Building on these discoveries, Linden et al. [1] identified a new source class -- TeV halos -- powered by the interactions of high-energy electrons and positrons that have escaped from the PWN, but which remain trapped in a larger region where diffusion is inhibited compared to the interstellar medium. Many theoretical properties of TeV halos remain mysterious, but empirical arguments suggest that they are ubiquitous. The key to progress is finding more halos. We outline prospects for new discoveries and calculate their expectations and uncertainties. We predict, using models normalized to current data, that future HAWC and CTA observations will detect in total $\sim$50--240 TeV halos, though we note that multiple systematic uncertainties still exist. Further, the existing HESS source catalog could contain $\sim$10--50 TeV halos that are presently classified as unidentified sources or PWN candidates. We quantify the importance of these detections for new probes of the evolution of TeV halos, pulsar properties, and the sources of high-energy gamma rays and cosmic rays.

Tables

Table 1. Table 1: Physical quantities for Geminga and its TeV halo.

Observed	${\dot{E}}^{𝒢}$ [erg s-1]	3.2$\times {10}^{34}$	Manchester et al. (2005)
	$d^{𝒢}$ [pc]	250	Verbiest et al. (2012)
	${\Phi }_{7TeV}^{𝒢}$ [TeV-1cm-2s-1]	4.87$\times$10-14	Abeysekara et al. (2017b)
Calculated	$F_{\mathrm{TeV}}^{𝒢}$ [cm-2s-1]	$3.5\times {10}^{-12}$
	$L_{\mathrm{TeV}}^{𝒢}$ [erg s-1]	1.1$\times$1032

Table 2. Table 2: The most important uncertainties in the number of TeV halos that are discussed in this work. For each uncertainty, we note an alternative model, and roughly indicate the effect that such a model would have on the predicted TeV halo source count.

Name of Uncertainty	Default	Alternative	Effect
pulsar population
$P_0$ distribution	Gaussian	Uniform	Increase, $\times 2$
$B_0$ distribution	$⟨{\log }_{10}{B}_0⟩$ = 12.65	$⟨{\log }_{10}{B}_0⟩$ = 13.10	Decrease, $\times 0.5$
$\gamma$-ray efficiency
$\dot{E}$ dependence	$L_{\gamma }\propto \dot{E}$	$L_{\gamma }\propto {\dot{E}}^{0.8}$	Decrease, $\times 0.5$
		$L_{\gamma }\propto {\dot{E}}^{1.2}$	Increase, $\times 2$
Age dependence	$L_{\gamma }/\dot{E}$ = const.	$L_{\gamma }/\dot{E}\propto {\left(T_{\mathrm{age}}\right)}^{0.5}$	Decrease, $\times 0.3$
		$L_{\gamma }/\dot{E}\propto {\left(T_{\mathrm{age}}\right)}^{-0.5}$	Increase, $\times 3$
Source-to-source scatter	None	${\log }_{10}\left(L_{\gamma }/\dot{E}\right)\sim N(1,{0.5}^2)$ (lognormal, $\sigma =0.5$)	Increase, $\times 2$

Equations

$$\dot{E}(t)=\frac{8{\pi }^4{B}_0^2{R}^6}{3c^3{P}_0^4}{\left(1+\frac{t}{{\tau }_{\mathrm{sd}}}\right)}^{-2},$$

$$f_{\mathrm{beam}}=\left[9{\left({\log }_{10}\frac{P}{10\mathrm{s}}\right)}^2+3\right]\%.$$

$${\Phi }_{7TeV}={\Phi }_{7TeV}^{\mathcal{G}}\left(\frac{\dot{E}}{{\dot{E}}^{\mathcal{G}}}\right){\left(\frac{d^{\mathcal{G}}}{d}\right)}^2.$$

Full text

I Introduction

Milagro observations revealed extended TeV $\gamma$-ray emission surrounding the nearby Geminga pulsar, now confirmed by the High Altitude Water Cherenkov (HAWC) observatory Abdo et al. (2009); Abeysekara et al. (2017a, b). Additionally, HAWC has detected similar emission surrounding another nearby pulsar, PSR B0656+14, commonly associated with the Monogem ring Thorsett et al. (2003), and which we refer to as the “Monogem pulsar.” These sources are bright ($\sim 10^{32}$ erg s*-1*), have hard spectra ($\sim E^{-2.2}$), and are spatially extended ($\sim$ 25 pc). In addition, the High Energy Stereoscopic System (HESS) has detected a number of TeV $\gamma$-ray sources coincident with pulsars or pulsar wind nebulae (PWNe) Abdalla et al. (2018a, b). Though they refer to these as “TeV PWN,” they find that many are significantly larger than expected from PWN theory Linden et al. (2017); Khangulyan et al. (2018); Gaensler and Slane (2006). The sources noted above appear morphologically and dynamically distinct from PWNe detected in X-ray and radio observations.

Linden et al. (2017) identified these sources as a new $\gamma$-ray source class (“TeV Halos”) and interpreted their emission as the result of electrons and positrons interacting with the ambient interstellar radiation field outside the PWN. The possibility of significantly extended leptonic emission was first predicted in Ref. Aharonian (1995), and its importance was further discussed in Refs. Aharonian et al. (1997); Aharonian (2004); Yüksel et al. (2009); Aharonian (2013). Moreover, Linden et al. (2017) showed that a large fraction of 2HWC catalog sources are coincident with pulsars, and predicted that TeV halos are a generic feature of pulsar emission.

In Fig. 1, we show how a TeV halo compares to other features at the site of a past core-collapse supernova explosion. For a given source, it may be that not all components are detectable or even present at the same time. A PWN, powered by the rotational energy of the central pulsar, is delimited by the contact discontinuity between the shocked pulsar wind and the ejecta or interstellar matter. An SNR, powered by the energy of the supernova explosion, is delimited by its interaction with the interstellar medium. A TeV halo is likely intermediate in size, is powered by cosmic rays diffusing away from the PWN, and does not have a well-defined boundary. The size of a PWN can be on the order of 0.1–1 pc, though some may range up to $\sim$10 pc Gaensler and Slane (2006); Kargaltsev and Pavlov (2008), and the size of an SNR may span $\sim$1–100 pc Badenes et al. (2010); Stafford et al. (2018), depending on their properties, evolutionary stages and environment. The typical size of a TeV halo is not known, but Geminga and Monogem observations indicate that it may be on the order of 10 pc for middle-aged pulsars. For the three types of object, differences in radii lead to larger differences in volumes that further support different physical origins.

The identification of TeV halos as a new source class is supported by the subsequent detection of two more TeV halos by HAWC Riviere et al. (2017); Brisbois et al. (2018), one of which was predicted by Ref. Linden et al. (2017). However, many details about TeV halos remain unknown and further observations will have broad implications. Apart from shedding light on the properties of the TeV halos themselves, these observations will reveal new aspects of pulsar formation and evolution Linden et al. (2017), and will probe sources of high-energy $\gamma$-rays Hooper et al. (2018); Linden and Buckman (2018); Hooper and Linden (2018a) and cosmic-ray electrons and positrons Hooper et al. (2017); Hooper and Linden (2018b); Evoli et al. (2018); Fang et al. (2018); Profumo et al. (2018); Bucciantini (2018); Tang and Piran (2018); Xi et al. (2018).

Here, we outline a multifaceted strategy to discover more TeV halos and to constrain their evolution. We quantify the role of Galactic source searches and diffuse measurements using water Cherenkov telescopes like HAWC. We also consider Galactic and extragalactic source surveys by imaging air Cherenkov telescopes, focusing on the Cherenkov Telescope Array (CTA). Further, we show that follow-up studies of existing TeV $\gamma$-ray sources in the HESS catalog, especially those classified as PWN or unidentified sources, could find more TeV halos. For our overall approach, we use standard methods for pulsar population synthesis and treat the Geminga TeV halo as a prototype. Our results go significantly beyond those of prior work, yielding new insights into both the prospects for future TeV halo discoveries, and the implications of TeV halo observations for our understanding of astrophysics.

In Sec. II, we briefly review the properties of TeV halos. In Sec. III, we present our methods for modeling TeV halo populations. In Sec. IV, we compare predictions with current observations and constrain model parameters. In Sec. V, we outline future directions to find more TeV halos. In Sec. VI, we present our conclusions.

II What are TeV halos?

TeV halos are defined as the non-thermal emission produced in regions outside a PWN, but within a region where pulsar activity dominates cosmic-ray diffusion (Linden et al. (2017); Aharonian (1995); Aharonian et al. (1997); Aharonian (2004); Yüksel et al. (2009)). Within this region, multi-TeV $\gamma$-rays are produced by the inverse-Compton scattering of ambient photons by $\sim$10 TeV electrons and positrons accelerated by the pulsar wind termination shock. Observations indicate that the TeV halo produces bright $\gamma$-ray emission with a hard spectrum.

We begin by examining the key features of the best-studied TeV halo, Geminga, which is about 340 kyr old Manchester et al. (2005) and believed to reside approximately 250 pc from Earth Verbiest et al. (2012). HAWC detects TeV $\gamma$-ray emission extending to an angular size of $\sim$5∘, corresponding to $\sim$25 pc in physical extent Abeysekara et al. (2017a). The differential $\gamma$-ray luminosity at 7 TeV is 2.9\times 10^{31}$$(d/250\ {\rm pc})^{2} erg s*-1*, with a local spectral index of $-$2.2.

The lack of gas-correlated emission indicates a leptonic origin. Within the context of an inverse-Compton model, several parameters regarding the electron population can be calculated Hooper et al. (2017); Abeysekara et al. (2017a). To produce the bright $\gamma$-ray luminosity, $\sim$10% of the total pulsar spin-down power must be converted into e± pairs. Furthermore, to produce the hard $\gamma$-ray spectrum, the electron population should be injected with a hard power-law spectrum between \sim$$-1.5 and $-$2.2.

The most notable feature of TeV halos is their size. The Geminga TeV halo is significantly larger than its X-ray PWN, which is confined within $3^{\prime}$ of the central pulsar Posselt et al. (2017). This indicates that the electrons and positrons responsible for TeV halo emission have already escaped the PWN and are interacting with the interstellar radiation field. The TeV halo morphology is consistent with cosmic-ray diffusion, rather than advection Abeysekara et al. (2017a). However, this diffusion must be inhibited. Assuming that the TeV halo medium is filled with the $\sim$1 eV cm*-3* interstellar radiation field and the $\sim$3 $\mu$G magnetic field typical of the Galactic Plane, we would expect 10 TeV e± to cool in $\sim$40 kyr. In the interstellar medium, electrons and positrons that propagate for $\sim$40 kyr should diffuse over a distance of $\sim$700 pc Trotta et al. (2011). However, the TeV halo power appears to be confined within $\sim$25 pc of the pulsar center.

TeV halo emission is not unique to Geminga. The HAWC collaboration has identified at least three other TeV halos with similar features: Monogem (111 kyr, 290 pc), PSR B0540+23 (253 kyr, 1.56 kpc), and PSR J0633+0632 (59 kyr, 1.35 kpc) Manchester et al. (2005); Riviere et al. (2017); Brisbois et al. (2018). In addition, Linden et al. (2017) listed 13 more TeV halo candidates in the 2HWC catalog. The 2HWC survey also provides a hint of TeV halo emission around millisecond pulsars Hooper and Linden (2018a).

In addition to HAWC, imaging air Cherenkov telescopes like HESS, MAGIC and VERITAS have detected a number of extended TeV $\gamma$-ray sources that are associated with pulsars or PWNe observed at other wavelengths. These systems are called “TeV PWN,” but many of them have an extension exceeding $\sim$10 pc Abdalla et al. (2018a, b), while hydrodynamical simulations predict a typical PWN size on the order of 1 pc van der Swaluw et al. (2004); Gaensler and Slane (2006); Gelfand et al. (2009); Bucciantini (2011); Martín et al. (2016); Ishizaki et al. (2018); Zhu et al. (2018); van Rensburg et al. (2018). More pointedly, they are usually much more extended than the size of X-ray PWN observed from the same system Kargaltsev and Pavlov (2010); Kargaltsev et al. (2013). These observations suggest that some of these $\gamma$-ray sources may be interpreted as TeV halos, instead of emission from confined particles inside PWNe. In particular, HESS J1825-137 has the largest radius ($\sim$50 pc Abdalla et al. (2019)) among “TeV PWN” Khangulyan et al. (2018). A TeV halo explanation for this source is already discussed in Refs. Aharonian (2004, 2013).

Despite the significant number of TeV halos that have been (or are potentially) detected in current surveys, many of their properties remain mysterious. In particular, we do not understand the evolution of the key observable TeV halo properties: their luminosity, spectrum, and spatial morphology. In recent work Linden et al. (2017), TeV halo predictions have been evaluated utilizing a “Geminga-like” model, where the ratio of the $\gamma$-ray flux to $\dot{E}/d^{2}$ is constant for all systems with an efficiency set to the best-fit value of Geminga, and the physical size of all TeV halos is $\sim$10 pc. On one hand, this model appears reasonably consistent with the data — choosing to instead normalize the TeV halo flux to the average $\gamma$-ray efficiency of all firmly identified TeV halos changes the normalization constant only by a factor of $\sim$2 compared to the “Geminga-like” model. On the other hand, there is nearly an order of magnitude variation in the efficiencies of individual candidate sources, the origin of which is not understood.

A key question is when a TeV halo first forms. Several considerations indicate that the “Geminga-like” model may not apply to young pulsars. Theoretically, high-energy cosmic rays are expected to be efficiently confined in young PWNe and quickly lose energy to adiabatic and synchrotron cooling in the strong PWN magnetic field Rees and Gunn (1974); Kennel and Coroniti (1984); de Jager et al. (2009); Tanaka and Takahara (2010); Torres et al. (2014); Vorster et al. (2013); Olmi et al. (2016). This may imply that particles do not escape into young TeV halos. Moreover, the creation of a halo may require cosmic-ray self-generated turbulence, which is produced through the resonant interactions of Alfvén waves with accelerated electrons and positrons. The growth-rate of self-generated turbulence is model dependent, but typically occurs on $>$ kyr timescales Evoli et al. (2018).

Observationally, the Crab pulsar (964 yr, 2 kpc) does not appear to produce TeV halo emission Aharonian et al. (2006); Albert et al. (2008); Abeysekara et al. (2017c), indicating that TeV halos may not be visible within the first kyr of pulsar evolution. An intriguing edge case is the Vela pulsar (11 kyr, 280 pc). Vela does not appear to produce a bright TeV halo (compared to the luminosity expected if the formation efficiency is Geminga-like). However, Vela does have dim, spatially-extended emission detected in radio and GeV-TeV $\gamma$-ray observations Abdo et al. (2010); Abramowski et al. (2012); Grondin et al. (2013); Tibaldo et al. (2018). This has historically been interpreted as a class of “relic PWN” that are left behind after the interaction of the expanding PWN and the SNR reverse shock, and which are powered by old electrons accumulated since the birth of the pulsar Blondin et al. (2001); de Jager et al. (2008); Hinton et al. (2011); Slane et al. (2018). Interestingly, the size of this extended emission is $\sim$10 pc, comparable to that of observed TeV halos. Thus, Vela could be interpreted as a transition case, where inefficient TeV halos first form. Further TeV observations around $\sim$1–10 kyr pulsars are needed to study the properties of young systems.

Because detailed examinations of young systems are beyond the scope of this study, in this paper we use a standard “Geminga-like” model, but introduce a new parameter $T_{\rm min}$, before which pulsars are assumed to exhibit no TeV halo activity. Observations of the Crab and Vela suggest $T_{\rm min}\gtrsim$ 1–10 kyr, while Monogem and Geminga TeV halos suggest $T_{\rm min}\lesssim$ 100–300 kyr.

Another key question is whether TeV halo activity is ubiquitous to all pulsars. Theoretically, the creation of halos requires strongly inhabited diffusion around pulsars, which might be expected for all pulsars if a steep cosmic-ray gradient around them efficiently excites self-generated turbulence Evoli et al. (2018). Observationally, Linden et al. (2017) listed seven middle-aged pulsars that should be detected by HAWC, under the assumption that every pulsar has a Geminga-like TeV halo, and find that five are in fact associated with the 2HWC sources. These are consistent with the expectation that a significant fraction of pulsars have TeV halos.

We operate under the assumption that all pulsars older than $T_{\rm min}$ produce TeV halos. This can be tested in future surveys. We do not consider the maximum age of TeV halos, because late-time sources are not important due to their small spindown power.

III TeV Halo Population Models

To model the population of TeV halos, we generate an ensemble of pulsars with randomly assigned initial spin periods ($P_{0}$) and magnetic fields ($B_{0}$). We assume $M$ = 1.4 $M_{\odot}$ and $R$ = 12 km for all pulsars Lattimer and Prakash (2007). We then assign each pulsar an age ($T_{\rm age}$) drawn from a uniform distribution spanning from 0 to 1 Gyr, and calculate the pulsar spindown power as:

[TABLE]

where $\tau_{\rm sd}=3Ic^{3}P_{0}^{2}/4\pi^{2}B_{0}^{2}R^{6}$ is the spindown timescale Shapiro and Teukolsky (1983); Hooper et al. (2018). We associate each pulsar with a randomly distributed position within the Milky Way, based on the pulsar distributions determined by Refs. Yusifov and Küçük (2004); Porter et al. (2017); Lorimer et al. (2006). Specifically, we adopt the radial distribution of Ref. Yusifov and Küçük (2004), and a scale height of 200 pc Porter et al. (2017), and calculate the pulsar position relative to Earth assuming a galactocentric distance of 8.5 kpc. Our results are only slightly affected if we use the alternative spatial distributions defined in Ref. Lorimer et al. (2006). We have verified that our models are reasonably consistent with the observations of nearby neutron stars, i.e., the seven isolated neutron stars and several pulsars younger than 1 Myr within around 500 pc Haberl (2005); Kaplan and van Kerkwijk (2009). Further, we calculate the probability that the pulsed radio emission from each pulsar is beamed towards Earth following the empirical relation defined in Ref. Tauris and Manchester (1998),

[TABLE]

Of all choices in our calculation, the most significant are those of $P_{0}$ and $B_{0}$, due to their strong dependence in Eq. 1: $B_{0}^{2}/P_{0}^{4}$ in the pre-factor and $P_{0}^{2}/B_{0}^{2}$ in $\tau_{\rm sd}$.

The $P_{0}$ distribution is poorly constrained Faucher-Giguère and Kaspi (2006); de Jager (2008); Watters and Romani (2011); Popov and Turolla (2012); Noutsos et al. (2013); Igoshev and Popov (2013); Cieślar et al. (2018); Gonthier et al. (2004); Gullón et al. (2014), because population statistics are not sensitive to it Gonthier et al. (2004); Gullón et al. (2014). Conventionally, pulsar population models adopt a Gaussian distribution with $\langle P_{0}\rangle$ = 300 ms and $\sigma_{P_{0}}$ = 150 ms, based on radio observations Faucher-Giguère and Kaspi (2006). However, studies of the $\gamma$-rays pulsar population hint at much smaller values $\langle P_{0}\rangle$ = 50 ms and $\sigma_{P_{0}}$ = 50/$\sqrt{2}$ ms Watters and Romani (2011). In what follows, we present results for both $P_{0}$ distributions. We also test an intermediate case of $\langle P_{0}\rangle$ = 120 ms and $\sigma_{P_{0}}$ = 60 ms. Finally, in Appendix A, we examine models that utilize a uniform, rather than a Gaussian distribution, for the initial spin period. In all cases, we set a minimum spin period at the Newtonian centrifugal breakup limit, $P_{0,\min}$ = $0.85$ $(M/1.4\ M_{\odot})^{1/2}$ $(R/12\ {\rm km})^{3/2}$ ms Lattimer and Prakash (2007).

For the $B_{0}$ distribution, we adopt a log-normal magnetic field distribution with mean $\langle\log_{10}B_{0}\rangle$ = 12.65 and standard deviation $\sigma_{\log_{10}B_{0}}$ = 0.55, which is derived from population studies of radio pulsars Faucher-Giguère and Kaspi (2006). Other studies predict magnetic fields that are larger by a factor of 2–4 Popov et al. (2010); Gullón et al. (2014). This uncertainty is discussed in Sec. IV.3. We do not include any term to account for the decay of the magnetic field strength, because it occurs on timescales of $>$Myr, much greater than the age of the majority of detectable TeV halos.

In Fig. 2, we show the evolution of the spindown power for six representative pulsars. Most of the integrated spindown power is spent before $\sim\tau_{\rm sd}$, which is 4 kyr for $P_{0}$ = 50 ms and 160 kyr for 300 ms in the fiducial case of $B_{0}$ = $10^{12.65}$ G (this would be more evident if we had plotted the power per log time, which would include multiplying by a factor $t$).

Using the ensemble of pulsars generated above, we assign a “Geminga-like” TeV halo to each pulsar that has an age older than $T_{\rm min}$. We normalize the differential $\gamma$-ray flux at 7 TeV ($\Phi_{\rm 7TeV}$) for each TeV halo relative to Geminga, using the spindown power ($\dot{E}$) and distance ($d$) as

[TABLE]

We adopt physical quantities for Geminga (superscript “$\mathcal{G}$”) as summarized in Table 1.

In addition to directly observable parameters such as the spindown energy and the 7-TeV $\gamma$-ray flux, our models also require us to derive parameters such as the integrated $\gamma$-ray flux ($F_{\rm TeV}$) and luminosity ($L_{\rm TeV}$) for each pulsar; we calculate these above 1 TeV. To normalize these parameters to Geminga, we follow the theoretical treatment of Ref. Linden and Buckman (2018), which calculates the inverse-Compton scattering $\gamma$-ray spectrum from Ref. Blumenthal and Gould (1970). We model the electron spectrum following Ref. Hooper et al. (2017), which derives the best-fit $\gamma$-ray spectrum from a combination of HAWC (7 TeV) and Milagro (35 TeV) observations of the Geminga TeV halo Abeysekara et al. (2017b); Abdo et al. (2009). Specifically, we assume that electrons are injected with a power-law index of $\alpha$ = 1.9 that cuts off exponentially at $E_{\rm cut}$ = 49 TeV. We adopt an energy-independent escape time of 1.8$\times$104 yr from the TeV halo emission region Hooper et al. (2017). The total e± luminosity is normalized to be $\eta\dot{E}$. We find the best-fit value of $\eta=0.12$ from the observed $\gamma$-ray flux. The derived values of $F_{\rm TeV}^{\mathcal{G}}$ and $L_{\rm TeV}^{\mathcal{G}}$ are provided in Table 1. We again calculate $F_{\rm TeV}$ and $L_{\rm TeV}$ for every other pulsar by scaling the best-fit Geminga values with $\dot{E}$ and $d$ as shown in Eq. (3).

The 2HWC catalog reports the photon index at 7 TeV (2.23 for Geminga). If we extrapolate this spectral index down to 1 TeV and use the 7-TeV differential flux, we derive values of $F_{\rm TeV}^{\mathcal{G}}$ that fall within $\sim$20$\%$ of the theoretically derived photon flux reported in Table 1. On the other hand, the values for $L_{\rm TeV}^{\mathcal{G}}$ are increased nearly by a factor of 2, and hence our calculated luminosity in Table 1 may be pessimistic.

In Fig. 3, we show the $\gamma$-ray luminosity function of Milky Way TeV halos for two different $P_{0}$ distributions and three different values of $T_{\rm min}$. We normalize the total number of Milky Way pulsars using a pulsar birth rate of 0.015 yr*-1* Lorimer et al. (2006). The upper panel shows the number weighting only, while the lower panel also includes the luminosity weighting.

If we do not set $T_{\rm min}$, the bright end of the number count (upper panel) has a slope of $\sim L^{-0.8}$, which is driven by the distribution of $P_{0}$. As pulsars get older (above $\tau_{\rm sd}$), they lose spindown power following $\dot{E}\propto t^{-2}$ and move to the left in this plot, producing a shallower slope of $\sim L^{-0.5}$ before the peak, where pulsars with average properties ($P_{0},B_{0},T_{\rm age}$) reside. We do not include the effect that old pulsars may terminate their activities below the radio death line Chen and Ruderman (1993), because late-time sources have small $\gamma$-ray luminosities and contribute negligibly to the source count. Furthermore, due to the shallow slopes of the number count, dim sources contribute negligibly to the total Galactic emission, as shown in Fig 3 (bottom).

IV Existing Model Constraints

The range of model parameters used in our predictions below can be constrained by current data. In Sec. IV.1, we predict the number of TeV halos that should be detected in the 2HWC source catalog. In Sec. IV.2, we estimate the contribution of unresolved TeV halos to the diffuse TeV $\gamma$-ray emission across the Galactic Plane, comparing our predictions with Milagro observations. In Sec. IV.3, we summarize model constraints and briefly discuss uncertainties.

IV.1 Sources in the 2HWC Catalog

The 2HWC catalog utilizes 507 days of HAWC data and identifies 39 high-significance sources within the field of view of $-20^{\circ}<$ decl. $<60^{\circ}$. The sensitivity depends on the photon spectral index and the source declination. We adopt the average of quoted values for spectral indices of $-2.5$ and $-2.0$ and the declination dependence given in Ref. Abeysekara et al. (2017b). The best sensitivity of 4.3$\times 10^{-15}$ TeV*-1* cm*-2* s*-1* (9% of the Geminga TeV halo flux) occurs at a declination of 20∘, and is degraded by a factor of $\sim$2 for declinations that differ by 30∘. We take into account the degradation of the flux sensitivity for sources that are larger than the size of the PSF, utilizing a model where the sensitivity decreases by a factor of $\theta_{\rm size}/\theta_{\rm PSF}$ compared to the point source sensitivity Hinton and Hofmann (2009). We assume a PSF size of $\theta_{\rm PSF}=0.2^{\circ}$ for HAWC Abeysekara et al. (2013). To determine the source size, we again utilize a Geminga-like model, assuming that all TeV halos have the same physical size as that of the Geminga halo ($\theta_{\rm size}=2^{\circ}$ at a distance of 250 pc). We ignore source confusion, where HAWC may identify neighboring or overlapping sources as one source, because our calculations show it to be unimportant.

We constrain our TeV halo models by requiring that they do not produce too many or too few systems that would be detected in the 2HWC catalog search. We set the maximum number of potential TeV halos in the 2HWC catalog at 36, because three sources (the Crab, Mrk501, and Mrk421) are associated with objects that are definitively not TeV halos. For the minimum number, we choose 2, because Geminga and Monogem were detected while the two other sources announced by Astronomer’s Telegrams Riviere et al. (2017); Brisbois et al. (2018) did not meet the flux threshold to be included in the 2HWC catalog. Both of these choices are conservative, as they do not take into account additional information concerning individual 2HWC objects.

We can additionally constrain the number of detectable TeV halos that would have radio beams oriented towards Earth. Such sources are especially compelling because the spatial coincidence points towards a TeV halo origin. We note that while the Monogem pulsar is a firmly detected radio pulsar, the Geminga pulsar has extremely dim radio emission and would not have been detected in blind radio searches Malofeev and Malov (1997). Hence, we conservatively assume that at least 1 TeV halo (Monogem) has been detected in the 2HWC catalog with a radio beam oriented towards Earth.

The lower limits on the number of beamed and unbeamed TeV halos would become much stronger if the TeV halo candidates that were first identified by Ref. (Linden et al., 2017) are confirmed by subsequent observations. Ref. (Linden et al., 2017) finds three additional 2HWC sources that are consistent with the position of middle-aged radio pulsars, and twelve additional 2HWC sources that are consistent with the positions of younger pulsars. They estimate that only 2.6 chance coincidences would be expected if the 2HWC sources were not associated with pulsar activity. We compare our model predictions with these candidate sources.

In Fig. 4 (top), we show the total number of detectable TeV halos produced by our model, regardless of whether the system has a radio beam that is oriented towards Earth. This prediction should thus be compared to the total number of detected 2HWC sources. We plot results for models with initial pulsar spin periods of $\langle P_{0}\rangle$ = 50 ms, 120 ms, and 300 ms, and find that our parameters allow us to vary the predicted number of detected halos by about one order of magnitude.

This variation translates into a constraint on the value of $T_{\rm min}$. If typical pulsars are born with relatively large spin periods (e.g., $\langle P_{0}\rangle$ = 300 ms), TeV halos produce about 10 sources in models where $T_{\rm min}$ = 0. In this case, the lower bound of $T_{\rm min}$ is not strongly constrained by 2HWC data. On the other hand, if we utilize our minimum value of $\langle P_{0}\rangle$ = 50 ms, TeV halos produce about 100 sources in models where $T_{\rm min}$ = 0. Because this exceeds the total number of 2HWC sources, this would require a simultaneous constraint of $T_{\rm min}$ $\gtrsim$ 50 kyr. In the remainder of the section, we adopt $T_{\rm min}$ = 50 kyr for the case of $\langle P_{0}\rangle$ = 50 ms as the most optimistic case, which predicts that most of the 39 sources in the 2HWC catalog are TeV halos. Models with $\langle P_{0}\rangle$ = 120 ms provide a critical case, approximately saturating the number of detectable TeV halos in models with $T_{\rm min}$ =0. Thus, in this case, the value of $T_{\rm min}$ is not strongly constrained at this point, but may be better constrained if future observations indicate that several 2HWC sources are not TeV halos.

To produce at least two detectable TeV halos, we need to set $T_{\rm min}\lesssim$ 300 kyr for $\langle P_{0}\rangle$ = 120ms. This constraint is not strong, because such a large value of $T_{\rm min}$ is already disfavored by the observations of Monogem (110 kyr) and Geminga (340 kyr). On the other hand, for $\langle P_{0}\rangle$ = 300ms, we can constrain $T_{\rm min}\lesssim$ 70 kyr.

In Fig. 4 (bottom), we show model predictions for the expected number of TeV halos in the 2HWC catalog that have radio beams aligned with Earth, compared with the age distribution of these TeV halo candidates. We first focus on middle-aged pulsars ($>$100 kyr). The $\langle P_{0}\rangle$ = 50 ms model predicts $\sim$9 sources, which slightly exceeds the number of TeV halo candidate systems identified in Ref. (Linden et al., 2017). This model is allowed, but if future observations rule out the TeV halo nature of several of these systems, it would be in tension with the data.

On the other hand, models with $\langle P_{0}\rangle$ = 300 ms produce $\lesssim$ 1 middle-aged TeV halo with a radio beam directed towards Earth, which approximately saturates the lower limit produced by the identification of Monogem. This model is allowed, but if future observations confirm the TeV halo origin of candidate sources, it would be disfavored. Intriguingly, though we adopted models with $\langle P_{0}\rangle$ = 50 ms and $\langle P_{0}\rangle$ = 300 ms based on previous pulsar studies Faucher-Giguère and Kaspi (2006); Watters and Romani (2011), they coincidentally also serve as reasonable estimates for the largest and smallest values allowed by the 2HWC data. The firm interpretation of existing 2HWC observations could potentially rule out either model.

We additionally show an intermediate case, with $\langle P_{0}\rangle$ = 120 ms, which predicts the observation of $\sim 2$ middle-aged TeV halos with radio beams oriented towards Earth. This model matches current observations well, and is likely remain consistent regardless of the interpretation of the 2HWC candidate sources.

Expanding our analysis to include all TeV halos with aligned radio beams regardless of the TeV halo age (including young sources), we find that the interpretations become trickier because models predict fairly similar TeV halo number counts. Models with $\langle P_{0}\rangle$ = 300 ms produce $\sim$ 3 TeV halos in case that $T_{\rm min}$ = 0. Meanwhile, models with $\langle P_{0}\rangle$ = 50 ms predict approximately $\sim$15 sources for $T_{\rm min}$ = 50 kyr. Our intermediate model with $\langle P_{0}\rangle$ = 120 ms also predicts the observation of $\sim$15 sources for $T_{\rm min}$ = 0. However, the age distribution of observed TeV halos differs markedly between models with and without a firm value of $T_{\rm min}$. Thus, future observations that correlate TeV halo activity with pulsars of known ages can more clearly distinguish between models of TeV halo formation, even in light of degeneracies between $\langle P_{0}\rangle$ and $T_{\rm min}$.

In Fig. 5, we show the cumulative flux distribution of all TeV halos within the HAWC field of view. The 50-ms model with no $T_{\rm min}$ produces $\sim$20 sources that have $\gamma$-ray fluxes larger than that of Geminga, while the 2HWC catalogue contains 5–12 such potential sources (depending on source extension), so the prediction is somewhat too high. Furthermore, this model predicts a few sources that are at least an order of magnitude brighter than Geminga, while no such source is reported, so the prediction is again somewhat too high, though consistent with Poisson fluctuations. Therefore, while this model is not ruled out by the flux distribution, it is in slight tension. All of the other models are consistent with data.

Due to the steep slope at the bright end of the luminosity function (Fig. 3), nearby sources are expected to dominate the source count. Indeed, in our estimate, about 50$\%$ of observable sources are located within $\simeq$ 3 kpc from Earth. This suggests that many observed TeV halos might have large angular sizes, indicating the importance of HAWC, which is suited for detecting extended sources.

IV.2 Diffuse TeV Gamma-Ray Emission Measurements with Milagro

Milagro measured the diffuse Galactic $\gamma$-ray flux above 3.5 TeV, finding $\phi(>$3.5 TeV) = (6.8$\pm$1.5$\pm$2.2) $\times 10^{-11}$ cm*-2* s*-1* sr*-1* within a region spanning $40^{\circ}<l<100^{\circ}$ and $|b|<5^{\circ}$ Atkins et al. (2005). We constrain our model by requiring that unresolved TeV halos do not overproduce this flux. Due to the hard TeV halo spectrum, alternative diffuse emission measurements at lower energies by ARGO-YBJ Bartoli et al. (2015) or a higher energy in a smaller region analyzed by Milagro Abdo et al. (2008) give comparable constraints.

To estimate the contribution from TeV halos to this emission, we include contributions from unresolved TeV halos with fluxes below that of Geminga. We also include the contribution from electrons and positrons that escape from unresolved TeV halos and provide a diffuse emission component that fills the interstellar medium. We remove contributions from any individual halo with a $\gamma$-ray flux exceeding Geminga because such a source would be detected by Milagro Linden and Buckman (2018). Because the number of such sources is small ($\sim$1 or fewer) and particles lose a significant fraction of their energy in the halo region, their contribution to the diffuse emission is not important. This treatment also allows us to remove unrealistically bright individual halos predicted for the $T_{\rm min}=0$ model, as seen in Fig. 3. Note that we show our model predictions in cumulative contributions from pulsars above 1 kyr. Thus, the results for $T_{\rm min}=0$ are identical for any model $T_{\rm min}\leq$ 1 kyr, and are not affected by these unrealistically bright halos that only occur for $T_{\rm min}\simeq 0$. We have verified that the number of sources that contribute is large enough that the result is not subject to statistical fluctuations.

In Fig. 6, we show that the current diffuse measurement does not strongly constrain our models. However, we stress that having a more precise measurement in the future could provide complementary constraints to future source surveys.

The diffuse TeV $\gamma$-rays are particularly important, because, as first shown in Ref. Prodanović et al. (2007), the Milagro measurements Atkins et al. (2005) of the diffuse flux from the Milky Way plane are significantly higher than expected from extrapolations of the GeV data (the “TeV excess”). In Ref. Linden and Buckman (2018), it was shown that TeV halos could provide an explanation of this long-standing mystery. Our results also show that unresolved TeV halos could significantly contribute to the diffuse TeV $\gamma$-ray flux. We note that the diffuse emission is dominated by bright sources (Fig. 3). The predicted contribution for $\langle P_{0}\rangle$ = 300 ms is smaller than that estimated by Ref. Linden and Buckman (2018), which adopted a harder electron spectrum of $\alpha$ = 1.7 and $E_{\rm cut}$ = 100 TeV. In other words, a better determination of the average electron injection spectrum could increase the predicted flux from unresolved TeV halos, producing tighter constraints on $\langle P_{0}\rangle$ and $T_{\rm min}$. Interestingly, in the case of $\langle P_{0}\rangle$ = 50 ms or 120 ms, unresolved TeV halos can explain a significant fraction of the Milagro diffuse data without changing the spectral shape from that of our best-fit Geminga model.

IV.3 Summary of Allowed Models and Uncertainties

Our results are primarily affected by the $P_{0}$ distribution, and the 2HWC source count allows us to constrain 50 ms $\lesssim\langle P_{0}\rangle\lesssim$ 300 ms. Both the 50-ms and 300-ms models are barely allowed, and further investigations of TeV halo candidates will place stronger constraints on $\langle P_{0}\rangle$. This indicates that TeV halo observations can provide an important new probe of the $P_{0}$ distribution, which is difficult to constrain by pulsar statistics.

2HWC data require $T_{\rm min}\gtrsim 50$ kyr for $\langle P_{0}\rangle$ = 50 ms and $T_{\rm min}\lesssim 70$ kyr for $\langle P_{0}\rangle$ = 300 ms. The value of $T_{\rm min}$ is not well constrained for $\langle P_{0}\rangle\gtrsim$120 ms, but the firm identification of TeV halos around Geminga and Monogem suggests $T_{\rm min}\lesssim$100–300 kyr. Further observations are needed to better constrain this parameter. We stress that 50 kyr is not a strict minimum age for a TeV halo. Rather, we found that, operating under the assumption that the initial spin period of pulsars has an average $\langle P_{0}\rangle$ = 50 ms, this minimum age was required to ensure that the TeV halo number was consistent with data. However, the true initial period may have a larger mean, which would eliminate the need for such a cutoff. Alternatively, there may be significant variations between individual objects that are not taken into account in this model. Furthermore, our calculations have several uncertainties noted below, which could relieve the constraint on $T_{\rm min}$ for the $\langle P_{0}\rangle$ = 50 ms model.

We have fixed the distribution of $B_{0}$ to follow a lognormal distribution with $\langle\log_{10}B_{0}\rangle$ = 12.65 and $\sigma_{\log_{10}B_{0}}$ = 0.55. Other studies that examined the magnetic field evolution of pulsars find best-fit mean values that are about 2–4 times larger Popov et al. (2010); Gullón et al. (2014). In these models, the larger magnetic field causes pulsars to spin down faster, producing a smaller spin-down power for pulsars with ages exceeding $\sim$1 kyr. Adopting an alternative model with $\langle\log_{10}B_{0}\rangle$ = 13.10 and $\sigma_{\log_{10}B_{0}}$ = 0.65 as derived in Ref. Gullón et al. (2014), we find the predicted number of detectable TeV halos are reduced by a factor of $\sim$ 2. This increases the tension between $\langle P_{0}\rangle$ = 300 ms and current HAWC observations, but relieves some tension between $\langle P_{0}\rangle$ = 50 ms and the HAWC data. In particular, for these stronger magnetic fields, $\langle P_{0}\rangle$ = 50 ms models become consistent with HAWC upper limits for much smaller values of $T_{\rm min}\gtrsim$10 kyr. Further examinations of the $B_{0}$ distribution will also be important for the study of TeV halo populations.

We also note that throughout this section we focus on “Geminga-like” TeV halos. We can also adopt different models to take into account deviations from this assumption. We first study the effect of variations in the $\gamma$-ray efficiencies. There might be nearly an order of magnitude variation among individual sources, as noted in Sec. II. We examine alternative models where $\gamma$-ray fluxes are multiplied by a factor of $10^{x}$. If $x$ is fixed to 0.5 (i.e., all TeV halos are about 3 times brighter than the Geminga halo), then the number of detectable sources is increased by a factor of $\sim$3. Similarly, if $x$ is fixed to -0.5, the source count decreases by a factor of $\sim$3. We then examine the case where $x$ is a random variable drawn from a normal distribution with mean 0 and standard deviation 0.5. The primary effect of such a dispersion would be to smooth out the falling number-count distribution (Fig 3), and increase the number of detectable sources. We find that the number of detectable sources is increased by a factor of $\sim$2 in the case of $\langle P_{0}\rangle$ = 120 ms.

Finally, we study the possibility that pulsars younger than $T_{\rm min}$ produce TeV halos with different properties. In particular, at early ages, TeV halos may have smaller $\gamma$-ray efficiencies, because most of the injected energy should be lost to synchrotron emission and there may be less particle energy escaping into the TeV halos, as discussed in Sec. II. Throughout this paper, this effect is simply treated by sharply cutting off contributions from pulsars younger than $T_{\rm min}$, but one could alternatively assume a smooth changes in the $\gamma$-ray efficiencies. This could lead to a detectable population that does not exceed 2HWC constraints. To be more quantitative on this point, we examine alternative models where the $\gamma$-ray fluxes are smoothly reduced by a factor of $(T_{\rm age}/{\rm 340\ kyr})^{\beta}$ for pulsars younger than Geminga. This replaces the sharp cutoff ($T_{\rm min}$) in our standard formalism. We find that the $\langle P_{0}\rangle$ = 50 ms model does not produce too many TeV halos for $\beta\gtrsim$0.7. Further studies are needed to more thoroughly examine this parameter space.

In Table 2, we show each major uncertainty, mention an alternative model, and roughly indicate the net effect of this model on the predicted number of TeV halo sources. In addition to models mentioned above, we further test several different scenarios, which are explained in Appendix B. The exact effect of different uncertainties depends on the standard model that we use for comparison, so we adopt a constant default model of $\langle P_{0}\rangle$ = 120 ms and $T_{\rm min}$ = 10 kyr in all cases, and show their age dependence in Fig. 12 in Appendix B.

All of the uncertainties noted above could change the number of detectable sources by a factor of $\sim$2. While these changes are important, they are subdominant to the effect of variations in the $P_{0}$ distribution, and support our assertion that the $P_{0}$ distribution dominates the uncertainty in our models. Note that different $B_{0}$ distributions may lead to smaller number counts, while source variations may increase the number of detectable systems, implying that our default case occupies a reasonable middle value. More TeV halo observations would allow us to better examine these models, and place stronger constraints on pulsar properties.

V Future Directions

Upcoming surveys have great power to detect TeV halos. In Sec. V.1, we quantitatively assess the prospects for Galactic source searches with HAWC and CTA. In Sec. V.2, we do the same for extragalactic searches with CTA. In Sec. V.3, we show that detailed morphological studies of existing HESS sources could potentially identify many TeV halos.

V.1 Extended Source Survey with HAWC and CTA

We begin by outlining methods to identify TeV halos. One straightforward way to claim that TeV emission is powered by a pulsar is to detect the radio beam from the pulsed emission or to find a compact PWN at the center of the $\gamma$-ray source. We may also detect a TeV halo component in a composite (TeV halo + PWN) system by examining if its $\gamma$-ray emission can be fit by two morphological components rather than one. In the case of bow-shock pulsars, we can more clearly discriminate TeV halos from PWN, whose size is clearly determined by the stand-off radius Bykov et al. (2017).

In some cases we may be able to detect extended emission around PWNe in other wavelengths, from synchrotron radiation produced by the same electrons and positrons that escape the compact PWN and produce TeV halo emission. Interestingly, Refs. Uchiyama et al. (2009); Bamba et al. (2010) potentially detected such emission in X-rays, suggesting the potential for identifying TeV halos in multi-wavelength observations.

In Fig. 7, we show expectations for the TeV halo population that could be uncovered by HAWC observations. We assume a 10-yr sensitivity that is improved by a factor of $\sqrt{5}$ compared to the quoted sensitivity of the 2HWC catalog, following the same declination dependence. This corresponds to a sensitivity that is approximately 4% of the Geminga flux for sources residing in optimal sky positions. These predictions are pessimistic, because HAWC has recently installed an upgrade and increased the instrumented area Joshi et al. (2017), an effect which is not included in our calculation.

The 10-yr HAWC survey promises to discover a significant TeV halo population. Even in the pessimistic case of $\langle P_{0}\rangle$ = 300 ms, we expect that $\sim$20 sources (including 4 already found) will be detected for $T_{\rm min}=0$. In the most optimistic case (e.g., $\langle P_{0}\rangle$ = 50 ms and $T_{\rm min}$ = 50 kyr), HAWC would be expected to detect $\sim$80 TeV halos. We note that these cases are nearly ruled out by existing TeV halo observations (Sec. IV.1). Our intermediate case ($\langle P_{0}\rangle$ = 120 ms) predicts $\sim$70 sources if $T_{\rm min}=0$, though this model prediction is only barely allowed from the 2HWC source count (see Fig 4). Such a large number of sources would allow us to significantly improve our constraints on the spectral, morphological, and evolutionary properties of TeV halos.

We stress that our predictions are based on the Geminga-like assumption provided in Eq. (3), combined with standard models of pulsar population synthesis established by previous studies (e.g., Faucher-Giguère and Kaspi (2006); Bates et al. (2014)). If HAWC detects a significantly smaller number of TeV halos, it would indicate that Eq. (3) cannot be applied to all pulsars, and that the observed Geminga-like halos must have unusually large TeV $\gamma$-ray efficiencies or may have specific properties or environment that generate halos. Conversely, if significantly more sources are observed than predicted, it would indicate that observed sources have relatively dim halos, compared to the average population. Either result would substantially enhance our understanding of these systems.

In Fig. 8, we show our prediction for the Galactic longitude distribution of detected TeV halos (top) and the diffuse flux from unresolved TeV halos (bottom) for our intermediate case of $\langle P_{0}\rangle$ = 120 ms. In addition to 10-yr HAWC observations, we make a prediction for hypothetical HAWC-like telescope that uniformly observes the sky with a sensitivity that is 3% of the Geminga flux. We also show the predicted contribution from TeV halos to the Milagro diffuse measurement, which probed the region within $30^{\circ}<l<110^{\circ}$ and $136^{\circ}<l<216^{\circ}$ Atkins et al. (2005); Abdo et al. (2008), and also to the HAWC diffuse measurement, for the region that falls fairly within the field of view. The source count (top panel) shows the sizeable impact of a HAWC-like water Cherenkov telescope at the Southern hemisphere, like the Southern Gamma-Ray Survey Observatory Mostafa and HAWC Collaboration (2017); Schoorlemmer et al. (2017); Mostafa et al. (2017); it would allow us to probe a region at the edge or outside the HAWC field of view, where we expect a significant number of detectable TeV halos. It also demonstrates that the effect of source confusion is not large; we expect at most $\sim$4 sources over $\Delta l=5^{\circ}$, which means that the typical intersource spacing is large compared to the angular resolution ($\sim 0.1^{\circ}$ in radius) and the typical source size (also $\sim 0.1^{\circ}$ with significant variations).

In Fig. 8 (bottom), we show that HAWC measurements will greatly improve our understanding of the diffuse TeV flux. In particular, in the region of $40^{\circ}<l<100^{\circ}$, where the diffuse TeV excess was first identified, our predictions indicate that more than half of the diffuse emission from TeV halos will be resolved into individual sources by future HAWC surveys. Apart from shedding light on the nature of the diffuse TeV emission, this would also put constraints on the population of unresolved sources that contribute to the remaining diffuse flux.

We also make a prediction for the future Galactic Plane survey with CTA Acharya et al. (2017). We assume that CTA will observe from $l=0^{\circ}$ to $360^{\circ}$ and $|b|<3^{\circ}$ with a sensitivity of 3 mCrab, where 1 Crab is defined as a $\gamma$-ray flux above 1 TeV of 2.26$\times 10^{-11}$ cm*-2s-1*, and 3 mCrab corresponds to 2$\%$ of the Geminga flux (defined as $F_{\rm TeV}^{\mathcal{G}}$ in Sec. III). We assume a PSF size of $\theta_{\rm PSF}=0.05^{\circ}$ to take take into account the degradation of the sensitivity for extended sources. Because the PSF of CTA is smaller, this effect is more important compared to HAWC observations.

Even in the pessimistic case of $\langle P_{0}\rangle$ = 300 ms, we predict that $\sim$30 TeV halos could be detected. In the case of $\langle P_{0}\rangle$ = 50 ms with $T_{\rm min}$ = 50 kyr, we predict that about 160 TeV halos could be detected. Our intermediate model, $\langle P_{0}\rangle$ = 120 ms, also predicts $\sim$150 for $T_{\rm min}$ = 0. These detections will be highly complementary to HAWC observations. HAWC is located in the Northern hemisphere, while CTA is expected to have better sensitivity in the Southern hemisphere. Moreover, while HAWC is suited for spatially extended sources, CTA can find more distant and dimmer sources. In our prediction for HAWC, the 10–90% containment fraction of TeV halo distances corresponds to roughly 1–10 kpc. In contrast, for CTA, the 10–90% containment window spans from 3–15 kpc. Together, these observations can map out much of the Galaxy in TeV halos. This also indicates that another water Cherenkov telescope at the Southern hemisphere would be critical for detecting nearby sources throughout the Galactic Plane.

Finally, we study the effect of source confusion for CTA. First, we find that CTA is expected to detect $\sim 7$ sources in a 5∘ wide bin even in the densest regions (near the Galactic Center). Second, 90% of CTA sources has distances of above 3 kpc, which translates to an angular extension less than 0.5∘ assuming the source size of 25 pc. These indicate that the source confusion has marginal effect on our results, but further understanding of luminosity function of other sources and CTA properties are needed to better quantify this effect.

V.2 Extragalactic Survey with CTA

Milky Way TeV halo searches must deal with large angular source size and distance uncertainties. One way to avoid these issues is to search for extragalactic sources. A good target is the Large Magellanic Cloud (LMC), which is nearby and face-on, and which will also be extensively observed by the CTA as part of its Key Scientific Program Acharya et al. (2017). CTA observations are expected to achieve an integrated energy flux sensitivity of 3 $\times$ 10*-14* erg cm*-2* s*-1* above 1 TeV. The collaboration expects to uncover $\sim$10 sources that are primarily SNRs, without including TeV halos.

We estimate the number of TeV halos that can be detected in the LMC by the CTA survey. We adopt a standard distance of 50 kpc, and count the number of halos with luminosities exceeding 9 $\times~{}10^{33}$ erg s*-1* above 1 TeV. The birth rate of pulsars in the LMC is normalized to be 0.005 yr*-1*, which is the lower value obtained by a previous study of LMC pulsar population modelling Ridley and Lorimer (2010). Though the interstellar infrared radiation field in the LMC is weaker compared to the Milky Way Israel et al. (2010), the predicted TeV halo flux is reduced only by a factor of 1.3 even if we set $\rho_{\rm IR}$ = 0, due to the contribution from the cosmic microwave background photons. Since this modification is degenerate with a number of uncertainties, we do not take this into account in what follows.

In Fig. 9, we show that CTA will likely detect at least $\sim$1, and potentially $\sim$30, extragalactic TeV halos in the LMC, substantially increasing the total number of sources detected with this survey. These observations will provide more important information than the source count alone, shedding insight into the brightest TeV halos in a region without significant distance uncertainties. Thus, CTA observations will provide complementary constraints to HAWC Milky Way observations. The differentiation of extragalactic TeV halos from PWNe will be challenging, because both will appear pointlike even with the unparalleled angular resolution of CTA. The best path forward will be to employ followup radio observations of synchrotron counterparts to further examine the emission morphology. Because the size of halos is poorly understood, there may also be a possibility that we could observe TeV halos extended beyond a radius of $\sim$ 100 pc, which could be detected as extended sources even in the LMC.

We note that the observed radio pulsars in the LMC appear to have a strikingly different distribution of spindown powers compared to expectations from pulsar evolution (Fig. 3). In particular, current observations detect two very bright ($\dot{E}>10^{38}$ erg s*-1*) pulsars, while the 11 other detected pulsars have low $\dot{E}$ ($\dot{E}<6\times 10^{34}$ erg s*-1*). There are no pulsars in between these ranges Manchester et al. (2005). This is most likely due to the combination of selection effects, weak correlations between $\dot{E}$ and radio luminosity Szary et al. (2014), and the randomness of pulse radiation beamed toward us. Future pulsar surveys may enable us to better examine the pulsar population in the LMC. Since TeV halo emissions are expected to be more isotropic and may better correlate with $\dot{E}$ than radio pulse emissions, they could provide complementary information regarding the population of bright pulsars, potentially resolving this tension, or confirming it. In the latter case, it would demand significant modifications to the theory of pulsar formation and evolution.

Finally, we note that we might also be able to observe a similar number of TeV halos in the Small Magellanic Cloud, because it has a distance and pulsar formation rate comparable to the LMC Ridley and Lorimer (2010). This would potentially provide information regarding the evolution of TeV halos in low metallicity environments.

V.3 Followup Study for HESS Sources

So far, we have focused on the existing survey catalog by HAWC, which is suited for extended source surveys. However, existing source catalogs from imaging air Cherenkov telescopes should also contain as-yet identified TeV halos.

Here we focus on the HESS Galactic Plane Survey (HGPS) catalog, which has detected 78 sources in total, 42 of which are associated with ATNF pulsars Abdalla et al. (2018a). Five of these tentative pulsar associations are known to be either an SNR or a binary, as well as the Arc and Galactic Center, while the remaining 37 sources are categorized as firmly identified PWN (including composite system), candidate PWN, or unidentified sources. We examine how many of these sources could be interpreted as TeV halos.

We make a prediction following the methodology for the CTA Galactic Plane survey in Sec. V.1. We include sources between Galactic longitudes spanning from $l=250^{\circ}$ to $65^{\circ}$ and latitude $|b|<3^{\circ}$. The sensitivity of HGPS is non-uniform across the Galactic Plane. Since our goal is not to make precise estimates of the HESS sensitivity, we simply assume a sensitivity of 1$\%$ Crab for point sources and utilize a PSF of $\theta_{\rm PSF}=0.08^{\circ}$. We further assume that any source of $\theta_{\rm size}>0.7^{\circ}$ is not observed, because HGPS is not able to detect such an extended source. This removes the contribution from any halos within about 700 pc of Earth.

In Fig. 10 (top), we show that 10–50 sources in the HGPS catalog could be TeV halos. This indicates that detailed morphological studies of HGPS sources could uncover many TeV halos in this catalog. The definitive identification of these sources would be important in constraining particle transport due to the unparalleled angular resolution of HESS.

In Fig. 10 (bottom), we compare the predicted number of TeV halos to the observed number of radio pulsar associations. Our model predicts that, among the 37 sources associated with ATNF pulsars, $\sim$6–20 sources could be TeV halos. Interestingly, the shape of our predicted age distribution matches observed data well. Because our predictions are based on the assumption that the $\gamma$-ray flux is proportional to $\dot{E}$, this agreement suggests that these HGPS sources are powered by pulsar activity, either PWNe or TeV halos, rather than SNR.

VI Conclusions

TeV halos are a new class of $\gamma$-ray sources Linden et al. (2017); Abdo et al. (2009); Abeysekara et al. (2017a); Aharonian et al. (1997); Yüksel et al. (2009); Linden and Buckman (2018); Hooper et al. (2018); Evoli et al. (2018); Hooper and Linden (2018a); Riviere et al. (2017); Brisbois et al. (2018). They are bright, have hard spectra, and are spatially extended. They are powered by electrons and positrons that escaped from the PWNe, but which remain confined in a region where diffusion is strongly suppressed. Empirical arguments suggest TeV halos are common around pulsars. However, among many candidate sources, only four TeV halos have been confirmed so far. The rest are likely undetected due to the diminishing sensitivity of TeV instruments to extended $\gamma$-ray sources.

In this work, for the first time, we theoretically quantify the role of future surveys to detect more TeV halos. We also study new implications for pulsar physics and existing $\gamma$-ray sources. We use standard methods for pulsar population synthesis and focus on a model where the TeV halo luminosity is calculated based on Geminga observations. Our analysis produced three main results.

• •

TeV halos could be the most important source class in future TeV $\gamma$-ray surveys. Within the context of our standard Geminga-like model, and utilizing the range of $P_{0}$ and $T_{\rm min}$ that are consistent with current datasets, we predict that HAWC will eventually detect $\sim$20–80 TeV halos, and future Galactic surveys by CTA will also find $\sim$30–160 halos, for a total of $\sim$50–240 halos. Further, CTA can potentially detect $\sim$10 TeV halos in the LMC and SMC. This indicates that TeV halos could be the dominant source class in the TeV $\gamma$-ray sky. Such a large number of sources would allow us to examine their properties and evolution in great detail.

• •

Further studies of unidentified TeV sources and “TeV PWN” are needed. We find that the HGPS catalog might contain $\sim$10–50 TeV halos, which are currently classified as either unidentified sources or PWNe. These results have three implications. First, imaging air Cherenkov telescopes like HESS, MAGIC and VERITAS can also play an important role in studying TeV halos. In particular, their high angular resolution would be critical in examining particle transport inside the halo. Second, it might be important for modelling of “TeV PWN” to take into account the emission from TeV halo regions. Third, X-ray and radio observations of “TeV PWN” may find many extended halos around compact PWNe. This synchrotron emission counterpart could help to identify TeV halos.

• •

TeV halos observations can constrain pulsar properties. Our predictions are primarily affected by the distribution of the initial spindown period, which is not well constrained by pulsar population studies. In other words, TeV halo observations can provide complementary constraints to existing radio surveys. Current 2HWC data allow 50 ms $\lesssim$ $\langle P_{0}\rangle$ $\lesssim$ 300 ms, and further studies will place tighter constraints.

We finally note that TeV halo observations may unlock new opportunities to study astrophysics.

• •

TeV halo observations would allow us to detect pulsars with radio emission not aligned toward Earth and hence which have been missed in previous blind searches Linden et al. (2017). Further, the angular size of halos could provide useful distance estimations for Galactic pulsars. Thus, many observations of TeV halos could allow us to map out pulsars in the Galaxy, including misaligned systems. This new method would work as an independent and complementary method compared to radio observations. In this regard, the Southern Gamma-Ray Survey Observatory would play an important role in detecting TeV halos across the Galactic Plane, especially in the inner Galaxy. In addition, next-generation telescopes like LHASSO Di Sciascio and LHAASO Collaboration (2016) can find more TeV halos.

• •

TeV halo observations would substantially improve our understanding of total galactic $\gamma$-ray emission. Usually, galactic $\gamma$-ray emission is assumed to be dominated by hadronic processes induced by diffusing protons and nuclei, especially in the GeV energy range. The bright and hard-spectrum emission from TeV halos suggests that leptonic emission mechanisms may be important for the diffuse emission in the TeV energy range. This has two implications for the cosmic background emission. First, the hard spectrum of TeV halo emission might make ordinary galaxies more important for TeV $\gamma$-ray background than expected only from hadronic emission, which falls off steeply. Second, the leptonic nature of TeV halo emission may make star-forming galaxies less important for the TeV neutrino background than expected from the assumption that all TeV galactic $\gamma$-ray emission is hadronic. In particular, predicting neutrino flux from galaxies by simply extrapolating their $\gamma$-ray flux could result in a substantial overestimate.

• •

TeV halos could help pinpoint the sources of IceCube neutrinos in our Galaxy. A promising way to search for neutrino emitters is to look into $\gamma$-ray source catalogs. However, if most $\gamma$-ray sources are TeV halos, there is less room for hadronic sources. This has both positive and negative implications for neutrino astronomy. It is unfortunate, since we can only expect high-energy neutrinos from a small fraction of identifiable $\gamma$-ray sources. On the other hand, if we can identify TeV halos in $\gamma$-ray source catalogs, we can ignore their neutrino contributions and reduce the trials factor in IceCube neutrino cross-correlations.

• •

Existing TeV halo observations indicate that pulsars contribute to the cosmic-ray electron and positron flux. However, future observations are necessary to understand the exact degree they contribute, which classes of systems are important, and the constraints which can be put on residual contributions by exotic physics such as dark matter annihilation. In particular, follow up observations of TeV halos with GeV $\gamma$-rays and other wavelengths can provide important complementary information capable of constraining the pulsar contribution to the positron excess, which is seen in GeV energy range.

Finally, it is remarkable that such an important source class escaped identification until the development of TeV $\gamma$-ray instruments, especially those that can detect extended sources. Although TeV halos may have corresponding emission in radio, X-ray, and GeV photons, this was not sufficiently obvious to recognize this source class. The significance of these observations thus predicts the importance of future TeV surveys in understanding the multi-wavelength sky.

Acknowledgements

We are grateful to Felix Aharonian, Aya Bamba, Chad Brisbois, Mattia Di Mauro, Fiorenza Donato, Henrike Fleischhack, Dan Hooper, Petra Huentemeyer, Wataru Ishizaki, Kelly Malone, Miguel Mostafa, Kohta Murase, Pat Slane, Andrew Smith, Terry Walker, and especially Katie Auchettl for discussions and helpful comments on the manuscript. We are grateful to Sophie Beacom for help with Fig. 1. We also thank two anonymous referees for helpful comments. TS is supported by Research Fellowship of Japan Society for the Promotion of Science (JSPS), and also supported by JSPS KAKENHI Grant Number JP 18J20943. JFB is supported by NSF grant PHY-1714479.

Appendix A Effect of Uniform Pulsar Period Distributions

One of the most important modeling assumptions in this paper is the spin-period of the pulsar at birth. In the main paper, we utilized several models that utilized a Gaussian distribution to describe the initial pulsar period distribution. We chose three average periods of 50 ms, 120 ms, or 300 ms with this Gaussian, along with variances $\sigma_{P_{0}}$ = 50/$\sqrt{2}$ ms, $\sigma_{P_{0}}$ = 60 ms and $\sigma_{P_{0}}$ = 150 ms, respectively. However, recent modelling of pulsar population indicates that a uniform period distribution of $0<P_{0}<500$ ms may be consistent with the pulsar statistics Gullón et al. (2015).

In Fig. 11, we show a version of Fig. 4 that shows the expected number of TeV halos in the 2HWC catalog as a function of the pulsar age for models with uniform initial spin-period distributions between the breakup limit and $2\langle P_{0}\rangle$. In addition to three cases of $\langle P_{0}\rangle$ =50, 120, and 300 ms, we show the case of a uniform $P_{0}$ distribution up to 500 ms, as suggested in Ref. Gullón et al. (2015).

We find that, in general, models with a uniform pulsar distribution of $P_{0}$ predict a greater number of pulsar than Gaussian models. This is because the uniform distribution produces a larger number of sources with very small $P_{0}$. In general, the effect of these models on the number of pulsars as a function of $T_{\rm age}$ is about a factor of $\sim$ 2, compared to models with a Gaussian distribution centered at $\langle P_{0}\rangle$.

Appendix B Model Uncertainties

Here, we describe all uncertainties in Table 2. In Fig. 12, we show how our results are changed for alternative models for the case of $\langle P_{0}\rangle$ = 120 ms.

• •

$P_{0}$** distribution**: In the main text, we adopt a Gaussian distribution following Refs. Faucher-Giguère and Kaspi (2006); Watters and Romani (2011). In Fig. 12 (top left), we compare this with a uniform distribution, motivated by Ref. Gullón et al. (2015).

• •

$B_{0}$** distribution**: In the main text, we adopt a log-normal distribution with mean $\langle\log_{10}B_{0}\rangle$ = 12.65 and $\sigma_{\log_{10}B_{0}}$ = 0.55 following Ref. Faucher-Giguère and Kaspi (2006). In Fig. 12 (top left), we compare this with another lognormal distribution with $\langle\log_{10}B_{0}\rangle$ = 13.10 and $\sigma_{\log_{10}B_{0}}$ = 0.65, as defined in Ref. Gullón et al. (2014).

• •

$\dot{E}$** dependence**: In the main text, we adopt $L_{\gamma}\propto\dot{E}$. In Fig. 12 (top right), we compare this with different models which adopt $L_{\gamma}\propto\dot{E}^{0.8}$ and $L_{\gamma}\propto\dot{E}^{1.2}$, normalizing the $\gamma$-ray flux with the Geminga halo.

• •

age dependence: In the main text, we assume that $\gamma$-ray efficiency $L_{\gamma}/\dot{E}$ is constant for all pulsars. In Fig. 12 (bottom left), we compare this with different models where $L_{\gamma}/\dot{E}$ depends on the age of pulsars as $L_{\gamma}/\dot{E}\propto(T_{\rm age})^{0.5}$ and $L_{\gamma}/\dot{E}\propto(T_{\rm age})^{-0.5}$, normalizing the $\gamma$-ray flux with the Geminga halo and adopting 340 kyr as the age of Geminga. In the case of $L_{\gamma}/\dot{E}\propto(T_{\rm age})^{0.5}$, we do not assume age dependence for sources older than 340 kyr, in order to avoid producing unrealistically bright late-time sources.

• •

source-to-source scatter: In the main text, we assume that $L_{\gamma}/\dot{E}$ is constant for all systems. In Fig. 12 (bottom right), we compare this with different models where $\log_{10}(L_{\gamma}/\dot{E})$ is a random variable drawn from a normal distribution with mean 1 and standard deviation 0.5.