Blog

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 E₀ = [0.1–1.0] TeV located between the ranges probed by Fermi-LAT and H.E.S.S. At high energies (E ≥ E₀), we allow the source spectral index to move inside the range β_TeV = [1.9–2.5] measured by H.E.S.S.²² 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

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 PWNe^23,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 work²⁰, PWNe distribution is described by:

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 plane²⁵. 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 Y_TeV(L_TeV) gives the source intrinsic luminosity distribution in the TeV energy domain. It is parameterized as a power-law:

that extends in the luminosity range \({L}_{{{{{{{{\rm{TeV}}}}}}}},\min }\,\le \,{L}_{{{{{{{{\rm{TeV}}}}}}}}}\,\le \,{L}_{{{{{{{{\rm{TeV}}}}}}}},\max }\)²⁶. 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 subject^2,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. data¹⁶. 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⁻¹, 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 × 10³ y (τ = 1.8 × 10³ y) for α = 1.8 (α = 1.5)²⁰. 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.8²⁰.

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:

The parameter R_Φ also determines the distribution of sources as a function of the flux they emit at GeV energies, according to:

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:

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 work²⁷ 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 E₀ 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.² 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 subject^2,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.². 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 errors². 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 E₀ = 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 E₀ 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.² for details. The red bands represent the predicted contribution of unresolved TeV pulsar wind nebulae for α = 1.8, E₀ = 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 E₀ and β_TeV are also allowed to vary. Blue lines represent the total gamma fluxes predicted as a function of the energy for α = 1.8, E₀ = 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 E₀ = 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-DR2³⁶ (orange band) and HGPS²² (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.²⁹ (gray points), Pothast et al.² and in Acero et al.⁴ (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 code³⁰ 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.² 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 Γ₁ 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 subject^3,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 E₀ 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 analyses^2,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 E₀ 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_Φ, E₀, β_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_Φ, E₀, β_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), E₀ = 1.0 TeV (E₀ = 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 catalog²².

In addition to the reference case α = 1.8 (third column) that is obtained by using the luminosity distribution index considered by previous analyses^2,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 emission²⁸. Black points show the results of our work that are compared to those given by Pothast et al.² (red points) and Acero et al.⁴ (orange points). We also show with gray points the results obtained by Peron et al.²⁹ 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.² and Peron et al.²⁹. 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 code³⁰ (dashed blue line). The theoretical CR distribution is shown e.g., in Fig. 8 of Acero et al.⁴ 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 spectrum^23,24 and it can be refined in the future by adopting the recent results by Fiori et al.³¹, 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 center^2,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.³², between total γ-ray emission measured by H.E.S.S., Milagro, Argo and HAWC and that obtained by implementing CR spectral hardening.