Evidences of low-diffusion bubbles around Galactic pulsars

TL;DR

This paper investigates the low-diffusion regions around Galactic pulsars, analyzing gamma-ray observations to estimate diffusion coefficients and assess the detectability of inverse Compton scattering halos with current and future gamma-ray telescopes.

Contribution

It provides the first detailed analysis of the morphology of ICS gamma-ray halos around pulsars, estimating low diffusion coefficients and predicting detection prospects with CTA, HAWC, and HESS.

Findings

Galactic pulsars are surrounded by low-diffusion bubbles of at least 80 pc.

The diffusion coefficient around these pulsars is estimated to be 2-30 x 10^{26} cm^2/s at 1 TeV.

Future gamma-ray experiments can significantly increase the number of detected ICS halos.

Abstract

Recently, a few-degrees extended $\gamma$-ray halo in the direction of Geminga pulsar has been detected by HAWC, Milagro and Fermi-LAT. These observations can be interpreted with positrons ($e^+$) and electrons ($e^-$) accelerated by Geminga pulsar wind nebula (PWN), released in a Galactic environment with a low diffusion coefficient ($D_0$), and inverse Compton scattering (ICS) with the interstellar radiation fields. We inspect here how the morphology of the ICS $\gamma$-ray flux depends on the energy, the pulsar age and distance, and the strength and extension of the low-diffusion bubble. In particular we show that $\gamma$-ray experiments with a peak of sensitivity at TeV energies are the most promising ones to detect ICS halos. We perform a study of the sensitivity of HAWC, HESS and the future CTA experiment finding that, with efficiencies of the order of a few %, the first two experiments should have already detected a few tens of ICS halos while the latter will increase the number of detections by a factor of 4. We then consider a sample of sources associated to PWNe and detected in the HESS Galactic plane survey and in the second HAWC catalog. We use the information available in these catalogs for the $\gamma$-ray spatial morphology and flux of these sources to inspect the value of $D_0$ around them and the $e^{\pm}$ injection spectrum. All sources are detected as extended with a $\gamma$-ray emission extended about $15-80$ pc. Assuming that most of the $e^{\pm}$ accelerated by these sources have been released in the interstellar medium, the diffusion coefficient is $2-30 \cdot 10^{26}$ cm$^2$/s at 1 TeV, i.e. two orders of magnitude smaller than the value considered to be the average in the Galaxy. These observations imply that Galactic PWNe have low-diffusion bubbles with a size of at least 80 pc.

Tables

Table 1. Table 1: List of the pulsars from the ATNF catalog in the HAWC field of view with the brightest predicted ICS halo flux at 10 TeV. We list the pulsar name, Galactic coordinates, distance, age and spin-down luminosity taken from the ATNF catalog. Then, we report the predicted extension θ 68 subscript 𝜃 68 \theta_{68} and ICS flux Φ γ 10 ​ TeV subscript superscript Φ 10 TeV 𝛾 \Phi^{10\,\rm{TeV}}_{\gamma} both calculated at 10 TeV and assuming D 0 = 7 ⋅ 10 25 subscript 𝐷 0 ⋅ 7 superscript 10 25 D_{0}=7\cdot 10^{25} cm 2 /s. Finally, we display the name as in 2HWC catalog and the classification given in TeVCat. Sources labeled as UNID are unidentified in the TeVCat catalog but are associated with potential ICS halo in our analysis since they have a pulsar within a small angular distance.

PSR	$l$	$b$	$d$	$T$	$\dot{E}$	${\Phi }_{\gamma }^{10\mathrm{TeV}}$	${\theta }_{68}$	Name	Class
	[deg]	[deg]	[kpc]	[kyr]	[erg/s]	[${\left(\mathrm{TeV}{\mathrm{cm}}^{2}s\right)}^{-1}$]	[deg]
J1826-1256	18.56	-0.38	1.55	14	$3.6\cdot {10}^{36}$	$2.5\cdot {10}^{-13}$	0.89	2HWC J1825-134	UNID
J2021+3651	75.22	0.11	1.80	17	$3.4\cdot {10}^{36}$	$1.6\cdot {10}^{-13}$	0.82	2HWC J2019+367	UNID
J1813-1246	17.24	2.44	2.63	43	$6.2\cdot {10}^{36}$	$8.6\cdot {10}^{-14}$	0.60	2HWC J1812-126	UNID
J1907+0602	40.18	-0.89	2.37	20	$2.8\cdot {10}^{36}$	$6.7\cdot {10}^{-14}$	0.64	2HWC J1908+063	UNID
J0633+1746	195.13	4.27	0.19	342	$3.3\cdot {10}^{34}$	$5.8\cdot {10}^{-14}$	6.54	GEMINGA PWN	TEV HALO
B0656+14	201.11	8.26	0.29	111	$3.8\cdot {10}^{34}$	$3.4\cdot {10}^{-14}$	4.71	2HWC J0700+143	TEV HALO
B1951+32	68.77	2.82	3.00	107	$3.7\cdot {10}^{36}$	$3.0\cdot {10}^{-14}$	0.46	undetected	undetected
J1811-1925	11.18	-0.35	5.00	23	$6.4\cdot {10}^{36}$	$2.8\cdot {10}^{-14}$	0.30	2HWC J1809-190	UNID
B1823-13	18.00	-0.69	3.61	21	$2.8\cdot {10}^{36}$	$2.6\cdot {10}^{-14}$	0.41	2HWC J1825-134	UNID
J1935+2025	56.05	-0.05	4.60	21	$4.7\cdot {10}^{36}$	$2.5\cdot {10}^{-14}$	0.32	SNR G054.1+00.3	PWN
J1954+2836	65.24	0.38	1.96	69	$1.1\cdot {10}^{36}$	$2.3\cdot {10}^{-14}$	0.77	2HWC J1955+285	UNID
J1809-1917	11.09	0.08	3.27	51	$1.8\cdot {10}^{36}$	$1.5\cdot {10}^{-14}$	0.47	2HWC J1809-190	UNID
J1838-0655	25.25	-0.20	6.60	23	$5.6\cdot {10}^{36}$	$1.3\cdot {10}^{-14}$	0.22	2HWC J1837-065	PWN
J1856+0245	36.01	0.06	6.32	21	$4.6\cdot {10}^{36}$	$1.2\cdot {10}^{-14}$	0.23	2HWC J1857+027	UNID
J1958+2846	65.88	-0.35	1.95	22	$3.4\cdot {10}^{35}$	$1.2\cdot {10}^{-14}$	0.79	2HWC J1955+285	UNID
J1740+1000	34.01	20.27	1.23	114	$2.3\cdot {10}^{35}$	$1.1\cdot {10}^{-14}$	1.15	undetected	undetected
J1913+1011	44.48	-0.17	4.61	169	$2.9\cdot {10}^{36}$	$9.1\cdot {10}^{-15}$	0.27	2HWC J1912+099	SHELL
J1837-0604	25.96	0.27	4.77	34	$2.0\cdot {10}^{36}$	$8.6\cdot {10}^{-15}$	0.32	2HWC J1837-065	UNID
J1907+0631	40.52	-0.48	3.40	11	$5.3\cdot {10}^{35}$	$6.9\cdot {10}^{-15}$	0.41	2HWC J1908+063	UNID
J1928+1746	52.93	0.11	4.34	83	$1.6\cdot {10}^{36}$	$6.5\cdot {10}^{-15}$	0.30	2HWC J1928+177	UNID
J0633+0632	205.09	-0.93	1.35	59	$1.2\cdot {10}^{35}$	$5.8\cdot {10}^{-15}$	1.14	HAWC J0635+070	TEV HALO
J1831-0952	21.90	-0.13	3.68	128	$1.1\cdot {10}^{36}$	$5.6\cdot {10}^{-15}$	0.39	2HWC J1831-098	PWN
J1828-1101	20.50	0.04	4.77	77	$1.6\cdot {10}^{36}$	$5.3\cdot {10}^{-15}$	0.28	2HWC J1831-098	UNID

Table 2. Table 2: Same as in Tab.  1 but for source not detected by HAWC so far.

PSR	$l$	$b$	$d$	$T$	$\dot{E}$	${\Phi }_{\gamma }^{10\mathrm{TeV}}$	${\theta }_{68}$
	[deg]	[deg]	[kpc]	[kyr]	[erg/s]	[${\left(\mathrm{TeV}{\mathrm{cm}}^{2}s\right)}^{-1}$]	[deg]
B1951+32	68.77	2.82	3.00	107	$3.7\cdot {10}^{+36}$	$3.0\cdot {10}^{-14}$	0.46
J1740+1000	34.01	20.27	1.23	114	$2.3\cdot {10}^{+35}$	$1.1\cdot {10}^{-14}$	1.15
J1755-0903	18.32	8.15	0.23	3870	$4.4\cdot {10}^{+33}$	$5.0\cdot {10}^{-15}$	5.48
J0729-1448	230.39	1.42	2.68	35	$2.8\cdot {10}^{+35}$	$4.0\cdot {10}^{-15}$	0.60
J0631+1036	201.22	0.45	2.10	44	$1.7\cdot {10}^{+35}$	$3.8\cdot {10}^{-15}$	0.76
B1929+10	47.38	-3.88	0.31	3100	$3.9\cdot {10}^{+33}$	$2.4\cdot {10}^{-15}$	4.11
J0538+2817	179.72	-1.69	1.30	618	$4.9\cdot {10}^{+34}$	$1.8\cdot {10}^{-15}$	1.02
J2043+2740	70.61	-9.15	1.48	1200	$5.6\cdot {10}^{+34}$	$1.6\cdot {10}^{-15}$	0.88
J1846+0919	40.69	5.34	1.53	360	$3.4\cdot {10}^{+34}$	$9.2\cdot {10}^{-16}$	0.86
J1900-09	25.46	4.73	0.30	1500	$1.2\cdot {10}^{+33}$	$7.9\cdot {10}^{-16}$	4.26
J2055+2539	70.69	-12.52	0.62	1240	$4.9\cdot {10}^{+33}$	$7.8\cdot {10}^{-16}$	1.97
B1702-19	3.19	13.03	0.75	1140	$6.1\cdot {10}^{+33}$	$6.6\cdot {10}^{-16}$	1.75
J0611+1436	195.38	-2.00	0.89	1070	$8.0\cdot {10}^{+33}$	$6.1\cdot {10}^{-16}$	1.47
J0357+3205	162.76	-16.01	0.83	540	$5.9\cdot {10}^{+33}$	$5.3\cdot {10}^{-16}$	1.59
B1930+22	57.36	1.55	10.90	40	$7.5\cdot {10}^{+35}$	$5.3\cdot {10}^{-16}$	0.12
B0450-18	217.08	-34.09	0.40	1510	$1.4\cdot {10}^{+33}$	$5.2\cdot {10}^{-16}$	3.21
B0950+08	228.91	43.70	0.26	17500	$5.6\cdot {10}^{+32}$	$4.9\cdot {10}^{-16}$	4.88
J2006+3102	68.67	-0.53	6.03	104	$2.2\cdot {10}^{+35}$	$4.5\cdot {10}^{-16}$	0.21
B0919+06	225.42	36.39	1.10	497	$6.8\cdot {10}^{+33}$	$3.5\cdot {10}^{-16}$	1.22
B1706-16	5.78	13.66	0.56	1640	$8.9\cdot {10}^{+32}$	$1.7\cdot {10}^{-16}$	2.17
J1921+0812	43.71	-2.93	2.90	622	$2.3\cdot {10}^{+34}$	$1.7\cdot {10}^{-16}$	0.45
J1816-0755	21.87	4.09	3.13	532	$2.5\cdot {10}^{+34}$	$1.6\cdot {10}^{-16}$	0.41
B1821-19	12.28	-3.11	3.70	573	$3.0\cdot {10}^{+34}$	$1.4\cdot {10}^{-16}$	0.33
J1848+0647	38.70	3.65	1.13	916	$2.7\cdot {10}^{+33}$	$1.3\cdot {10}^{-16}$	1.18

Table 3. Table 3: List of the pulsars considered in our analysis. See the text for more information on the criteria we use to select them. We list the pulsar name, Galactic coordinates, pulsar distance, age and spin-down luminosity, association name, extension as given in HGSP or Abeysekara et al. ( 2017a ) ( θ gauss HESS subscript superscript 𝜃 HESS gauss \theta^{\rm{HESS}}_{\rm{gauss}} ). We also show the angular size ( θ gauss subscript 𝜃 gauss \theta_{\rm{gauss}} ) found by fitting with a gaussian function the source surface brightness derived with the HESS flux maps with R c = 0.1 ∘ subscript 𝑅 𝑐 superscript 0.1 R_{c}=0.1^{\circ} (see the text for further details). Finally, we report predicted size of the ICS halo at 1 TeV using D 0 = 7 ⋅ 10 25 subscript 𝐷 0 ⋅ 7 superscript 10 25 D_{0}=7\cdot 10^{25} cm 2 /s ( θ 68 subscript 𝜃 68 \theta_{68} ), and classification as in TeVcat. The first (second) block corresponds to o ​ l ​ d 𝑜 𝑙 𝑑 old ( y ​ o ​ u ​ n ​ g 𝑦 𝑜 𝑢 𝑛 𝑔 young ) sources.

PSR	$l$	$b$	$d$	$T$	$\dot{E}$	Name	${\theta }_{\mathrm{gauss}}^{\mathrm{HESS}}$	${\theta }_{\mathrm{gauss}}$	${\theta }_{68}$	Type
	[deg]	[deg]	[kpc]	[kyr]	[erg/s]		[deg]	[deg]	[deg]
J1016-5857	284.08	-1.88	3.16	21	2.6$\cdot {10}^{36}$	HESS J1018-589 B	$0.15\pm 0.03$	$0.14\pm 0.03$	0.14	PWN
J1028-5819	285.06	-0.50	1.42	90	8.3$\cdot {10}^{35}$	HESS J1026-582	$0.13\pm 0.04$	$0.18\pm 0.04$	0.27	PWN
J1459-6053	317.89	-1.79	1.84	65	9.1$\cdot {10}^{35}$	HESS J1458-608	$0.37\pm 0.03$	$0.37\pm 0.10$	0.44	PWN
J1632-4757	. 336.30	0.08	4.84	240	5.0$\cdot {10}^{34}$	HESS J1632-478	$0.18\pm 0.02$	$0.25\pm 0.04$	0.14	PWN
J1718-3825	348.95	-0.43	3.49	90	1.3$\cdot {10}^{36}$	HESS J1718-385	$0.12\pm 0.01$	$0.13\pm 0.02$	0.09	PWN
J1809-1917	11.18	-0.35	3.27	51.7	1.8$\cdot {10}^{36}$	HESS J1809-193(2HWC J1809-190)	$0.40\pm 0.05$	$0.35\pm 0.03$	0.37	UNID
J1813-1246	17.24	2.44	2.63	43	6.2$\cdot {10}^{36}$	HESS J1813-126(2HWC J1812-126)	$0.21\pm 0.03$	$0.20\pm 0.09$	0.33	UNID
B1823-13	18.00	-0.69	3.61	21	2.8$\cdot {10}^{36}$	HESS J1825-137(2HWC J1825-134)	$0.46\pm 0.03$	$0.36\pm 0.02$	0.28	HALO
J1831-952	21.90	-0.13	3.68	128	1.1$\cdot {10}^{36}$	HESS J1831-098(2HWC J1831-098)	$0.15$	$0.19\pm 0.05$	0.21	PWN
J1838-0655	25.25	-0.20	6.60	23	5.6$\cdot {10}^{36}$	HESS J1837-069(2HWC J1837-065)	$0.36\pm 0.03$	$0.31\pm 0.02$	0.14	PWN
J1841-0524	27.02	-0.33	4.12	30.2	1.0$\cdot {10}^{36}$	HESS J1841-055	$0.40\pm 0.03$	$0.50\pm 0.08$	0.21	UNID
J1856+0245	36.01	0.06	6.32	21	4.6$\cdot {10}^{36}$	HESS J1857+026(2HWC J1857+027)	$0.26\pm 0.06$	$0.23\pm 0.04$	0.14	UNID
J1857+0143	35.17	-0.57	4.57	71	4.5$\cdot {10}^{35}$	HESS J1858+020	$0.08\pm 0.02$	$0.12\pm 0.04$	0.21	UNID
J1907+0602	40.18	-0.89	2.37	20	2.8$\cdot {10}^{36}$	HESS J1908+063(2HWC J1908+063)	$0.49\pm 0.03$	$0.50\pm 0.10$	0.26	UNID
J1913+1011	44.48	-0.17	4.61	169	2.9$\cdot {10}^{36}$	HESS J1912+101(2HWC J1912+099)	$0.49\pm 0.04$	$0.42\pm 0.06$	0.16	SHELL
B0833-45	263.55	-2.79	0.28	11.3	6.9$\cdot {10}^{36}$	HESS J0835-455(Vela X)	$0.58\pm 0.05$	$0.48\pm 0.06$	0.66	PWN
J1301-6305	304.10	-0.24	10.72	11	1.7$\cdot {10}^{36}$	HESS J1303-631	$0.18\pm 0.02$	$0.20\pm 0.02$	0.10	PWN
J1357-6429	309.92	-2.51	3.10	7.3	3.1$\cdot {10}^{36}$	HESS J1356-645	$0.23\pm 0.02$	$0.22\pm 0.04$	0.28	PWN
J1420-6048	313.54	0.23	5.63	13	1.0$\cdot {10}^{37}$	HESS J1420-607	$0.08\pm 0.01$	$0.13\pm 0.02$	0.14	PWN
J1617-5055	332.50	-0.28	4.74	8.1	1.6$\cdot {10}^{36}$	HESS J1616-508	$0.23\pm 0.03$	$0.22\pm 0.02$	0.13	PWN
J1640-4631	338.32	-0.02	12.75	3.4	4.4$\cdot {10}^{36}$	HESS J1640-465	$0.18\pm 0.02$	$0.13\pm 0.01$	0.10	PWN
B1706-44	343.10	-2.69	2.60	18	3.4$\cdot {10}^{36}$	HESS J1708-443	$0.28\pm 0.03$	$0.24\pm 0.07$	0.23	PWN
J1813-1749	12.82	-0.02	4.70	5.6	5.6$\cdot {10}^{37}$	HESS J1813-178(2HWC J1814-173)	$0.049\pm 0.004$	$0.019\pm 0.003$	0.10	PWN
J1826-1256	18.56	-0.38	1.55	14	3.6$\cdot {10}^{36}$	HESS J1826-130(2HWC J1825-134)	$0.15\pm 0.02$	$0.23\pm 0.06$	0.28	UNID
J1833-1034	21.50	-0.89	4.10	4.9	3.4$\cdot {10}^{37}$	HESS J1833-105		$0.08\pm 0.02$	0.30	PWN
J0633+1746	195.13	4.27	0.19	342	3.3$\cdot {10}^{34}$	GEMINGA(2HWC J0635+180)	$5.5\pm 0.7$	$1.4\pm 0.2$	6.54	HALO
B0656+14	201.11	8.26	0.29	111	3.8$\cdot {10}^{34}$	MONOGEM(2HWC J0700+143)	$4.8\pm 0.6$	$2.7\pm 0.4$	4.71	HALO

Table 4. Table 4: Results of our analysis for D 0 subscript 𝐷 0 D_{0} and η 𝜂 \eta . We report the source name, the size of extension of the ICS halo found using the function in Eq.  18 ( θ ICS / 2 subscript 𝜃 ICS 2 \theta_{\rm{ICS}}/2 ), half of the the size of the ICS halo calculated using θ ICS / 2 subscript 𝜃 ICS 2 \theta_{\rm{ICS}}/2 , the best fit and 1 ​ σ 1 𝜎 1\sigma error for D 0 subscript 𝐷 0 D_{0} and the e ± superscript 𝑒 plus-or-minus e^{\pm} spectral index.

Name	${\theta }_{\mathrm{ICS}}/2$	Size	$D_0$	${\gamma }_e$
	[deg]	pc	${10}^{25}$ [cm2/s]
HESS J1018-589B	$0.27\pm 0.10$	$15\pm 5$	${2.2}_{-0.9}^{+1.5}$(${2.4}_{-0.9}^{+1.4}$)	$2.5\pm 0.1$
HESS J1026-582	$0.25\pm 0.09$	$6\pm 2$	${0.48}_{-0.20}^{+0.38}$(${0.31}_{-0.13}^{+0.23}$)	$1.6\pm 0.1$
HESS J1458-608	$1.20\pm 0.73$	$38\pm 23$	${4.7}_{-2.8}^{+12.3}$(${4.5}_{-2.5}^{+8.0}$)	$2.7\pm 0.1$
HESS J1632-478	$1.4\pm 0.4$	$119\pm 34$	${8.7}_{-3.9}^{+9.1}$(${9.1}_{-3.7}^{+7.9}$)	$1.9\pm 0.1$
HESS J1718-385	$0.21\pm 0.05$	$13\pm 3$	${1.4}_{-0.5}^{+0.8}$(${1.4}_{-0.4}^{+0.7}$)	$1.2\pm 0.1$
HESS J1809-193	$0.76\pm 0.12$	$44\pm 7$	${7.3}_{-1.5}^{+2.1}$(${6.8}_{-1.3}^{+1.8}$)	$2.3\pm 0.1$
HESS J1813-126	$0.63\pm 0.40$	$30\pm 24$	${2.4}_{-1.7}^{+6.5}$(${2.8}_{-1.6}^{+6.0}$)	$1.9\pm 0.2$
HESS J1825-137	$1.15\pm 0.09$	$73\pm 6$	${21}_{-3}^{+3}$(${23}_{-2}^{+3}$)	$2.8\pm 0.1$
HESS J1831-098	$0.52\pm 0.26$	$34\pm 16$	${6.0}_{-3.2}^{+10.1}$(${7.1}_{-3.7}^{+9.3}$)	$1.2\pm 0.1$
HESS J1837-069	$0.77\pm 0.09$	$89\pm 11$	${41}_{-6}^{+8}$(${43}_{-6}^{+8}$)	$2.6\pm 0.1$
HESS J1841-055	$2.50\pm 0.42$	$180\pm 30$	${93}_{-35}^{+84}$(${85}_{-28}^{+53}$)	$2.8\pm 0.1$
HESS J1857+026	$0.58\pm 0.16$	$64\pm 18$	${23}_{-7}^{+11}$(${23}_{-7}^{+9}$)	$2.9\pm 0.1$
HESS J1858+020	$0.25\pm 0.08$	$20\pm 6$	${2.8}_{-1.4}^{+2.7}$(${2.9}_{-1.3}^{+2.2}$)	$1.8\pm 0.2$
HESS J1908+063	$2.2\pm 1.7$	$91\pm 60$	${32}_{-16}^{+56}$(${32}_{-15}^{+41}$)	$2.7\pm 0.1$
HESS J1912+101	$1.05\pm 0.38$	$85\pm 30$	${43}_{-20}^{+46}$(${49}_{-21}^{+43}$)	$1.8\pm 0.1$
HESS J0835-455	$1.65\pm 0.27$	$8.1\pm 1.3$	${0.84}_{-0.19}^{+0.27}$(${0.85}_{-0.17}^{+0.23}$)	$2.4\pm 0.1$
HESS J1303-631	$0.47\pm 0.04$	$88\pm 8$	${48}_{-7}^{+9}$(${45}_{-6}^{+8}$)	$2.4\pm 0.1$
HESS J1356-645	$0.52\pm 0.17$	$28\pm 9$	${12}_{-4}^{+7}$(${13}_{-4}^{+6}$)	$2.8\pm 0.1$
HESS J1420-607	$0.30\pm 0.04$	$29\pm 4$	${8.7}_{-1.6}^{+2.0}$(${9.1}_{-1.7}^{+2.1}$)	$2.5\pm 0.1$
HESS J1616-508	$0.55\pm 0.09$	$46\pm 7$	${19}_{-4}^{+5}$(${21}_{-3}^{+4}$)	$2.9\pm 0.1$
HESS J1640-465	$0.17\pm 0.01$	$39\pm 3$	${18}_{-2}^{+3}$(${19}_{-1}^{+2}$)	$2.9\pm 0.1$
HESS J1708-443	$0.49\pm 0.23$	$22\pm 10$	${5.4}_{-2.9}^{+6.1}$(${5.2}_{-2.5}^{+5.0}$)	$2.6\pm 0.1$
HESS J1813-178	$0.19\pm 0.02$	$15\pm 2$	${5.0}_{-0.9}^{+1.0}$(${4.6}_{-0.8}^{+0.9}$)	$2.6\pm 0.1$
HESS J1826-130	$1.13\pm 0.46$	$31\pm 13$	${4.9}_{-2.2}^{+4.8}$(${5.5}_{-2.3}^{+4.7}$)	$2.4\pm 0.2$
HESS J1833-105	$0.21\pm 0.07$	$15\pm 5$	${4.6}_{-1.2}^{+3.4}$(${4.7}_{-1.7}^{+3.0}$)	$3.0\pm 0.2$
Geminga	$5.5\pm 0.7$	$24$	${5.0}_{-1.0}^{+2.0}$(${2.1}_{-0.7}^{+1.0}$)	2.3
Monogem	$4.8\pm 0.6$	$24$	${25}_{-2.1}^{+3.3}$	2.1

Equations

$$Q(E,t)=L(t){\left(\frac{E}{E_0}\right)}^{-\gamma }\exp \left(-\frac{E}{E_c}\right),$$

$$L(t)=\frac{L_0}{{\left(1+\frac{t}{{\tau }_0}\right)}^2},$$

$$E_{\mathrm{tot}}={\int }_0^{T}dt{\int }_{E_1}^{\infty }dEEQ(E,t),$$

$$W_0={\tau }_0\dot{E}{\left(1+\frac{T}{{\tau }_0}\right)}^2.$$

$${\mathcal{N}}_e(E,\mathbf{r},t)$$

$$\Delta \tau (E,E_s)\equiv {\int }_E^{E_s}\frac{d{E}^{\prime }}{b(E^{\prime })}=t-t_{\mathrm{obs}}.$$

$${\lambda }^2={\lambda }^2(E,E_s)\equiv 4{\int }_E^{E_s}d{E}^{\prime }\frac{D(E^{\prime })}{b(E^{\prime })},$$

$${\Phi }_{e^{\pm }}(E)=\frac{c}{4\pi }{\mathcal{N}}_e(E,|\mathbf{r}-{\mathbf{r}}_{\mathbf{s}}|=d,t=T).$$

$$\displaystyle D(E,r)=\left\{\begin{array}[]{rl}&D_{0}(E/1{\rm\,GeV})^{\delta}{\rm\;for\;}0<r<r_{b},\\ &D_{2}(E/1{\rm\,GeV})^{\delta}{\rm\;for\;}r\geq r_{b}.\end{array}\right.$$

$${\Phi }_{\gamma }(E_{\gamma },\Delta \Omega )={\int }_{m_e{c}^2}^{\infty }dE\mathcal{M}(E,\Delta \Omega ){\mathcal{P}}^{\mathrm{IC}}(E,E_{\gamma }).$$

$$\mathcal{M}(E,\Delta \Omega )={\int }_{\Delta \Omega }d\Omega {\int }_0^{\infty }ds{\mathcal{N}}_e(E,s,T).$$

$${\Phi }_{\gamma }^{68\%}(E_{\gamma })=2\pi {\int }_0^{{\theta }_{68}}\frac{d{\Phi }_{\gamma }}{d\theta }(E_{\gamma })\sin \theta d\theta ,$$

$${\theta }_{\mathrm{motion}}(E_{\gamma })=\mathrm{atan}\left(\frac{{\mathrm{d}}_{\mathrm{motion}}({\mathrm{E}}_{\gamma })}{\mathrm{d}}\right),$$

$$d_{\mathrm{motion}}(E_{\gamma })=\frac{v_T{E}_e(E_{\gamma })}{b(E_e(E_{\gamma }))}.$$

$${\theta }_{\mathrm{motion}}(E_{\gamma })[\mathrm{deg}]=\mathrm{atan}\left(0.324\frac{\frac{{\mathrm{v}}_{\mathrm{T}}[km/s]}{\mathrm{b}[10^{-16}GeV/s]({\mathrm{E}}_{\mathrm{e}}[\mathrm{GeV}]])}}{\mathrm{d}[\mathrm{kpc}]}\right).$$

$$\frac{d{\Phi }_{\gamma }}{d\theta }\sim \frac{1}{{\theta }_{\mathrm{ICS}}(\theta +0.06\cdot {\theta }_{\mathrm{ICS}})}{e}^{-{\left(\frac{\theta }{{\theta }_{\mathrm{ICS}}}\right)}^2},$$

Full text

I Introduction

A $\gamma$-ray emission at TeV energies and of a few-degrees extension size in the direction of Geminga and Monogem pulsar wind nebulae (PWNe) has been detected by HAWC Abeysekara et al. (2017a) and Milagro Abdo et al. (2009). The presence of a $\gamma$-ray halo around Geminga has been recently confirmed by Di Mauro et al. (2019) with an analysis of Fermi-LAT data above 8 GeV, whose extension reaches about 15 degrees at 10 GeV. The $\gamma$-ray halos detected around Geminga and Monogem are interpreted as photons produced by electrons ($e^{-}$) and $e^{+}$ ($e^{+}$) accelerated by their PWNe and inverse Compton scattering (ICS) low-energy photons of the interstellar radiation fields (ISRFs). These observations may give us the possibility to shed light on the origin of the $e^{+}$ excess in cosmic rays (CR), firstly detected by Pamela Adriani et al. (2013), then by Fermi-LAT Ackermann et al. (2012) and recently, with an unprecedented precision, by AMS-02 Aguilar et al. (2019). The extension of detected $\gamma$-ray halos suggests that the diffusion around these PWNe is about two orders of magnitude less intense than the value assumed to fit the latest CR data measured by AMS-02 (see, e.g. Kappl et al. (2015); Genolini et al. (2015, 2019)). The inferred diffusion coefficient is in fact of about $10^{27}$ cm2/s at 100 GeV Abeysekara et al. (2017a); Di Mauro et al. (2019).

The ICS halos detected around Geminga and Monogem are called by some authors “TeV halos”, since they have been mainly detected at very-high-energy (VHE) (see, e.g., Linden et al. (2017)). However, we will refer to them as “ICS halos” because of the recent detection of the Geminga halo at Fermi-LAT energies, and because we prefer to characterize this emission with the physical process that generates it, and not with the energy at which it is detected. It is still unclear if these halos are generated by $e^{\pm}$ accelerated by PWNe and diffusing in the interstellar medium (ISM), or by $e^{\pm}$ propagating in a region still dominated by the PWN environment. Very recently, Ref. Giacinti et al. (2019) investigated this point by using a sample of Galactic PWNe taken from the HESS survey of the Galactic plane (HGPS) Abdalla et al. (2018a). They have estimated the $e^{\pm}$ density at the location of the source VHE $\gamma$-ray emission. Comparing this density with the one of the ISM, they concluded that for most of these sources, except for Geminga and Monogem, the $e^{\pm}$ are probably still confined in the PWN. Therefore, they call these sources $e^{\pm}$ halo, rather than TeV halo. Their calculation is based on a series of assumptions, such as the shape of the $e^{\pm}$ injection spectrum, the energy range for accelerated $e^{\pm}$, and no time dependence considered for the spin-down luminosity. Also, the size of the ICS halos is taken directly from the HESS catalog. Changing some of these assumptions their results might change significantly and many of the sources in their sample could have a density of $e^{\pm}$ of the same order of the ISM. This would imply that these cosmic particles might not be confined in the PWN. We will discuss in Sec. VI how their results would change assuming the size of the ICS emission as estimated in this paper.

The detection of ICS halos around Geminga and Monogem can provide key information about the acceleration mechanisms of $e^{\pm}$ from PWNe, and their propagation in the Galactic environment. For example Ref. Di Mauro et al. (2019) used the flux and morphology of the ICS halo detected from Geminga and found that this source contributes at most $10\%$ to the $e^{+}$ excess. They have also found evidences for a low-diffuse bubble located around the pulsar, with a size of around 100 pc and a value of the diffusion coefficient at 1 GeV of about $2.3\times 10^{26}$ cm2/s, i.e., two orders of magnitude lower than the average of the Galaxy. Several references (see, e.g., Hooper et al. (2017); Abeysekara et al. (2017a); Xi et al. (2018); Tang and Piran (2018); Fang et al. (2018); Di Mauro et al. (2019)) have studied the flux of $e^{+}$ from PWNe in light of the Milagro and HAWC data, and have drawn conclusions on the contribution of this source population to the $e^{+}$ excess. Reference Abeysekara et al. (2017a) uses the low diffusion found around Geminga and Monogem PWNe to propagate particles in the entire Galaxy, and claims their contribution is negligible. On the other hand the authors of Hooper et al. (2017); Fang et al. (2018) claim Geminga explains most of the $e^{+}$ data. Finally, references Tang and Piran (2018); Di Mauro et al. (2019) agree on the fact that the contribution of Geminga is at the $10\%$ level. Although, most of these papers suggest PWNe are likely the main contributors to the $e^{+}$ flux, they use the results based on only those two PWNe. Indeed, we still do not have a large enough sample of ICS halos and we have not collected evidences if such a low-diffusion bubble is present or not around a significant sample of Galactic pulsars.

In addition to Geminga and Monogem, many more ICS $\gamma$-ray halos could have been already detected in the direction of other Galactic pulsars by Imaging Atmospheric Cherenkov Telescopes (IACTs), HAWC, MILAGRO and Fermi-LAT. The HAWC Collaboration has recently released the 2HWC catalog Abeysekara et al. (2017b) which contains 39 sources detected close to the Galactic plane. Many of them have an extended $\gamma$-ray morphology, and are spatially close to powerful Galactic pulsars. The HESS Collaboration has recently published the results of a the HGPS catalog which is the most comprehensive survey of the Galactic plane in VHE $\gamma$ rays. This publication includes Galactic sky maps and the catalog with the properties of the 78 sources Abdalla et al. (2018a). Many of these sources have been detected as extended, and are probably associated to PWNe. Therefore, the 2HWC and HGPS catalogs represent two rich datasets for investigating the acceleration mechanism of $e^{+}$ from PWNe, and their diffusion around those sources. Following the detection of Geminga and Monogem extended halos, the possible presence of similar objects around other Galactic pulsars has been explored in Linden et al. (2017), by assuming a ”Geminga-like” TeV halo for each considered pulsar. Moreover, the expected number of ICS halos detectable by current and future observatories has been estimated in Sudoh et al. (2019). As discussed in the rest of the present paper, we significantly extend the current literature by performing a complete calculation of the ICS flux for each considered source. Furthermore, we present a novel analysis of the data provided by HGPS catalog to characterize the observed gamma-ray emission around many Galactic pulsars, in the light of the presence a possible ICS halo.

In the first part of this paper we will inspect how the extension of the ICS halo in PWNe depends on the age and distance of the host pulsar, and on the intensity of the diffusion coefficient present around them. The ICS halo size is a key parameter for IACTs which have a limited instantaneous field of view of $4-5^{\circ}$. Then, we will show how the ICS halo size depends on the extension of the low-diffusion bubble and the pulsar proper motion. In fact, pulsars have an average proper motion of 100 km/s Hobbs et al. (2005) and, as we have shown in Di Mauro et al. (2019), this effect distorts the ICS $\gamma$-ray morphology. In the present study, we argue that the most promising energy range for searching for ICS halos is above 100 GeV, where IACTs, and HAWC and Milagro operate. We then use the ICS flux to predict the brightest pulsars around which HAWC should detect an ICS halo. Finally, we predict the number of ICS halos detectable by HESS, HAWC and in the future by the Cherenkov Telescope Array (CTA) Acharya et al. (2017).

In the second part, we consider the PWNe already detected by IACTs. In particular, we use a sample of sources associated to PWNe or PWNe candidates taken from the 2HWC and HGPS catalogs. We use their measured size and flux to determine the diffusion coefficient around each source, and to estimate the minimal dimension of the low-diffusion bubble.

The paper is organized as follows. In Sec. II we present our model for the acceleration of $e^{\pm}$ from PWNe, $e^{\pm}$ propagation in the Galaxy and the flux of $\gamma$ rays for ICS. In Sec. III we investigate how the ICS halos size depends on the pulsar distance, age and proper motion, and how it changes according to the diffusion coefficient. In Sec. IV we study the detectability of ICS halos at IACTs and rank the pulsars in ATNF catalog Manchester et al. (2005) according to their expected ICS halo brightness. Sec. V contains the methodology employed for the derivation of the diffusion coefficient around the sources, whose results are presented in Sect.  VI. We draw our conclusions in Sec. VII.

II Model for the $e^{\pm}$ and $\gamma$-ray emission from a PWN

We recall here the basics for modeling the $e^{\pm}$ and the consequent ICS $\gamma$-ray emission from PWNe. We follow the formalism detailed in Di Mauro et al. (2019).

PWNe are thought to accelerate and inject $e^{\pm}$ in the ISM up to VHE (see, e.g., Chi et al. (1996); Amato (2014); Gaensler and Slane (2006)). A rapidly spinning neutron star, or pulsar, formed after a supernova explosion, is likely the engine of this process. The rotation of the pulsar induces an electric field that extracts $e^{-}$ from the star surface. These $e^{-}$ lose energy via curvature radiation while propagating far from the pulsar along the magnetic field lines, and the energetic emitted photons create a wind of $e^{\pm}$ pairs in the intense neutron star magnetic field.

For the sake of completeness, we here briefly recall the basic understandings of the PWN evolution, which is then treated effectively. According to Chi et al. (1996); Amato (2014); Gaensler and Slane (2006), the initial phase of the PWN evolution, called free expansion phase, occurs in the first few thousands of years. At this stage, the pulsar wind is surrounded by the expanding shell of the supernova remnant (SNR), which moves at a speed of about $5-10\cdot 10^{3}$ km/s, while the pulsar located at the center of the SNR has a velocity of the order of $400-500$ km/s. The expansion velocity of the pulsar wind increases constantly with time, and the size $R$ of the PWN goes as $R\propto t^{1.2}$ Chevalier (1977). During the free expansion, the pulsar wind expands very fast while the SNR ejecta interacts with the ISM creating a forward and reverse shock. The PWN reaches, at this stage, a size of about 10 pc.

After a few thousands years, the reverse shock moves inward and interacts with the outward moving PWN shock. This interaction constitutes a termination shock, and its bulk energy is dissipated into a relativistically, magnetized fluid, which shines as a PWN. The total energy of the SNR exceeds the one of the PWN by one or two orders of magnitude, so that the PWN can be compressed by up to a factor of 10 Gelfand et al. (2009). During this process the PWNe experiences a series of contractions and expansions until a steady balance is reached. Once the reverberations between the PWN and the SNR reverse shock have faded, the pulsar can again power a bubble steadily expanding as $R\propto t^{1.2}$ for $t<\tau_{0}$ and $R\propto t^{0.3}$ for $t>\tau_{0}$, where $\tau_{0}$ is the pulsar decay time Reynolds and Chevalier (1984); van der Swaluw et al. (2001). Therefore, at a time larger than $\tau_{0}$ the PWN size is not expected to have a strong evolution with the pulsar age. The $e^{\pm}$ pairs produced in the pulsar magnetosphere reach the termination shock and, due to the severe energy losses, their energy is at most a few tens of GeV. The termination shock is the place where particle acceleration eventually occurs, and a relatively large fraction (up to few tens of percent) of the wind bulk energy is converted into accelerated pairs. They then radiate into a photon spectrum extending from radio frequencies to TeV $\gamma$-rays, through synchrotron and ICS processes.

Given the initial velocity, the distance traveled by the pulsar from the explosion site after few tens of kyr can be comparable to or even larger than the radius of an equivalent spherical PWN around a stationary pulsar. The pulsar thus can abandon its original wind bubble, leaving behind itself a relic PWN, and generating a new, smaller PWN around its current position, which is called bow shock. Observationally, this appears as a central, possibly distorted PWN visible in radio and X-ray and powered by freshly accelerated $e^{\pm}$. The relic PWN is powered by $e^{\pm}$ injected along its formation history.

The PWNe considered in this paper are older than a few thousands of year. Therefore, these PWNe have probably already interacted with the reverse shock of the SNR. Moreover, the $e^{\pm}$ accelerated by younger sources could be still confined inside the PWN or the SNR, while for older sources they have been probably injected from the relic and bow shock components of the PWN, and released in the ISM environment. In order to inspect any dependence of our results by the presence of the SNR and PWN environment, we select PWNe powered by pulsars of different ages from a few to hundreds of kyr.

We consider a model in which $e^{\pm}$ are continuously injected with a rate that follows the pulsar spin-down energy, i.e. a continuous injection scenario. This scenario is indeed required to generate the TeV photons detected by Milagro and HAWC for Geminga and Monogem Yuksel et al. (2009); Abeysekara et al. (2017a); Di Mauro et al. (2019). A common alternative is to consider a burst like scenario, where all the particles are emitted from the source at a time equal to the age of source $T$. In our model, the injection spectrum $Q(E,t)$ for the accelerated $e^{\pm}$ pairs is assumed to effectively describe the particles that are produced during the acceleration process and released in the ISM, while no attempt is made to describe the dynamical evolution during the first thousands of years of the PWN, or possible modification in the spectrum of particles during the release processes. The injection spectrum $Q(E,t)$ can be effectively described by a power law with an exponential cutoff:

[TABLE]

where the magnetic dipole braking $L(t)$ (assuming a magnetic braking index of 3) is defined as:

[TABLE]

and $\tau_{0}$ is the characteristic pulsar spin-down timescale. The cutoff energy $E_{c}$ is fixed to $10^{3}$ TeV. We set $\tau_{0}=12$ kyr if not stated otherwise, following Aharonian et al. (1995); Malyshev et al. (2009); Abeysekara et al. (2017a); Di Mauro et al. (2019). A smaller value of $\tau_{0}$, as derived for example from fits to radio, X-ray and $\gamma$-ray data to few, very young pulsars Torres et al. (2014), would have the consequence to lower the $\gamma$-ray flux at high-energy and so higher efficiency values would be found. The injection spectrum of $e^{\pm}$ is usually measured with a broken power-law spectrum (see, e.g., Malyshev et al. (2009)) with a break at energies between a few up to hundreds of GeV. The $\gamma$-ray energies we will consider are beyond a few hundreds of GeV so for $e^{\pm}$ well beyond the TeV energies. Therefore, the injection spectrum of particles from PWNe can be modeled in our care with a simple power-law and neglecting the presence of a break.

The total energy emitted by the source in $e^{\pm}$ is given by:

[TABLE]

where we fix $E_{1}=0.1$ GeV Buesching et al. (2008); Sushch and Hnatyk (2014). $E_{{\rm tot}}$ is related to the pulsar total spin down energy $W_{0}$ with the following relation $E_{{\rm tot}}=\eta W_{0}$, where $\eta$ is the fraction of the pulsar’s spin-down luminosity which goes into $e^{\pm}$ particles. $W_{0}$ can be computed from catalogued quantities as the pulsar age $T$, the decay time $\tau_{0}$, and the spin-down luminosity $\dot{E}$:

[TABLE]

The normalization of the injection spectrum ($L_{0}$, see Eq. 1) is found using Eq. 3 and assuming that $E_{{\rm tot}}=\eta W_{0}$.

The source term in Eq. 1 is inserted in a diffusion-loss equation to compute the $e^{\pm}$ number density $\mathcal{N}_{e}(E,\mathbf{r},t)$ per unit volume and energy, and at an observed energy $E$, a position $\mathbf{r}$, and time $t$. We account for space-indipendent energy losses $b(E)$ by means of synchrotron and inverse Compton processes. The interstellar photon populations at different wavelengths have been taken from Vernetto and Lipari (2016). The Galactic magnetic field intensity has been assumed $B=3.6\;\mu$G, as resulting from the sum (in quadrature) of the regular and turbulent components (Sun et al., 2007). For more details on the propagation model, we address to Di Mauro et al. (2019) and references therein. In the continuous injection scenario and with a homogeneous diffusion in the Galaxy, the solution for the $e^{\pm}$ number density $\mathcal{N}_{e}(E,\mathbf{r},t)$ at an observed energy $E$, position $\mathbf{r}$ and time $t$ is given by Yuksel et al. (2009):

[TABLE]

where the integration over $t^{\prime}$ accounts for the PWN releasing $e^{\pm}$ continuously in time. The energy $E_{s}$ is the initial energy of $e^{\pm}$ that cool down to $E$ in a loss time $\Delta\tau$:

[TABLE]

The $b(E)$ is the energy loss function, $\mathbf{r_{s}}$ indicates the source position, and $\lambda$ is the typical propagation scale length defined as:

[TABLE]

with $D(E)$ the diffusion coefficient. As a matter of fact, the tipycal propagation time, defined as $\lambda^{2}/4D(E)$ sets grossly in the range $10^{3}-10^{5}$ kyr, decreasing with $E$. The flux of $e^{\pm}$ at Earth is given by:

[TABLE]

Recent results Abeysekara et al. (2017a); Di Mauro et al. (2019) suggest that the diffusion coefficient around Geminga and Monogem PWNe is $\sim 10^{26}$ cm2/s at 1 GeV, i.e. about two orders of magnitude smaller than the value derived for the entire Galaxy through a fit to AMS-02 CR data Kappl et al. (2015); Genolini et al. (2015, 2019). A phenomenological description for this discrepancy proposes a two-zone diffusion model, where the region of low diffusion is contained around the source, and delimited by an empirical radius $r_{b}$ Profumo et al. (2018); Tang and Piran (2018). The inhibition of diffusion near pulsars has been recently discussed in Evoli et al. (2018), and a possible theoretical interpretation is provided. This paper predicts a very strong dependence of the diffusion coefficient as a function of the pulsar age with $D_{0}\sim 10^{26}$ cm2/s at 1 GeV for sources with $T\sim 20$ kyr and values close to the average of the Galaxy for $T>100$ kyr. Nevertheless, a conclusive understanding of this phenomenon is not yet achieved, and the analysis we present in this paper can give new insights on the theoretical models.

In this paper we include the phenomenological two-zone diffusion model as in Tang and Piran (2018); Di Mauro et al. (2019) to account for these recent observations, for which the diffusion coefficient is defined as:

[TABLE]

Here $r_{b}$ is the size of the low-diffusion bubble while $D_{0}$ and $D_{2}$ are the diffusion coefficients inside and outside the bubbles, respectively.

The $e^{\pm}$ accelerated by PWNe can produce photons whose energy covers a wide range (see, e.g., Slane (2017) for a recent review). From radio to X-ray energies, photons are produced by $e^{\pm}$ through synchrotron radiation caused by the magnetic fields. On the other hand, at higher energies $\gamma$ rays are produced from VHE $e^{\pm}$ escaped from the PWN by the ICS off the ISRF. We are interested here in the extended halo emission of the size of at least tens of arcminutes (i.e., around tens of parsec) generated by $e^{\pm}$ injected by PWNe in the Galactic environment, and not to the small-scale structures extended between few arcseconds to arcminutes and observed in the nebula, as for example jets and torii (see e.g. Posselt et al. (2017)).

The ICS photon flux emitted by a PWN, at a $\gamma$-ray energy $E_{\gamma}$ and within a solid angle $\Delta\Omega$, can be computed as Blumenthal and Gould (1970); Cirelli et al. (2011):

[TABLE]

The term $\mathcal{M}(E,\Delta\Omega)$ represents the spectrum of $e^{+}$ and $e^{-}$ of energy $E$ propagating in the Galaxy and from a solid angle $\Delta\Omega$:

[TABLE]

$\mathcal{N}_{e}(E,s,T)$ is the energy spectrum of $e^{\pm}$ taken from Eq. 5, $s$ is the line of sight, and $\mathcal{P}^{\rm IC}(E,E_{\gamma})$ is the power of photons emitted by a single $e^{\pm}$ by ICS, defined as in Blumenthal and Gould (1970); Delahaye et al. (2010). The ICS occurs off the CMB, described by a blackbody energy density ($T_{{\rm CMB}}=2.753$ K), the infrared light (peaked at $T_{{\rm IR}}=3.5\cdot 10^{-3}$ eV) and by the starlight ($T_{{\rm SL}}=0.3$ eV) Vernetto and Lipari (2016); Porter et al. (2006); Popescu et al. (2017). Our results are obtained for the ISRF energy density in the local Galaxy reported in Vernetto and Lipari (2016), but we have explicitly checked that they do not get modified by using the model in Porter et al. (2006).

As shown in Di Mauro et al. (2019), the proper motion of the pulsar could alter the morphology of the $\gamma$-ray ICS halo. The proper motion affects significantly the morphology of the $\gamma$-ray emission for pulsars older than about 100 kyr and moving with a velocity of at least 100 km s*-1*. This is particularly true for Geminga, that is a very close pulsar ($d=250$ pc), has a transverse proper motion of $v_{T}\approx 211$ km s*-1* Faherty et al. (2007) and $T=340$ kyr. We implement this effect in Eq. 5 by replacing $\bf{r_{s}}$ with ${\bf{v_{T}}}t$, where $\bf{v_{T}}$ is the vector of the transverse velocity.

As for the ICS photon flux emitted by a PWN, our benchmark is the one-zone diffusion model in which, effectively, $r_{b}\rightarrow\infty$, and the $D_{0}$ corresponds to the low diffusion coefficient around the PWN. Using the one-zone diffusion model for the ICS is appropriate since most of the $\gamma$-ray emission is generated close to the pulsar where the low diffusion probably acts. On the other hand, for the calculation of the $e^{+}$ flux at Earth the two-zone model must be considered since the size of the low-diffusion zone around the PWN is much smaller than the propagation volume from the source to the Earth. We have already applied these choices in Di Mauro et al. (2019). We will also discuss some examples in which the ICS photon flux is computed in a two-zone diffusion model.

III Angular size of the $\gamma$-ray ICS halos

In this section we study the size of ICS halos, defined through the $\gamma$-ray flux, for different values of $E_{\gamma}$, and as a function of the strength ($D_{0}$) and size ($r_{b}$) of the low-diffusion bubble, the age and distance of the host pulsar, and of its proper motion, in order to motivate the selection of pulsars used in Sec. V.1.

The size is a key parameter for the detectability of ICS halos. Indeed, IACTs have a few degrees instantaneous field of view and a very extended halo would be difficult to detect. It is also challenging to detect an halo with a size larger than about $10^{\circ}$ with Fermi-LAT data, because below 100 GeV the interstellar emission is by far the major contributor of the observed flux, and an imperfect modeling of this component could produce spurious residuals and unreliable results.

We define the size of an ICS halo as the angle $\theta_{68}$ which contains the $68\%$ of the flux $\Phi_{\gamma}$:

[TABLE]

where $d\Phi_{\gamma}/d\theta$ is the surface brightness and is computed from Eqs. 12 and 13, where $\Delta\Omega$ depends on the angle $\theta$ from the center of the pulsar. This formulation of the ICS halo size follows the definition of the 68% containment radius, used by $\gamma$-ray experiments to define the size of extended sources (see, e.g., Acero et al. (2015)).

We first investigate how $\theta_{\rm{68}}$ changes according to the pulsar distance and age, the diffusion coefficient and $\gamma$-ray energy. We calculate $\theta_{\rm{68}}$ for a grid of pulsar distance and age values between $d\in[0.1,10]$ kpc and $T\in[10,10^{4}]$ kyr, repeated for three values of $D_{0}$: $6\cdot 10^{25}$ cm2/s, $2\cdot 10^{26}$ cm2/s and $1\cdot 10^{27}$ cm2/s. The first two values of $D_{0}$ are inspired to the results for Geminga found in Abeysekara et al. (2017a); Di Mauro et al. (2019) while the third one has been set to a value closer to the average Galactic diffusion. Finally, we repeat this exercise for $E_{\gamma}=10$ GeV, which is relevant for Fermi-LAT data, and $E_{\gamma}=1$ TeV, where the IACTs have their peak of sensitivity.

We show our results for $\theta_{\rm{68}}$ in Fig. 1 and 2, where we superimpose the ATNF catalog pulsars. We notice that $\theta_{\rm{68}}$ is significantly smaller at $E_{\gamma}=1$ TeV than at 10 GeV for sources older than about 200 kyr. Indeed, for such old sources VHE $e^{\pm}$ lose energy very quickly, so that the ICS $\gamma$-ray emission is much closer to the pulsar location. This trend is confirmed by the recent detection of the Geminga ICS halo with a size of about $5^{\circ}$ above 5 TeV Abeysekara et al. (2017a) and about $15^{\circ}$ at 10 GeV Di Mauro et al. (2019). On the other hand, sources younger than about 200 kyr have extension at 1 TeV that is slightly larger than the one at 10 GeV because for these ages 1 TeV $e^{\pm}$ have a propagation length $\lambda$ (see Eq. 7) that is larger than the one at 10 GeV.

We also observe that the larger is $D_{0}$ the larger is $\theta_{\rm{68}}$. For example, for a source as Geminga with $d=0.25$ pc and $T=340$ kyr and at $E_{\gamma}=10$ GeV ($E_{\gamma}=1$ TeV) the size of $\theta_{\rm{68}}$ is $15^{\circ}$, $25^{\circ}$ and $30^{\circ}$ ($10^{\circ}$, $18^{\circ}$ and $25^{\circ}$) for $D_{0}$ equal to $6\cdot 10^{25}$ cm2/s, $2\cdot 10^{26}$ cm2/s and $1\cdot 10^{27}$ cm2/s, respectively. A higher diffusion coefficient makes the particle travel a larger distance in the Galaxy before losing most of its energy.

IACTs have an instantaneous field of view between $3.5-5^{\circ}$, thus if $D_{0}=10^{27}$ cm2/s only sources farther than about 3 kpc would have a detectable ICS halo. On the other hand, if $D_{0}\sim 10^{26}$ cm2/s, as detected for Geminga in Abeysekara et al. (2017a); Di Mauro et al. (2019), most of the ATNF catalog pulsars would be good targets for ICS halo searches by IACTs. Instead, in the Fermi-LAT energy range most of Galactic pulsars would generate very extended halos. More precisely, fixing $D_{0}=6\cdot 10^{25}$ cm2/s ($D_{0}=2\cdot 10^{26}$ cm2/s), the size of $\theta_{\rm{68}}$ would be smaller than two degrees only for $d\geq 10^{0.58\log_{10}{(T\,[\,\rm{kyr}])}-1.2}$ kpc ($d\geq 10^{0.57\log_{10}{(T\,[\,\rm{kyr}])}-0.9}$ kpc). This means that a source with an age of 100 kyr should be farther than about 0.9 kpc (1.7 kpc) if $D_{0}=6\cdot 10^{25}$ cm2/s ($D_{0}=2\cdot 10^{26}$ cm2/s) to be detected with an extension smaller than two degrees.

These results may change if a two-zone diffusion model is considered. In this model, the pulsar is located into a bubble of low diffusion where $e^{\pm}$ are more effectively confined. In general, assuming a two-zone diffusion model has the effect of increasing $\theta_{\rm{68}}$. In Fig. 3 we show the surface brightness $d\Phi_{\gamma}/d\theta$ calculated for a pulsar with $d=1$ kpc and $T=100$ kyr at $E_{\gamma}=1$ TeV, and assuming either a one or a two-zone diffusion model, where $D_{0}=2\cdot 10^{26}$ cm2/s only within $r_{b}=50$ pc. It is clear from the figure that the two-zone diffusion model has a much wider profile, which results into a more extended ICS flux. This effect depends on the value of $r_{b}$.

In Fig. 4 we study $\theta_{\rm{68}}$ as a function of $r_{b}$, for a pulsar with $d=1$ kpc and $T=100$ kyr, and an other one with $d=2$ kpc and $T=60$ kyr. For these two cases we set $D_{0}=6\cdot 10^{25}$ cm2/s, $2\cdot 10^{26}$ cm2/s and $1\cdot 10^{27}$ cm2/s. For $r_{b}\geq 0.1$ kpc, $\theta_{\rm{68}}$ tends to the value obtained with the one zone model (see Fig. 1). Indeed, for such a large low-diffusion zone bubble most of the $e^{\pm}$ lose completely their energy before reaching the high-diffusion zone. Therefore, they are completely trapped inside the low-diffusion bubble. This effect could be a result of the confinement of the CRs inside the PWN and/or the PWN. For example, in the case of the pulsar with $d=1$ kpc and $T=100$ kyr, $\gamma$ rays with energies of $E_{\gamma}=1$ TeV are on average produced by $e^{\pm}$ of 10 TeV energy. These $e^{\pm}$ in a diffusion environment with $D_{0}=6\cdot 10^{25}$ cm2/s have a propagation length of about 30 pc. Therefore, if $r_{b}$ is larger than this length, $\theta_{\rm{68}}\approx 2^{\circ}$, similar to the value found for the one-zone model. On the other hand, for smaller values of $r_{b}$, $e^{\pm}$ exit the low-diffusion zone before losing most of their energy and produce a significant ICS flux in the high-diffusion zone. Since outside the low-diffusion bubble $e^{\pm}$ travel significant larger distances, the ICS halo can become very extended. We also notice that the lower is $D_{0}$ the lower is the value of $r_{b}$ at which we observe the transition between small and large values of $\theta_{\rm{68}}$. This is due to the fact that with a less intense $D_{0}$, $e^{\pm}$ travel shorter distances before losing most of their energies. We conclude that for $D_{0}\sim 10^{26}$ cm2/s values of $r_{b}\geq 80$ pc do not alter significantly $\theta_{\rm{68}}$. In other words, diffusion coefficient values of the order of $D_{0}\sim 10^{26}$ cm2/s with $\theta_{\rm{68}}$ at the degree scale implies $r_{b}\gtrsim 80$ pc.

III.1 Pulsar proper motion

An other element that can affect the spatial morphology of $\gamma$ rays produced for ICS is the pulsar proper motion. The analysis presented in Di Mauro et al. (2019) shows that the proper motion of the Geminga pulsar, which moves with a transverse velocity $v_{T}=$ 211 km/s Faherty et al. (2007), shapes significantly the $\gamma$-ray ICS template below 100 GeV. In particular we have shown that at 10 GeV the ICS halo has a distortion of about $10^{\circ}$ in the opposite direction of the proper motion. This is due to the fact that 10 GeV $\gamma$ rays are produced by $e^{\pm}$ emitted by the pulsar tens of kyr ago. Therefore, a significant fraction of the $\gamma$-ray flux is detected in the direction where the pulsar was in the past.

The ICS power $\mathcal{P}^{ICS}$ has a peak at around $E_{e}=1.5$ TeV for $E_{\gamma}=10$ GeV and $E_{e}=60$ TeV for $E_{\gamma}=10$ TeV. We use here the ISRF model as in Vernetto and Lipari (2016). Very similar results are found with the model presented in Refs. Porter et al. (2006); Popescu et al. (2017). An electron of energy of 1.5 TeV (60 TeV) loses most of its energy after about 300 kyr (20 kyr). In this time lapse the Geminga pulsar has travelled across the sky for 60 pc (4 pc). Therefore, we expect that the size of extension of the ICS halo is distorted by about $12^{\circ}$ ($0.9^{\circ}$) in the opposite direction of the proper motion (see Fig. 10 in Di Mauro et al. (2019)).

We generalize this calculation and derive the source distance and age values for which the proper motion is a relevant effect in the ICS flux. The angular size $\theta_{\rm{motion}}$ by which the ICS halo is distorted due to the pulsar proper motion can be parametrized as:

[TABLE]

where $d$ is the actual distance of the source from Earth and $d_{\rm{motion}}$ is:

[TABLE]

Here $E_{e}(E_{\gamma})$ is the energy of the electron for which the ICS power $\mathcal{P}^{ICS}$ has its peak for a given $\gamma$-ray energy and $v_{T}$ is the transverse velocity of pulsar. We can now put together Eq. 15 and 16 finding:

[TABLE]

In Fig. 5 we show the value of $\theta_{\rm{motion}}$ for $\gamma$-ray energy of 10 GeV and 1 TeV. Here we assume energy losses for ICS and synchrotron radiation parametrized as $b(E)=5\times 10^{-17}\rm{GeV/s}\,(E_{e}[\,\rm{GeV}])^{2}$. The angular distortion at 1 TeV is significantly smaller with respect to the 10 GeV case. Indeed, at 1 TeV $\theta_{\rm{motion}}>1^{\circ}$ only for pulsars with velocities larger than about 300 km/s and closer than a few hundred pc. For all other $v_{T}$-$d$ combinations the angular distortion is not significant.

In Abdalla et al. (2018b) the HESS Collaboration found that the offset between the PNW $\gamma$-ray emission and the central pulsar is between $0.2-0.4^{\circ}$. From Fig. 5 this would be consistent with pulsar proper motion with velocities $v_{T}$ smaller that few hundred km/s. Indeed, most of the pulsars have velocities of the order of 100 km/s (see, e.g., Hobbs et al. (2005) for a compilation of pulsar proper motion measurements).

On the other hand, at $E_{\gamma}=10$ GeV even moderate pulsar velocities affect the morphology of the ICS emission, implying $\theta_{\rm{motion}}$ of at least a few degrees. This represents a limiting factor for detecting ICS halos in Fermi-LAT data, since $v_{T}$ is known only for a limited number of pulsars (about 230 over almost 3000 detected so far). Indeed, performing a search for ICS emission from a pulsar with unknown $\bf{v_{T}}$ is challenging, since the intensity and direction of the motion can create a significant asymmetry in the morphology. This issue is probably alleviated by the fact that the most promising pulsars for the ICS halo search are also the better observed and studied and for many of them the proper motion has been already measured.

IV Inverse Compton Scattering halos at TeV energies

In this section we illustrate how the $\gamma$-ray flux selects the most promising pulsars with a detectable ICS halo. First, we predict the number of ICS halos that could be detected by HAWC, HESS and CTA as a function of the efficiency $\eta$. The number of expected ICS halos detections with HAWC, HESS and CTA has been recently calculated in Sudoh et al. (2019). Their model uses different assumptions with respect to ours. In particular, instead of computing the extended ICS flux for each source, they rescale the observed Geminga gamma-ray flux to all the sources, assuming they are ”Gemiga-like” systems. Moreover, their results for the cumulative number of detections vary by about one order of magnitude according to the choice of the pulsar rotational period and the alignment of the pulsar jet. Therefore, their results are not easily comparable with ours.

IV.1 IACTs detectability of extended ICS halos

HAWC. The 2HWC catalog Abeysekara et al. (2017b) reports the sensitivity for the detection of a point source as a function of the declination. The lowest detectable flux at 7 TeV is $6\cdot 10^{-15}({\rm TeV\,cm^{2}\,s})^{-1}$ for declination angles in the range $10^{\circ}-30^{\circ}$ and a point source with a spectral slope of $-2.5$. However, this value is not appropriate for our scope, because we are interested in the detection of extended ICS halos with a size of a fraction of the degree (see Tab. 1). We estimate the HAWC sensitivity to ICS halos by taking the publicly available data of the 2HWC Survey111https://data.hawc-observatory.org/datasets/2hwc-survey/index.php. This on-line resource provides - at each direction in the HAWC field of view - the significance for the presence of a source, for different spatial and spectral assumptions. In particular, it provides the significance, the flux measurement and the $95\%$ CL flux upper limit at 7 TeV for a point like source with a spectral index of $-2.7$, or for an extended source sizing 0.5*∘, 1.0∘* and 2.0*∘, with a spectral index of $-2.0$. We estimate the average flux at different sky directions for the detection at about $5\sigma$ significance to be $[8,9,10,20]\cdot 10^{-15}$ $({\rm TeV\,cm^{2}\,s})^{-1}$ for a point like, or extended source of size 0.5∘, 1.0∘* and 2.0*∘*, respectively. Most of the pulsars are predicted to have an ICS halo with an angular extension from a fraction of a degree to a few degrees (see Fig. 1), so we fix the flux sensitivity at 7 TeV to be $1\cdot 10^{-14}\;({\rm TeV\,cm^{2}\,s})^{-1}$, that is valid for a $1^{\circ}$ extended source. We note that we are not including any declination dependence of the sensitivity.

HESS. In order to estimate the flux sensitivity for HESS we use the information published in the HGPS catalog Abdalla et al. (2018a). The HESS Collaboration has calculated that the flux sensitivity for a point source with a spectral index of $-2.3$ is around $1\%$ of the Crab flux, i.e. about $3\cdot 10^{-13}\;({\rm TeV\,cm^{2}\,s})^{-1}$ at 1 TeV. This has been calculated for the optimistic case of an isolated source, while the real sensitivity is probably higher. They also show that there is a dependence of the flux sensitivity with the Galactic longitude. In the range between $l\in[40^{\circ},300^{\circ}]$ the sensitivity for point sources is roughly constant and at its lowest level, while it increases outside these directions. This result cannot be directly used for extended sources. We select therefore the sources detected as extended with $\theta_{68}\sim 0.1^{\circ}-0.4^{\circ}$ with the faintest flux at 1 TeV. We find a dependence for the flux of those sources with the size of extension. For sources detected with $\theta_{68}\approx 0.1^{\circ}$ the faintest detected flux is $5\cdot 10^{-13}$ $({\rm TeV\,cm^{2}\,s})^{-1}$, for $\theta_{68}\approx 0.2$ it is $1\cdot 10^{-12}$ $({\rm TeV\,cm^{2}\,s})^{-1}$ and for $\theta_{68}\approx 0.4$ it is $2\cdot 10^{-12}$ $({\rm TeV\,cm^{2}\,s})^{-1}$. We make the simplistic assumption of neglecting the dependence with $\theta_{\rm{68}}$, and fix the sensitivity to $1\cdot 10^{-12}$ $({\rm TeV\,cm^{2}\,s})^{-1}$ at 1 TeV. Moreover, we neglect the longitude dependence which is present in a minor portion of the HESS field of view.

CTA is the next generation ground-based observatory for $\gamma$-ray astronomy at VHE Mazin (2019). With more than 100 telescopes located in the northern and southern hemispheres, CTA will be the world’s largest and most sensitive high-energy $\gamma$-ray observatory. Ref. Ambrogi et al. (2018) ha calculated the flux sensitivity for the detection at the $5\sigma$ CL of an extended source with 50h observation time and different sizes of extension. The sensitivity flux at 10 TeV is $7\cdot 10^{-16}$ $({\rm TeV\,cm^{2}\,s})^{-1}$ ($1.2\cdot 10^{-15}$ $({\rm TeV\,cm^{2}\,s})^{-1}$) for a $0.1^{\circ}$ ($0.5^{\circ}$) extension. We will use $1\cdot 10^{-15}$ $({\rm TeV\,cm^{2}\,s})^{-1}$ in the rest of this section.

In Fig. 6 we show the number $N$ of ICS halos detectable by HAWC, HESS and CTA as a function of the efficiency $\eta$ (see Eq. 1). We calculate $\Phi_{\gamma}$ using Eq. 12 for all the ATNF catalog pulsars. If the flux is above the sensitivity of each experiment, it contributes to this number. The design of CTA is very promising for the detection of ICS halos. Indeed, with an efficiency a slow as a few %, this future experiment could detect about $100-130$ ICS halos. On the other hand, HAWC and HESS might have already detected around $25$ and $35$ halos, respectively. This is a realist number, given that 2HWC and HGPS catalogs contain 39 and 78 sources, and only a fraction of them are probably associated to ICS halos. We can revert the reasoning and use the number of sources detected in 2HWC and HGPS catalogs to find a rough upper limit for the average efficiency, which reads about 0.07 for 2HWC and 0.25 for HGPS.

IV.2 Ranking of the brightest expected ICS halos

We can also use our predictions for the ICS flux in order to outline the most promising targets among Galactic pulsars for the detection of a possible ICS halos. We pick the distance, age and spin down energy of pulsars from the ATNF catalog222http://www.atnf.csiro.au/research/pulsar/psrcat/ and calculate, using Eq. 12, the ICS flux ($\Phi_{\gamma}$) at 1 TeV, which is relevant for HESS, and at 10 TeV, where HAWC and the future CTA experiment have their peak of sensitivity. We rank the sources according to $\Phi_{\gamma}$ assuming that all the PWNe have the same efficiency $\eta=0.01$. We note that the efficiency acts as a mere normalization for the ICS flux, and does not influence the relative ranking of the sources. In Tab. 1 we report the list of the 23 highest pulsars in the HAWC field of view ranked according to the brightest predicted ICS halo flux at 10 TeV. We select only sources with $DEC\in[-20^{\circ},40^{\circ}]$ since this is the constrain of the HAWC field of view. We also report the predicted extension $\theta_{68}$ at 10 TeV calculated using Eq. 14 and for $D_{0}=7\cdot 10^{25}$ cm2/s. $\theta_{68}$ falls in the range between $0.40^{\circ}-0.80^{\circ}$ for most of the sources, while for Geminga and Monogem (2HWC J0700+143), which are very close sources, $\theta_{68}$ is about $7^{\circ}$ and $5^{\circ}$. This implies that $D_{0}$ should be of the order of $\sim 10^{26}$ cm2/s at 1 GeV if the $\gamma$-ray emission is due to ICS. Only two out of these 23 have not already been detected by HAWC. These two sources are associated to the pulsars PSR B1951+32 and PSR J1740+1000 and will very likely be reported in future HAWC catalogs. The 2HWC Survey reports for these sources a significance of 1.3$\sigma$ and 2.3$\sigma$, respectively. The fact that most of the sources in Tab. 1 have been already detected in 2HWC, demonstrates that the ICS flux is a very efficient indicator to select promising Galactic $\gamma$-ray sources.

A list of sources detectable (or already detected) by HAWC have been also presented in Linden et al. (2017). Indeed, some of the sources reported in this paper are also among the most promising ones in our list in Tab. 1. However, the complete list in Tab. 1 contains differences with respect to Linden et al. (2017). This is explained by the different estimation of the ICS flux. Indeed, Ref. Linden et al. (2017) uses a simplified model, that is based on a mere rescaling of the Geminga ICS flux observed by HAWC, defined through the distance and spin-down luminosity of the considered sources. The calculation is neglecting different ingredients which can vary the ICS flux, such as the source age. Moreover, the authors focus on sources with $T>100$ kyr. Instead, we perform for each source the complete calculation of their ICS flux, including also younger sources which can still exhibit an ICS halo.

HAWC is planning to operate at least until 2023 and to upgrade the detector and the data analysis (see, e.g., Joshi and Jardin-Blicq (2018)). These improvements and the increase of statistics will improve the sensitivity by a factor of at least 2. Since the results presented so far in the 2HWC catalog consider only 2 years of data, we can expect that it could be able to detect many more ICS halos in the near future.

According to the ICS flux at 10 TeV, we compile a list of pulsars promisingly detectable in the direction where HAWC could reasonably have the sensitivity to detect an ICS halo. We list these sources in Tab. 2, including the two non detected sources in Tab. 1. The $\theta_{68}$ and $\Phi_{\gamma}$ are computed as in Tab. 1. The fluxes are in the range between $10^{-16}-10^{-14}$ $({\rm TeV\,cm^{2}\,s})^{-1}$. As reported above, the HAWC sensitivity for the detection of an extended source is about $1\cdot 10^{-14}$ $({\rm TeV\,cm^{2}\,s})^{-1}$. With the future HAWC improvements, the first sources of Tab. 2 could be detected by HAWC. In case our efficiency, here fixed at $\eta=0.01$ would be underestimated, several other sources could be potentially detectable with HAWC.

V Derivation of $D_{0}$ in ICS halos

The main goal of our analysis is to estimate the diffusion coefficient $D_{0}$ around the pulsars under the hypothesis that the VHE $\gamma$-ray emission is due to ICS.

V.1 Selection of source sample

In this section we build a sample of sources in order to study the physical properties ($D_{0}$ and $r_{b}$) of ICS halo candidates. We focus on the detected emissions around pulsars at VHE since, as we have seen in the previous section, their angular extension is much smaller than at lower energies and makes the detection feasible for IACTs. Moreover, at these energies the pulsar proper motion does not effect significantly the ICS morphology.

We compute the ICS $\gamma$-ray flux for all the ATNF pulsars, and select the ones with the highest predicted ICS $\gamma$-ray flux, and having an extended counterpart already detected by HESS. Indeed, we will use the flux maps, which have been released in the HGPS catalog333https://www.mpi-hd.mpg.de/hfm/HESS/hgps/. We also add Geminga and Monogem for which the HAWC Collaboration has released the surface brightness Abeysekara et al. (2017a). We report in Tab. 3 the list of pulsars corresponding to these criteria with their age, distance and position in the sky. We also indicate the spatial extension as measured by HESS using a gaussian function ($\theta^{\rm{HESS}}_{\rm{gauss}}$).

We divide our sample in old and young pulsars fixing an age limit of 20 kyr. Indeed, as we described in Sec. II, $e^{\pm}$ are believed to be accelerated in PWNe to very high energies at the termination shock. This happens in the Sedov phase, i.e. in a time between a few up to twenty thousands of years after the supernova explosion (Chi et al., 1996; Blasi and Amato, 2011). After this stage, accelerated $e^{\pm}$ produce photons from radio, through synchrotron emission, up to VHE $\gamma$ rays by ICS. The size of extension thus depends on the PWN evolution. We consider separately the old and young PWN samples to inspect any dependence on the PWN evolution.

The list of sources in Tab. 3 exhibits an observed extended emission with $\theta^{\rm{HESS}}_{\rm{gauss}}$ $\sim[0.1^{\circ},0.5^{\circ}]$, which translates into a physical size of $\sim[8-35]$ pc. This size has been calculated by HESS using a spatial gaussian function, with the size of extension as the standard deviation parameter. We report also the predicted size of ICS emission calculated, for each source, with $\theta_{\rm{68}}$, i.e. as the 68% containment radius (see Eq. 14), at 1 TeV and for $D_{0}=7\cdot 10^{25}$ cm2/s. We apply the following procedure to calculate $\theta_{\rm{68}}$. We calculate the surface brightness for different angular distances from the source. Then we calculate, interpolating between the angle values considered, the distance that contains the 68% of the flux following the definition in Eq. 14. Overall, we find a good match between the measured and predicted size of extension, implying that the morphology of the $\gamma$-ray emission from these sources should be consistent with a diffusion environment with $D_{0}\sim 10^{26}$ cm2/s.

Most of the sources in our sample are located in the inner 4 kpc from the Earth and are younger than 100 kyr. 7 of them are classified in the TeVCat as unidentified, since no PWN structure has been identified in radio or X rays. However, a very powerful pulsar is found close to them, making the presence of a PWN a viable possibility.

We add here few comments about the association of few sources in Tab. 3. HESS J1858+020 is positionally compatible with the ATNF catalog pulsars PSR J1857+0143, J1857+0210 and B1855+02. However, assuming the same efficiency for all three, PSR J1857+0143 would have an ICS flux higher than a factor of 50 (100) with respect to J1857+0210 (B1855+02). For our purposes, we thus assume that HESS J1858+020 is associated to PSR J1857+0143. HESS J1303-631 position is compatible with PSR J1301-6305 and PSR J1301-6310. Computing ICS flux with the same efficiency for both, PSR J1301-6305 overclasses PSR J1301-6310 by a factor of about 50. Moreover, PSR J1301-6310 has a small distance from us and is relatively old, so the ICS flux is expected much more extended than $\theta^{\rm{HESS}}_{\rm{gauss}}=0.18^{\circ}$. Therefore, we associate HESS J1303-631 to PSR J1301-6305. Finally, HESS J1831-098 is found to have $TS=59$ in the main HGPS analysis, but only $TS=17$ in the cross-check analysis made using an alternative calibration, reconstruction, and gamma-hadron separation method, and is therefore considered as a source candidate Abdalla et al. (2018a).

V.2 Analysis technique

$D_{0}$ is derived through a fit to the surface brightness, source by source. We use the HESS $\gamma$-ray flux maps to derive the observed spatial distribution of the $\gamma$-ray emission around HESS sources. In order to extract the surface brightness as a function of the angular distance, we take the flux, the sensitivity and the significance maps from the HGPS catalog444https://www.mpi-hd.mpg.de/hfm/HESS/hgps/. These maps contain the flux integrated above 1 TeV, around a circular region defined by the correlation radius $R_{c}$. They are provided for $R_{c}=0.1^{\circ}$ and $0.2^{\circ}$, and with a pixel size of $0.02^{\circ}$. Therefore, each pixel contains information partially present also in the closest pixels. In order to limit this oversampling, we use the case with $R_{c}=0.1^{\circ}$ as our benchmark case, but we will comment on the results found with $R_{c}=0.2^{\circ}$. We show in Fig. 7 the flux maps for four sources in our sample.

We select a region of interest (ROI) around each source between $0.7^{\circ}$ and $1.1^{\circ}$ of radius depending on the extension of the source, as given in HGPS. We choose the size of the ROI in order to limit the contribution of background sources and include mainly the emission of the central sources. For example, HESS J1708-443 has a measured size of about $0.3^{\circ}$, so we choose a ROI of $0.8^{\circ}$ which contains entirely the $\gamma$-ray flux from that source. For sources extended $0.1^{\circ}$ ($0.4^{\circ}$) we typically select ROIs with $0.6^{\circ}$ ($1.1^{\circ}$) radius.

We do not include in our analysis a $\gamma$-ray background component, which could be modeled with the interstellar emission and flux from faint sources. Indeed, assuming that the background is isotropic, it should act as a mere normalization without changing significantly the angular profile of the TeV surface brightness. We also check whether there are other sources from the HGPS catalog in the ROI. If any other source is present in the same ROI, we remove the quadrant where this is located. For example, the source HESS J1616-508 is located at longitude and latitude $l_{S}=332.48^{\circ}$ and $b_{S}=-0.17^{\circ}$ and is close to HESS J1614-518 ($l=331.47^{\circ}$ and $b=-0.60^{\circ}$). Therefore, we remove from the analysis the region given by $l<l_{S}$ and $b<b_{S}$ in order to avoid any contamination from HESS J1614-518. We apply the same method to the following sources: HESS J1026-582, HESS J1303-631, HESS J1420-607, HESS J1458-608, HESS J1616-508, HESS J1632-478, HESS J1718-385, HESS J1825-137, HESS J1826-130, HESS J1831-098, HESS J1833-105, HESS J1841-055, HESS J1857+026, HESS J1858+020.

We assume for the ICS flux the one-zone diffusion model since the surface brightness data is not enough precise to constrain also the $r_{b}$ parameter. In order to use a two-zone diffusion model one would need to fit directly the HESS data around each source, but this is not available at the moment. Moreover, a two-zone diffusion model and the specific value of $r_{b}$ should affect mostly the $\gamma$-ray flux in the outer part of the halo, while we fit mainly the core of the $\gamma$-ray emission for each source.

This is the procedure we use to calculate the surface brightness of each source using the HESS flux maps. The maps are given as the $\gamma$-ray flux integrated over the solid angle (and have units of GeV/cm2/s/deg2). We calculate txhe total flux in concentric annuli and we divide it by their solid angle. We use as a benchmark case an annuli size bin of $0.08^{\circ}$. We also test larger and smaller values, finding very similar results.

Before using this technique to extract $D_{0}$ for each source, we have to verify if the flux maps extracted from the HGPS catalog represent well the flux distribution around the sources in our sample. In order to achieve this goal, we perform a fit to the surface brightness data assuming, as in the in the HGPS catalog, a gaussian function ($\propto\exp{(-\theta^{2}/(2\cdot\theta_{\rm{gauss}}^{2}))}$). Then, we compare the best fit values for the size of extension ($\theta_{\rm{gauss}}$) with the ones reported in the HGPS catalog ($\theta^{\rm{HESS}}_{\rm{gauss}}$). The best fit values and $1\sigma$ errors for $\theta^{\rm{HESS}}_{\rm{gauss}}$ and $\theta_{\rm{gauss}}$ are reported in Tab. 3 for $R_{c}=0.1^{\circ}$ . The source extensions we derive from the flux maps are compatible with the values reported in the HESS catalog. We find similar results using the flux maps provided for $R_{c}=0.2^{\circ}$. We are thus confident that, regardless the oversampling, the flux maps released by HESS can be used as a viable proxy to study the source spatial extension of the $\gamma$-ray flux.

VI Results for the diffusion around PWNe

We now perform a fit to the surface brightness data to find the diffusion coefficient around each PWN in our sample in Tab. 3. This is performed by using the ICS flux calculation (see Sec. II), by leaving $D_{0}$ and $\eta$ (see Eq. 1) as free parameters of the fit. The efficiency $\eta$ acts as a normalization, while the diffusion coefficient at 1 GeV $D_{0}$ modifies the angular profile of the ICS flux. We show in Tab. 4 and in Fig. 9 our results and the best fit and $1\sigma$ error for $D_{0}$. In Fig. 9, the diffusion coefficient (see Eq. 11) has been rescaled to 1 TeV, which is the typical energy scale of this analysis since we are considering VHE $\gamma$ rays. We also show in Tab. 4 the size of the ICS halo found implementing the empirical function Abeysekara et al. (2017b):

[TABLE]

where $\theta_{\rm{ICS}}/2$ is the angle that contains the $80\%$ of the observed flux. We find that this functional form indeed better describes for many sources in our sample the surface brightness data with respect to the gaussian function. In Fig. 8 we report the surface brightness data together with the best fit to the ICS model, found with $D_{0}$ as the free parameters. The best fit reproduces the observed surface brightness profile. Indeed, this model predicts the proper angular decrease of the surface brightness through the description of leptons diffusion around the source. We also show the fit with a mere gaussian template which, for these and several other sources, is a worst fit than the physical ICS model. This does not apply to all the sources in our sample but for most of them the ICS model is at least as good as the gaussian template.

The best fit values for $D_{0}$, source by source, are distributed in the $10^{25}-10^{27}$ cm2/s range. In particular, the mean value and the standard deviation over the entire sample are $D_{0}=9.1^{+17.4}_{-6.0}\cdot 10^{25}$ cm2/s. We find very similar values if we use a size for the annuli of $0.1^{\circ}$: $D_{0}=8.2^{+20.9}_{-5.9}\cdot 10^{25}$ cm2/s.

The results we find for the young PWN sample could be affected by the presence of the SNR, and could thus be systematically different from the ones found for the old sample, for which on the other hand the SNR has lost its power (see discussion in Sec. II). Therefore, we compute $D_{0}$ for the young and old PWN sample separately, finding $D_{0}=8.9^{+17.1}_{-5.9}\cdot 10^{25}$ cm2/s and $D_{0}=7.8^{+23.2}_{-5.8}\cdot 10^{25}$ cm2/s, respectively. There is no clear difference for the two PWN samples. Therefore, we do not see any evolution of $D_{0}$ with respect to the age as predicted by Evoli et al. (2018). This is visible in Fig. 9, where we show the value of the diffusion coefficient at 1 TeV ($D(1\,\rm{TeV})$) as a function of the PWN age. We remind that we show $D(1\,\rm{TeV})$ because we use $\gamma$-ray data above hundreds GeV that are produced for ICS by $e^{\pm}$ at TeV energies. We find for our sample $D(1\,\rm{TeV})=8.2^{+20.9}_{-5.9}\cdot 10^{26}$ cm2/s. There is a variation in the values of $D(1\,\rm{TeV})$ of about 1 order of magnitude. Our results for $D(1\,\rm{TeV})$ are compatible with the ones found for Geminga and Monogem with HAWC Abeysekara et al. (2017a) and Fermi-LAT data Di Mauro et al. (2019). We also show in Fig. 9 the results for the diffusion coefficient (scaled to 1 TeV considering the uncertainties on the normalization and the slope $\delta$) derived in Kappl et al. (2015); Genolini et al. (2019) from a fit to AMS-02 CR data. These numbers are representative of the average diffusion coefficient in the Galaxy. The intensity of $D(1\,\rm{TeV})$ we find with our analysis is about two orders of magnitude smaller than the one derived for the ISM. We also run our analysis on the HESS flux maps derived assuming $R_{c}=0.2^{\circ}$. We find a diffusion coefficient at 1 GeV for the entire sample of $D_{0}=13.6^{+33.1}_{-9.6}\cdot 10^{25}$ cm2/s while for the young and old PWN sample separately is $D_{0}=14.5^{+25.3}_{-9.2}\cdot 10^{25}$ cm2/s and $D_{0}=13.0^{+37.8}_{-9.7}\cdot 10^{25}$ cm2/s, respectively. These values are consistent within $1\sigma$ with the ones reported above for $R_{c}=0.1^{\circ}$.

An important characteristic of the ICS emission around PWNe is that their extension is connected to the size of the low-diffusion zone located around these sources (see Sec. III). In particular the size of the low-diffusion bubble must be at least large as the ICS region. We estimate the ICS halo size by considering the parameter $\theta_{\rm{ICS}}/2$ in Eq. 18. Then, we convert the angular scale into a physical size using $d\cdot\tan{(\theta_{\rm{ICS}}/2)}$. The average size of the ICS halo is $34^{+43}_{-19}$ pc for the entire sample, and $29^{+30}_{-15}$ pc and $38^{+52}_{-22}$ pc for the young and old sub-samples, respectively. We show in Fig. 10 the ICS halo size as a function of the age of the pulsar, together with the evolution model reported in Sec. II. In particular we use: $R\propto t^{1.2}$ for $t<3$ kyr, $R\propto t^{0.73}$ for $12<t<3$ kyr and $R\propto t^{0.3}$ for $t>12$ kyr. This model is compatible with the observed sizes and ages, and our results are comparable with the ones released for PWN by HESS Abdalla et al. (2018b). However, there is a large scatter in the data that prevents us to refine the model for the expansion rate as a function of time. The scatter we find is probably due to the fact that every pulsar has a different environment and a different evolution that makes the size of ICS flux significantly different for PWN with a similar age. Since the size of the ICS halos is of the order of 35 pc for the PWNe of our sample, this implies that the low-diffusion bubble should be at least large as this size. In particular for this average ICS halos size, $r_{b}$ should be at least of the order of 80 pc (see discussion in Sec. III and Fig. 4). However, some of the sources, e.g. HESS J1632-478, HESS J1825-137, HESS J1837-069, HESS J1841-055, HESS J1912+101 and HESS J1303-631, have a much more extended ICS halos size implying that the size of the low-diffusion bubble could even exceed 100 pc.

These results have been obtained within the one-zone diffusion model (see Sect. II). We now explore the possibility that a low diffusive regime is present in a region around the source within a radius $r_{b}$, according to Eq. 11. In order to show the effect of $r_{b}$ on the surface brightness we consider the very bright HESS J1825-137 source, which has surface brightness data with relatively small uncertainties (see Fig. 8). We calculate the best fit for $D_{0}$ for a two-zone diffusion model with $r_{b}$ variable between $10-200$ pc. The range of $D_{0}$ and $r_{b}$ that best represents the data are: $D_{0}\in[2.5,15]\cdot 10^{26}$ cm2/s and $r_{b}>60$ pc. The best fit is for $D_{0}=6\cdot 10^{26}$ cm2/s and $r_{b}=80$ pc, but the $\chi^{2}$ distribution is flat for $r_{b}>70$ pc and gives comparably good fit for increasing values of $r_{b}$ and decreasing values of $D_{0}$. In particular, for $r_{b}=[60,80,100,120]$ pc the best for $D_{0}$ is $D_{0}\in[15.9,6.3,2.5,2.2]\cdot 10^{26}$ cm2/s. Therefore, for $r_{b}>80$ pc the best fit for $D_{0}$ tends to the value we find with the one-zone diffusion model (see Tab. 4). We show in Fig. 11 the contour plot for the $\chi^{2}$ values as a function of $D_{0}$ and $r_{b}$. This exercise demonstrates that surface brightness data could be used in principle to bound the size of the low-diffusion bubble. However, it is prohibitive to run this analysis for all the sources in our sample, because the surface brightness data for most of the sources have large uncertainties.

In Ref. Giacinti et al. (2019), a sample of PWNe and PWN candidates from the HGPS catalog have been considered to estimate the density of $e^{\pm}$ in the ICS halo. The $e^{\pm}$ density has been calculated with different methods finding that, for most of the sources, it is larger than the one of the ISM. This implies that the $e^{\pm}$ that produce the ICS halos are probably confined in the PWN. One of the main assumption in that paper is the size of the ICS region, which is taken directly as the source extension provided in the HGPS catalog, i.e. as the standard deviation for a Gaussian spatial distribution of $\gamma$ rays. These sizes are probably an underestimate of the halo size. Indeed for many sources the sizes they assume are much smaller than the values we report in Tab. 4 with $\theta_{\rm{ICS}}/2$. In particular, this happens for the following sources: HESS J1718-385, HESS J1809-193, HESS J1813-178, HESS J1825-137, HESS J1858+020, HESS J1908+063, HESS J1303-631, HESS J1356-645, HESS J1420-607 and HESS J1833-105. The differences in the halo size is for most of the sources about a factor of 2 thus bringing a difference in halo volume of almost 1 order of magnitude. If this factor is considered in their calculation, many of their sources would have a $e^{\pm}$ density comparable to the one of the ISM. This would change significantly their conclusion because they could not exclude that, for most of their sources, the $e^{\pm}$ are probably not confined in the PWN and actually are traveling in the ISM. For example the source HESS J1825-137 has a $e^{\pm}$ density of about 0.2 eV/cm3 in Giacinti et al. (2019), about twice the one of the ISM, calculated using a size of the halo of 48.3 pc. On the other hand, we find for the same source that the size is about 73 pc. Using this number, the $e^{\pm}$ density becomes 0.06 eV/cm3, i.e. smaller than the ISM one.

We can now estimate the efficiency $\eta$ and the power-law index $\gamma_{e}$ using the measurement of the differential flux at 1 TeV and of the $\gamma$-ray flux spectral index published in the HGPS catalog. Specifically, we fit the value of $\gamma_{e}$ to the observed $\gamma$-ray slope and then we find the efficiency which reproduces the flux data at 1 TeV. The $\gamma_{e}$ values are derived assuming Eq. 1 for the $e^{\pm}$ injection spectrum. We report in Tab. 4 the results for $\gamma_{e}$, together with the size of the halo and the diffusion coefficient. Indeed, $\gamma_{e}$ turns out to be well determined for each single source, but showing a great spread over the whole sample. Most of the values of $\gamma_{e}$ are very soft and in the range $2.5-3.0$. Only 7 of them are harder, with values between $1.2-1.9$. Given that our study is devoted to energies well above the TeV, we do not introduce any further break at lower energies as instead assumed to model low-energy data from Fermi-LAT and X-ray telescopes Torres et al. (2014).

We calculate the efficiency for the conversion of pulsar spin-down energy into $e^{\pm}$ using Eq. 1. We assume for each source the $\gamma_{e}$ values reported in Tab. 4. We find very high values of $\eta$, sometimes even exceeding 1. These high values for $\eta$ are likely due to the choice not to set a break into Eq. 1. The $e^{\pm}$ injection spectrum is indeed usually modeled with a broken power law with a break around hundreds GeV and with an index above and below the break of about 1.4 and 2.2, respectively (see, e.g., Torres et al. (2014).) The bias here is that we are extrapolating very soft indexes below the energy break where actually $\gamma_{e}$ hardens. Indeed, we remind that the efficiency is calculated from an $e^{\pm}$ energy of 0.1 GeV while this analysis is constraining the injection spectrum for TeV energies. We can revert the sentence inferring that the $e^{\pm}$ injection spectrum is probably harder below the energy range covered by HESS. For all the reasons reported above we decide to not show the values of $\eta$ that we have found. In order to constrain more efficiently the efficiency and $e^{\pm}$ injection spectrum $\gamma$-ray data at GeV energies must be considered. Fermi-LAT data are ideal to this scope, as we have already shown in Di Mauro et al. (2019). We are planning to perform, in a followup paper, a combined analysis of HESS and Fermi-LAT data from the sources considered in this paper in order to derive $\eta$ and $\gamma_{e}$.

VII Conclusions

The detection of low-diffusion regions, few tens of pc wide, found around Geminga and Monogem pulsars analyzing Fermi-LAT Di Mauro et al. (2019) and HAWC Abeysekara et al. (2017a) $\gamma$-ray data raises the question if this is a peculiarity or a general property of Galactic pulsars.

In this paper, trying to answer this question, we have studied the physical properties of these halos, believed to be generated by $e^{\pm}$ accelerated by PWN and ICS with the ISRF. We have studied the size of ICS halos as a function of the strength and size of the low-diffusion bubble, the age and distance of the host pulsar, and of its proper motion. We find that current IACTs are able to probe diffusion coefficients $\leq 10^{27}$ cm2/s around most of the pulsars closer than 10 kpc from the Earth. We show that, at VHE, the pulsar proper motion has a limited effect on the ICS spatial morphology.

We then rank ATNF pulsars according to the ICS flux and demonstrate that this parameter is very efficient to indicate the detectability of each source. Indeed, out of 23 pulsars in the HAWC field of view and predicted by our model to have the brightest ICS halo fluxes, 21 have been included in the 2HWC catalog. We provide in Tab. 2 the list of sources not yet detected by HAWC, and ranked by their ICS $\gamma$-ray flux. Given the ICS emission is the process producing the VHE photons similarly with Geminga and Monogem, we predict these sources to be the next-to-be-discovered as ICS halos in HAWC data. As a further prediction, we also list the angular size of the ICS halo of each source. The number of ICS halos potentially already detected by HAWC and HESS ranges between 25-50 assuming a conversion efficiency $\eta$ at the % level. As for CTA, an efficiency as low as 0.01 could lead to the detection of at least one hundred ICS halos.

We employ the flux maps provided in the HGPS catalog and the Geminga and Monogem surface brightness published by the HAWC Collaboration in order to derive the diffusion around a sample of 27 PWNe and PWN candidates. We demonstrate that the $e^{\pm}$, released from the sources in our sample, propagate in a low-diffusion Galactic environment with a diffusion coefficient which is about two orders of magnitude lower than the value recently derived for the entire Galaxy through a fit to AMS-02 CR data. The mean value and the standard deviation over the entire sample are at 1 GeV $D_{0}=9.1^{+17.4}_{-6.0}\cdot 10^{25}$ cm2/s. We do not register any dependence of this numbers on the age of pulsar, meaning that probably the effect of confinement of the SNR and PWN is not very strong even for the younger sources in our sample. The characterization of the pulsar environment by a low diffusion region turns out to be a general trend for all the analyzed sources. The size of the ICS halos have been found to be on average 35 pc implying that the low-diffusion bubbles should be larger than this size. For some of the sources in our sample, e.g. HESS J1632-478, HESS J1825-137, HESS J1837-069, HESS J1841-055, HESS J1912+101 and HESS J1303-631, the low-diffusion bubble size could exceed 100 pc. These numbers should be used as an estimate for $r_{b}$ in the two-zone diffusion model employed to propagate $e^{\pm}$ from the pulsar to the Earth. Since, as we have explained in Sec. V.1, the PWNe consider in this paper are also the highest ranked according to the ICS flux at TeV energies, we do not believe our results are biased towards objects that have smaller $D_{0}$ and so have a more concentrated $\gamma$-ray emission. The consequences of the present results for the interpretation of the $e^{+}$ flux data in terms of Galactic PWNe and for the propagation of cosmic rays will be investigating in a forthcoming paper.

Acknowledgements.

The authors thank Luigi Tibaldo, Ke Fang, Andrew James Smith, Regina Caputo and Roger Romani for insightful discussions. MDM acknowledges support by the NASA Fermi Guest Investigator Program 2019 Cycle 12 through the Fermi Program N. 121119 (P.I. MDM) entitled “Detecting $\gamma$-ray halos around PWNe and interpretation of the positron excess”. The work of FD and SM is supported by the ”Departments of Excellence 2018 - 2022” Grant awarded by the Italian Ministry of Education, University and Research (MIUR) (L. 232/2016). FD and SM acknowledge financial contribution from the agreement ASI-INAF n.2017-14-H.0 and the Fondazione CRT for the grant 2017/58675.