The contribution of Galactic TeV pulsar wind nebulae to Fermi large area telescope diffuse emission
Pulsar wind nebulae are expected to contribute to γ observations both in the GeV and TeV energy domains. We indicate with ΦGeV (ΦTeV) the integrated source flux in the energy range 1–100 GeV (1–100 TeV) probed by Fermi-LAT (H.E.S.S.). We assume that all the sources in the considered population have approximately the same emission spectrum, described by a broken power-law with different spectral indexes βGeV and βTeV in the GeV and TeV energy domain and with a transition energy E0 = [0.1–1.0] TeV located between the ranges probed by Fermi-LAT and H.E.S.S. At high energies (E ≥ E0), we allow the source spectral index to move inside the range βTeV = [1.9–2.5] measured by H.E.S.S.22 for identified PWNe, see Supplementary Table 1. The index βGeV is instead determined by requiring realistic values for the parameter RΦ, defined as the ratio
$${R}_{{{\Phi }}}\equiv \frac{{{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}}{{{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}}$$
(1)
between fluxes emitted by a given source in different energy domains. As it is discussed in the “Methods” section, we obtain a consistent description of the HGPS and the Fermi-LAT Fourth Source Catalog (4FGL-DR2) for \({R}_{{{\Phi }}}=\left[250{-}1500\right]\) that corresponds to βGeV inside the global range βGeV = [1.06–2.19] (see Eq. (8) for the general relationship between the spectral parameters in our analysis). The assumed source spectrum can be further validated by considering the average observational properties of PWNe observed by Fermi-LAT and H.E.S.S. in the GeV/TeV domain, see the “Methods” section for details. Moreover, the corresponding spectral shapes are consistent with theoretical predictions for γ-ray emission from PWNe23,24.
The PWNe population in the TeV domain
The properties of the considered source population can be constrained by observation in the TeV energy domain. Following a previous work20, PWNe distribution is described by:
$$\frac{dN}{{d}^{3}r\,d{L}_{{{{{{{{\rm{TeV}}}}}}}}}}=\rho \left({{{{{{{\bf{r}}}}}}}}\right){Y}_{{{{{{{{\rm{TeV}}}}}}}}}\left({L}_{{{{{{{{\rm{TeV}}}}}}}}}\right)$$
(2)
where r indicates the source distance from the Galactic Center. The function ρ(r) describes the spatial distribution of the sources and it is conventionally normalized to one when integrated in the entire Galaxy. It is assumed to be proportional to the pulsar distribution in the Galactic plane25. The source density along the direction perpendicular to the Galactic plane is assumed to scale as \(\exp \left(-\left|z\right|/H\right)\) where H = 0.2 kpc represents the thickness of the Galactic disk.
The function YTeV(LTeV) gives the source intrinsic luminosity distribution in the TeV energy domain. It is parameterized as a power-law:
$${Y}_{{{{{{{{\rm{TeV}}}}}}}}}({L}_{{{{{{{{\rm{TeV}}}}}}}}})=\frac{R\,\tau \,(\alpha -1)}{{L}_{{{{{{{{\rm{TeV}}}}}}}},\max }}{\left(\frac{{L}_{{{{{{{{\rm{TeV}}}}}}}}}}{{L}_{{{{{{{{\rm{TeV}}}}}}}},\max }}\right)}^{-\alpha}$$
(3)
that extends in the luminosity range \({L}_{{{{{{{{\rm{TeV}}}}}}}},\min }\,\le \,{L}_{{{{{{{{\rm{TeV}}}}}}}}}\,\le \,{L}_{{{{{{{{\rm{TeV}}}}}}}},\max }\)26. This functional form, that is generically adopted in population studies, is naturally obtained for a population of fading sources, such as PWNe or TeV Halos, produced with a constant rate R and having intrinsic luminosity that decreases over a time scale τ, see the “Methods” section for details.
Previous analyses on the subject2,4 have been performed under the assumption that the index of the luminosity distribution is α = 1.8 because this leads to a good description of observational data. We conform to this choice for our reference case and we note that the value α = 1.8 is also obtained for a population of pulsar-powered sources, if the efficiency of TeV emission is correlated to spin-down power as suggested by H.E.S.S. data16. However, to show the dependence of our results on the performed assumptions, we also consider the alternative hypothesis α = 1.5 that is obtained by postulating that the efficiency of TeV emission is constant in time.
By fitting the flux, latitude and longitude distribution of bright sources in the HGPS catalog and assuming that the PWNe birth rate is equal to that of core-collapse SN explosions, i.e., R = 0.019 yr−1, one obtains \({L}_{{{{{{{{\rm{TeV}}}}}}}},\max }=6.8\times 1{0}^{35}\,{{{{{{{\rm{erg}}}}}}}}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}\) (\({L}_{{{{{{{{\rm{TeV}}}}}}}},\max }=4.9\times 1{0}^{35}\,{{{{{{{\rm{erg}}}}}}}}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}\)) and τ = 0.5 × 103 y (τ = 1.8 × 103 y) for α = 1.8 (α = 1.5)20. Note that our analysis is only sensitive to the product Rτ, indeed a smaller PWN formation rate can be balanced by a higher value of τ (and viceversa) with no consequences for the present discussion. The fit to HGPS sources permits us to constrain the cumulative flux \({{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}^{{{{{{{{\rm{tot}}}}}}}}}\) produced by PWNe population in the TeV domain with~ 30% statistical accuracy, being it equal to \({{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}^{{{{{{{{\rm{tot}}}}}}}}}=\left(5.{9}_{-1.5}^{+1.8}\right)\times 1{0}^{-10}\,{{{{{{{{\rm{cm}}}}}}}}}^{-2}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}\) for α = 1.8. To see a comparison of the observed cumulative flux with the one obtained in our model see Supplementary Note 2. This estimate does not critically depend on the adopted assumptions (e.g., the source space distribution, the Galactic disk thickness, the source physical dimensions, etc.). The largest effect is obtained by modifying the luminosity index α. We get e.g., \({{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}^{{{{{{{{\rm{tot}}}}}}}}}=\left(3.{8}_{-1.1}^{+1.2}\right)\times 1{0}^{-10}\,{{{{{{{{\rm{cm}}}}}}}}}^{-2}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}\) for α = 1.5 that corresponds to ~35% reduction with respect to the reference case α = 1.820.
The above results have been obtained by assuming βTeV = 2.3 which corresponds to the average spectral index of sources observed by H.E.S.S.. It should be remarked, however, that the adopted value of βTeV has no effects on the cumulative flux \({{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}^{{{{{{{{\rm{tot}}}}}}}}}\), neither on the distribution dN/dΦTeV of sources as function of flux in the TeV domain. These quantities are thus directly determined by observational data, independently on assumptions for the source emission spectrum.
The total and unresolved emission in the GeV domain
The total flux produced at Earth by TeV PWNe population in the GeV domain depends on the parameter RΦ and it is given by:
$${{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{tot}}}}}}}}}={R}_{{{\Phi }}}\,{{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}^{{{{{{{{\rm{tot}}}}}}}}}.$$
(4)
The parameter RΦ also determines the distribution of sources as a function of the flux they emit at GeV energies, according to:
$$\frac{dN}{d{{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}}=\frac{1}{{R}_{{{\Phi }}}}\frac{dN}{d{{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}}\left({{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}/{R}_{{{\Phi }}}\right)\,.$$
(5)
Faint sources cannot be individually resolved by Fermi-LAT and contribute to the large-scale diffuse emission observed by this experiment. The unresolved contribution can be calculated as:
$${{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{NR}}}}}}}}}=\int\nolimits_{0}^{{{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{th}}}}}}}}}}d{{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}\,{{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}\,\frac{dN}{d{{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}}$$
(6)
where \({{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{th}}}}}}}}}\) is the Fermi-LAT detection threshold. For objects contained in the Galactic plane, this is estimated as \({{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{th}}}}}}}}}=1{0}^{-9}\,{{{{{{{{\rm{cm}}}}}}}}}^{-2}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}\) by looking at the turnover of the observed source number as a function of the photon flux above 1 GeV (see Fig. 24, panel (a) of the recent work27 by Acero et al.). By considering that the flux distribution scales as \(dN/d{{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}\propto {{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}^{-\alpha }\) for ΦTeV → 0 (see the “Methods” section), we expect that \({{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{NR}}}}}}}}}\propto {R}_{{{\Phi }}}^{\alpha -1}\). We remark that the total and unresolved fluxes, \({{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{tot}}}}}}}}}\) and \({{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{NR}}}}}}}}}\), only depend on RΦ and are independent on the assumed values for the spectral parameters E0 and βTeV.
In the last line of Table 1, we give the flux \({{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{NR}}}}}}}}}\) produced by PWNe that are not resolved by Fermi-LAT for the two extreme values RΦ = 250 and 1500. These fluxes are compared with the large-scale diffuse emission associated with interstellar gas \({{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{diff}}}}}}}}}\) detected by Fermi-LAT (see second column in Table 1) in the 1–100 GeV energy range and determined in Pothast et al.2 by using 9.3 years of Fermi-LAT Pass 8 data. The energy integrated fluxes have been obtained by interpolating the experimental points and integrating in the energy range 1–100 GeV. We see that unresolved emission by PWNe corresponds to a fraction ~3% (for RΦ = 250) and ~11% (for RΦ = 1500) of the diffuse gamma-ray emission associated with interstellar gas. The above results are obtained by assuming that the source luminosity distribution index is α = 1.8 to conform with previous analyses on the subject2,4 that have been performed under this hypothesis. Results for α = 1.5 are smaller and are reported in the last two columns of Table 1.
In order to probe the radial dependence of the PWNe contribution, we repeat our calculations by considering the Galactocentric rings adopted in Pothast et al.2. The flux produced by unresolved TeV PWNe in each ring is compared with the Fermi-LAT diffuse emission from the same region. As we see from Table 1, the unresolved contribution becomes more relevant in the central rings, due to the fact that the source density (and the average distance from the Sun position) is larger. In the most internal region (1.7 ≤ r ≤ 4.5 kpc), unresolved sources account for ~9% (~36%) of the Fermi-LAT diffuse emission associated with interstellar gas for RΦ = 250 (RΦ = 1500) and α = 1.8. We do not consider the central region r ≤ 1.7 because it is affected by large systematic errors2. This clearly shows that this component is not negligible and cannot be ignored in the interpretation of Fermi-LAT diffuse emission data.
Spectral analysis
The effect of the unresolved TeV PWNe population on the determination of CR diffuse emission spectral index is displayed in Fig. 1. The purpose of this figure is not to discuss comprehensively the effects of parameter variations in our calculation. Rather, our goal is to illustrate our approach and to explain why, despite the extremely large range of variation of the RΦ parameter (determining the PWNe integrated flux in the GeV domain), one still gets a prediction for the spectral index of CR diffuse emission. For this reason, we fix the spectral parameters to the values that better reproduce the cumulative spectral energy distribution of the PWNe observed both by Fermi-LAT and H.E.S.S. (i.e., βTeV = 2.4 and E0 = 0.8 TeV, see Fig. 2) and we only vary the flux ratio in the range RΦ = 250−1500. On the other hand, the final results of our analysis reported in Table 2 and displayed in Fig. 3, also take into account the effects of possible variations of E0 and βTeV.
Black data points show the total diffuse γ-ray emission associated with interstellar gas measured by Fermi-LAT in each galactocentric ring, the error bars represent the statistical error see Pothast et al.2 for details. The red bands represent the predicted contribution of unresolved TeV pulsar wind nebulae for α = 1.8, E0 = 0.8 TeV, βTeV = 2.4 and RΦ ranging between 250 and 1500. Green lines show the diffuse cosmic ray emission inferred by fitting the data with (solid) and without (dashed) including the pulsar wind nebulae contribution. The dark green bands show the systematic error produced by variations of the flux ratio RΦ. Light green bands show the total systematical uncertainty obtained when E0 and βTeV are also allowed to vary. Blue lines represent the total gamma fluxes predicted as a function of the energy for α = 1.8, E0 = 0.8 TeV, βTeV = 2.4 and RΦ ranging between 250 and 1500. The gray lines show a power-law with an index of 2.7 for comparison.
The cumulative spectral energy distribution of the Pulsar Wind Nebulae observed both by Fermi-LAT and H.E.S.S. (black thick line). We show with a red line a broken power-law spectrum with an energy break E0 = 0.8 TeV and spectral indexes βTeV = 2.4 and βGeV = 1.89. The shaded bands are obtained by propagating statistical and systematic uncertainties in the source spectra given in 4FGL-DR236 (orange band) and HGPS22 (blue band) catalogs. The dashed black line represents the average theoretical spectrum obtained by integrating over the spectral parameter space (see Sect. Robustness of our results and comparison with previous studies). It corresponds to the central values of ΓBF given in Table 2 and it is normalized in order to produce the same number of photons in the TeV energy domain as the observed cumulative Pulsar Wind Nebulae spectrum.
The gamma-ray emissivity (a) and the cosmic ray spectral index (b) obtained in this work (black points) compared with the ones in Peron et al.29 (gray points), Pothast et al.2 and in Acero et al.4 (orange points). The error bars for black points are obtained by summing in quadrature statistical and systematical uncertainties, see Table 2. In particular, thin error bars show the systematic uncertainties conservatively estimated by varying the spectral parameters in our phenomenological model, while thick error bars only include statistical uncertainties in the fitting procedure. The blue dashed line represents the CR distribution predicted by the GALPROP code30 normalized at 8.5 kpc at the emissivity value obtained in this work in the ring 8–10 kpc.
Black data points in Fig. 1 represent the total γ − ray flux associated with interstellar gas observed by Fermi-LAT in each galactocentric ring in 25 log-spaced energy bins between 0.34−228.65 GeV and in the latitude window ∣b∣ < 20.25∘. These data have been previously fitted in Pothast et al.2 with a single power-law \(\propto {E}^{-{{{\Gamma }}}_{1}}\), obtaining the green dashed lines reported in Fig. 1. The decrease of the best-fit spectral indexes Γ1 in the inner rings with respect to the locally observed value, see second column of Table 2, has been considered as the evidence of a progressive large-scale hardening of CRs spectrum toward the Galactic Center. The same conclusion was obtained by previous analyses on the subject3,4 performed by using a similar approach. One can get a visual perception of the situation by comparing the green dashed lines with the gray solid lines in Fig. 1 that describe power laws with spectral index fixed at the local value, i.e., ~2.7, suitably normalized to reproduce the observed flux at 2 GeV.
The above conclusion is only valid if the unresolved source contribution is negligible, so that the total observed emission can be identified with the “truly” diffuse component produced by CRs interaction with interstellar matter. This assumption is, however, not adequate in the inner Galaxy, as it is shown with red solid lines in Fig. 1 that give the unresolved PWNe contribution as a function of energy for the reference case α = 1.8 and the two extreme values RΦ = 250 and 1500. The red shaded area can be considered as the systematical uncertainty associated to the parameter RΦ. The effects of possible variations of E0 and βTeV on the unresolved PWNe emission are shown in Supplementary Note 3, see Supplementary Fig. 2. The relevant point to note in this figure is that the GeV source spectral index βGeV and the flux ratio RΦ are correlated, as it is discussed in the “Methods” section (see Eq. (8)). As a result of this, PWNe unresolved emission in the inner Galaxy is either relatively large or has a hard spectrum, providing a contribution at E ~ 100 GeV that is almost independent on RΦ. This is the natural consequence of the fact that the source emission above 1 TeV is observationally fixed by HGPS data. This important piece of information cannot be neglected and it is included in our work.
If we take unresolved PWNe emission into account, the evidence for CR spectral hardening in the inner Galaxy may be considerably weakened, as it is shown by the green thick solid lines in Fig. 1 that represent the component of the total diffuse gamma-ray flux that can be ascribed to CR interactions. This is still parameterized as a single power-law \(\propto {E}^{-{{{\Gamma }}}_{{{{{{{{\rm{BF}}}}}}}}}}\) (the number of degrees of freedom in the fit is not changed) but the total flux, described by blue lines in Fig. 1, is obtained as the sum of CR diffuse emission plus the unresolved PWNe contribution. The best fit spectral indexes ΓBF of the truly diffuse emission obtained in each ring are larger and closer to the value measured at the Sun position with respect to those obtained in the previous analyses2,3,4 that do not take unresolved sources into account, see Table 2. Our results are mildly dependent on the flux ratio RΦ, as it can be understood by looking at the dark green bands in Fig. 1 that show the effects of varying this parameter in the range RΦ = [250–1500]. The light green bands also take into account possible variations of the spectral parameters E0 and βTeV and provide a conservative estimate of the total systematical error for CR diffuse emission, as it discussed in the next paragraph.
Robustness of our results and comparison with previous studies
In order to estimate the uncertainties in our approach, we repeat our analysis for different combinations of the spectral parameters (RΦ, E0, βTeV). The final results of our analysis are given in Table 2. Here, the first errors describe systematic uncertainties, evaluated as the maximal variations of ΓBF that are obtained when (RΦ, E0, βTeV) are simultaneously varied in the 3-dim parameter space defined by the ranges \({R}_{{{\Phi }}}=\left[250{-}1500\right]\), \({E}_{0}=\left[0.1{-}1.0\right]\,{{{{{{{\rm{TeV}}}}}}}}\) and \({\beta }_{{{{{{{{\rm{TeV}}}}}}}}}=\left[1.9{-}2.5\right]\). We see that our conclusions are stable and not challenged by possible systematic effects connected with the assumed source spectrum. It is important to remark that our estimates for systematic uncertainties are very conservative. The smaller (larger) values for ΓBF are e.g., obtained by assuming that all sources in the considered population have βTeV = 1.9 (βTeV = 2.5), E0 = 1.0 TeV (E0 = 0.1 TeV) with a marginal dependence on the assumed RΦ, i.e., they correspond to a physical situation that is extremely unlikely. TeV PWNe in our Galaxy are indeed expected to have a distribution of spectral properties with compensating effects among extreme assumptions. The central values for ΓBF given in Table 2 are obtained by integrating over the whole parameters space. We assume logarithmic uniform distributions for the spectral break position and for the flux ratio, while for βTeV we consider a Gaussian distribution centered in βTeV = 2.4 and with dispersion 0.15 as reported in the HGPS catalog22.
In addition to the reference case α = 1.8 (third column) that is obtained by using the luminosity distribution index considered by previous analyses2,4, we also display the results obtained by assuming the alternative value α = 1.5 (last column). In this case, one obtains smaller effects on ΓBF, coherently with the fact that unresolved PWNe emission in the GeV domain is smaller, see Table 1. It would be important to have further phenomenological and/or theoretical constraints on the α parameter for future analyses.
The results of our reference case (α = 1.8) are compared with those given by other analyses in Fig. 3, where we show the γ-ray emissivity per H atom at 2 GeV (a) which is a proxy of the CR spatial distribution in the Galaxy, and the CR proton spectral index (b), obtained by adding 0.1 to the spectral indexes of the truly diffuse gamma emission28. Black points show the results of our work that are compared to those given by Pothast et al.2 (red points) and Acero et al.4 (orange points). We also show with gray points the results obtained by Peron et al.29 by studying γ-ray emission in the direction of giant molecular clouds. The thin error bars (for the black points) show the systematic uncertainties conservatively estimated as discussed above while the thick error bars only include statistical uncertainties. We see that the emissivity calculated in this work is in good agreement with that obtained by Pothast et al.2 and Peron et al.29. This is not surprising because we don’t expect any significant effect at 2 GeV due to the presence of unresolved sources. The three data sets agree quite well with theoretical expectations for the CR distribution from GALPROP code30 (dashed blue line). The theoretical CR distribution is shown e.g., in Fig. 8 of Acero et al.4 where the specific GALPROP configuration is also given. The inclusion of unresolved PWNe affects the CR spectral index that can be increased up to 0.18 in the central ring adjusting it to the locally observed value, i.e., ~2.8. The cosmic ray reconstructed spectrum still shows a residual difference with the local value in the other rings. We see, however, that unresolved PWNe naturally accounts for a large part of the spectral index variation as a function of r that has been reported by previous analyses, weakening considerably the evidence for CR spectral hardening in the inner Galaxy.
Conclusions
The TeV Galactic sky is dominated by a population of bright young PWNe whose properties are constrained by present H.E.S.S. Galactic Plane Survey (HGPS) data. We predict the cumulative emission produced by this population in the GeV domain within a phenomenological model that is based on the average spectral properties of identified PWNe. Our phenomenological description is also consistent with theoretical predictions of PWNe emission spectrum23,24 and it can be refined in the future by adopting the recent results by Fiori et al.31, appeared during the reviewing procedure of this work. We argue that a relevant fraction of the TeV PWNe population cannot be resolved by Fermi-LAT. The γ-ray flux due to unresolved TeV PWNe and the truly diffuse emission, due to CR interactions with the interstellar gas, add up contributing to shape the radial and spectral behavior of the total diffuse γ-ray emission observed by Fermi-LAT.
The spatial distribution of TeV PWNe, peaking around r = 4 kpc from the Galactic Center, combined with the detector flux threshold modulates the relative contribution of unresolved sources in different Galactocentric rings. In particular, the relevance of this component increases in the inner rings where the total diffuse emission has a different spectral distribution with respect to the local one. Previous analyses neglected the contribution due to unresolved PWNe and interpreted the observed spectral behavior of the total diffuse emission as indirect evidence for CR spectral hardening toward the Galactic center2,3,4. We have shown that the emergence of PWNe unresolved component in the central region, which is characterized by an average spectral index βGeV < 2, can affect this conclusion, by naturally accounting for (a large part of) the spectral index observed variation as a function of r. Our results could also solve the tension, discussed in Cataldo et al.32, between total γ-ray emission measured by H.E.S.S., Milagro, Argo and HAWC and that obtained by implementing CR spectral hardening.