The contribution of Galactic TeV pulsar wind nebulae to Fermi large area telescope diffuse emission

Abstract

The large-scale diffuse γ − ray flux observed by Fermi Large Area Telescope (Fermi-LAT) in the 1–100 GeV energy range, parameterized as ∝ E^−Γ, has a spectral index Γ that depends on the distance from the Galactic center. This feature, if attributed to the diffuse emission produced by cosmic rays interactions with the interstellar gas, can be interpreted as the evidence of a progressive cosmic ray spectral hardening towards the Galactic center. This interpretation challenges the paradigm of uniform cosmic rays diffusion throughout the Galaxy. We report on the implications of TeV Pulsar Wind Nebulae observed by the High Energy Stereoscopic System (H.E.S.S.) Galactic Plane Survey in the 1–100 TeV energy range for the interpretation of Fermi-LAT data. We argue that a relevant fraction of this population cannot be resolved by Fermi-LAT in the GeV domain providing a relevant contribution to the large-scale diffuse emission, ranging within ~4%–40% of the total diffuse γ-ray emission in the inner Galaxy. This additional component may account for a large part of the spectral index variation observed by Fermi-LAT, weakening the evidence of cosmic ray spectral hardening in the inner Galaxy.

Introduction

Cosmic rays (CRs) with energy below ~1 PeV are believed to originate in the Milky Way and to spread in the entire Galaxy due to diffusion in local magnetic fields¹. The diffuse γ-ray emission, produced by the interaction of CRs with the gas contained in the galactic disk, carries information on the energy distribution of CRs in different regions of the Galaxy.

Recent observations at GeV energies performed by Fermi-LAT suggest that the hadronic diffuse gamma-ray emission, parameterized as ∝E^−Γ, has a spectral index Γ in the inner Galaxy which is smaller by an amount ~−0.2 than the value observed at the Sun position². This feature can be considered as indirect evidence of a progressive CR spectral hardening towards the Galactic center^3,4. This conclusion, however, challenges standard implementations of the CR diffusion paradigm, in which uniform diffusion throughout the Galaxy is assumed, and would require a more complex description of CR transport^5,6. It is thus extremely important to consider any possible alternative explanations of Fermi-LAT results⁷.

An essential step for the observational identification of CR diffuse emission is the evaluation of the cumulative flux produced by sources that are too faint to be resolved by Fermi-LAT. These sources are not individually detected but give rise to a large-scale diffuse flux superimposed on that produced by CR interactions.

To investigate the role of this additional component recent works^2,4 performed a source population study concluding that the diffuse flux associated to unresolved sources is not large enough to explain the spectral anomaly being below 3% at 1 GeV (20% at ≃100 GeV) of the total observed diffuse emission. Both studies are tuned on the 3FGL catalog. As a consequence, they reproduce the population of Galactic sources observed in the GeV energy domain which is largely dominated by Pulsars. These objects have γ-ray spectra with exponential cutoff at few GeV and are expected to provide a negligible contribution to observed emission at E ≥ 10 GeV.

In the last decade, Imaging Atmospheric Cherenkov Telescopes, like H.E.S.S.⁸, MAGIC⁹, and VERITAS¹⁰, and air shower arrays, such as Argo-YBJ¹¹, Milagro¹², and HAWC^13,14,15, provided a detailed description of Galactic γ-ray emission in the energy range 0.1–100 TeV. The emerging picture is that TeV Galactic sky is dominated by a population of bright sources powered by pulsar activity, such as pulsar wind nebulae (PWNe)¹⁶ or TeV halos^17,18,19, whose properties can be effectively constrained by observations at TeV energies^20,21. These objects are clearly expected to emit also in the GeV energy domain where, however, population studies are more difficult because different kinds of sources dominate the observed emission.

In this paper, we took advantage of the constraints provided by H.E.S.S. Galactic Plane Survey (HGPS) to discuss the implications of TeV PWNe for the interpretation of Fermi-LAT data in the GeV domain. We quantify the contribution of unresolved TeV PWNe to large-scale diffuse emission observed by Fermi-LAT at different distances from the Galactic center. We show that the inclusion of this additional component can affect the reconstructed CR energy distribution from Fermi-LAT data, weakening the evidence of a progressive hardening of the cosmic-ray spectrum toward the Galactic center.

Results and discussion

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

$${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 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:

$$\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 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:

$${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 }$²⁶. 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^−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 × 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:

$${{{\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 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.

Table 1 The cumulative gamma fluxes due to unresolved Pulsar Wind Nebulae.

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.

Fig. 1: The diffuse gamma-ray emission as a function of energy in different galactocentric rings.

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.

Fig. 2: Pulsar wind nebulae cumulative spectrum.

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.

Table 2 Gamma ray spectral indexes.

Fig. 3: Gamma-rays emissivity and cosmic ray proton spectral index in different galactocentric rings.

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.

Methods

Flux and luminosity ratios

The sources considered in this work are expected to contribute to observations both in the GeV and TeV energy domains. We indicate with Φ_GeV (Φ_TeV) and L_GeV (L_TeV) the integrated source flux and luminosity in the energy range 1–100 GeV (1–100 TeV) probed by Fermi-LAT (H.E.S.S.). We assume for simplicity that all the sources in the considered population have approximately the same emission spectrum. This automatically implies that the ratio R_Φ ≡ Φ_GeV/Φ_TeV between fluxes emitted in different energy domains by a given source is fixed. The relationship between intrinsic luminosity and flux produced at Earth is generically written as:

$${{{\Phi }}}_{{{{{{{{\rm{X}}}}}}}}}=\frac{{L}_{{{{{{{{\rm{X}}}}}}}}}}{4\pi {r}^{2}{E}_{{{{{{{{\rm{X}}}}}}}}}}$$

(7)

where r is the source distance, E_X is the average energy of emitted photons and X = GeV, TeV indicates the considered energy range. It is also useful to define the integrated emissivity F_X ≡ L_X/E_X that corresponds to the total amount of photons emitted per unit time by a given source in the X = GeV, TeV energy domain.

Source spectrum

The source emission spectrum φ(E) can have a different behavior at GeV and TeV energies. We take this into account by parameterizing it with a broken power-law with different spectral indexes β_GeV and β_TeV in the GeV and TeV energy domain and with a transition energy E₀ located between the ranges probed by Fermi-LAT and H.E.S.S.. Even if our approach is completely phenomenological, the postulated spectral behavior is expected from a theoretical point of view. We are indeed considering the hypothesis that most of the bright TeV sources are young PWNe and/or TeV halos¹⁸. In this scenario, the observed gamma-ray emission is produced by inverse Compton scattering of high-energy electrons and positrons on background photons (cosmic microwave background, starlight, infrared). In the Thompson regime, this naturally produces hard gamma-ray emission with spectral index β ~ (p + 1)/2 where p is the electron/positron spectral index. At TeV energy, it produces instead a softer gamma-ray spectrum either due to the Klein-Nishina regime β ~ (p + 1) or to electron/positron energy losses^23,24,33. In the assumption of a broken power-law for the gamma ray spectrum, the flux ratio R_Φ, the energy break E₀ and the two spectral indexes β_GeV and β_TeV are not independent and are related by the following expression:

$${R}_{{{\Phi }}}=\frac{1-{\beta }_{{{{{{{{\rm{TeV}}}}}}}}}}{1-{\beta }_{{{{{{{{\rm{GeV}}}}}}}}}}\,\frac{\left[{({\epsilon }_{{{{{{{{\rm{GeV}}}}}}}}}^{\sup })}^{1-{\beta }_{{{{{{{{\rm{GeV}}}}}}}}}}-{({\epsilon }_{{{{{{{{\rm{GeV}}}}}}}}}^{\inf })}^{1-{\beta }_{{{{{{{{\rm{GeV}}}}}}}}}}\right]}{\left[{({\epsilon }_{{{{{{{{\rm{TeV}}}}}}}}}^{\sup })}^{1-{\beta }_{{{{{{{{\rm{TeV}}}}}}}}}}-{({\epsilon }_{{{{{{{{\rm{TeV}}}}}}}}}^{\inf })}^{1-{\beta }_{{{{{{{{\rm{TeV}}}}}}}}}}\right]}$$

(8)

where ${\epsilon }_{{{{{{{{\rm{GeV}}}}}}}}}^{\inf }\equiv (1.0\,{{{{{{{\rm{GeV}}}}}}}}/{E}_{0})$ and ${\epsilon }_{{{{{{{{\rm{GeV}}}}}}}}}^{\sup }\equiv (100\,{{{{{{{\rm{GeV}}}}}}}}/{E}_{0})$ (${\epsilon }_{{{{{{{{\rm{TeV}}}}}}}}}^{\inf }\equiv (1.0\,{{{{{{{\rm{TeV}}}}}}}}/{E}_{0})$ and ${\epsilon }_{{{{{{{{\rm{TeV}}}}}}}}}^{\sup }\equiv (100\,{{{{{{{\rm{TeV}}}}}}}}/{E}_{0})$) are the lower and upper bounds of the GeV (TeV) energy domains. The above relationship implements mathematically the fact that, if the source spectral behavior at high energies is known (i.e., β_TeV and E₀ are fixed), then the flux ratio R_Φ is an increasing function of β_GeV. In other words, the harder is the spectrum at GeV energies, the smaller is the integrated flux in the GeV domain. In our analysis, we vary the parameters (R_Φ, E₀, β_TeV) 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]$. The spectral index β_GeV of GeV emission is determined as a function of (R_Φ, E₀, β_TeV) by inverting Eq. (8) and it globally spans the range β_GeV = 1.06 –2.19. By repeating our analysis for different combinations (R_Φ, E₀, β_TeV), we determine the stability of our results against the assumed source spectrum and we estimate the systematic uncertainties produced by the scatter of spectral properties in the PWNe population.

The assumed source spectrum can be validated by considering the ensemble of PWNe firmly identified both in the 4FGL-DR2 and HGPS catalogs (12 objects, reported in Supplementary Table 1 of the Supplementary Note 1). Within this sample, the average values for R_Φ and β_GeV are 1122 and 1.89, respectively. These values fall inside the ranges of variation for these parameters considered in our analysis. The spectral break position E₀, estimated as the crossing point of the spectral fits given by Fermi-LAT and H.E.S.S. in the GeV and TeV domain respectively, falls inside the range 0.1 TeV–few TeV for all the sources. Finally, the characteristic spectral energy distribution (SED) of the PWNe population is estimated by calculating the cumulative spectrum of sources included both in 4FGL-DR2 and HGPS catalogs (black line in Fig. 2). This is obtained by considering the best-fit spectra of each source, as given by Fermi-LAT and H.E.S.S. for the respective energy ranges. We have included sources with a known distance D whose spectra have been weighted proportionally to their intrinsic luminosity (namely, they have been scaled by a factor (D/1kpc)²). The shaded bands are obtained by propagating the statistical errors on the spectral parameters of each source and by assuming a 20% (30%) systematic uncertainty for the flux normalization and 0.1 (0.2) systematic uncertainty for the spectral indexes in 4FGL-DR2 (HGPS) catalog, respectively. The cumulative spectrum is well reproduced by a broken power-law with an energy break at E₀ ≃ 0.8 TeV, spectral indexes β_TeV = 2.4 and β_GeV = 1.89, and R_Φ ≃ 770 (red line in Fig. 2). This functional form is allowed within the parameter space considered in this paper.

Luminosity distribution

In the following, we focus on the TeV-luminosity function since this can be effectively constrained by HGPS observational results²⁰. The function Y_TeV(L_TeV) is parameterized as described in Eq. (3). This distribution is naturally obtained for a population of fading sources with intrinsic luminosity that decreases over a time scale τ according to:

$${L}_{{{{{{{{\rm{TeV}}}}}}}}}(t)={L}_{{{{{{{{\rm{TeV}}}}}}}},\max }{\left(1+\frac{t}{\tau }\right)}^{-\gamma }$$

(9)

where t indicates the time passed since source formation. In this assumption, the exponent of the luminosity distribution is given by α = 1/γ + 1.

The above description can be applied to potential TeV sources in the Galaxy, such as PWNe³⁴ or TeV Halos¹⁷, which are connected with the explosion of core-collapse SN and the formation of a pulsar. The birth rate of these objects is similar to that of SN explosions in our Galaxy, i.e., $R\simeq {R}_{{{{{{{{\rm{SN}}}}}}}}}=0.019\,{{{{{{{{\rm{yr}}}}}}}}}^{-1}$³⁵. If gamma-ray emission is powered by pulsar activity, the TeV-luminosity can be connected to the pulsar spin-down power, i.e.,:

$${L}_{{{{{{{{\rm{TeV}}}}}}}}}=\lambda \,\dot{E}$$

(10)

where λ ≤ 1 and:

$$\dot{E}={\dot{E}}_{0}{\left(1+\frac{t}{{\tau }_{{{{{{{{\rm{sd}}}}}}}}}}\right)}^{-2}$$

(11)

for energy loss dominated by magnetic dipole radiation (braking index n = 3). This implies that the fading timescale is determined by the pulsar spin-down time scale, i.e., τ = τ_sd. Moreover, if the efficiency of TeV emission does not depend on time ($\lambda \sim {{{{{{{\rm{const}}}}}}}}$), the exponent in Eq. (9) is γ = 2, that corresponds to a source luminosity function ${Y}_{{{{{{{{\rm{TeV}}}}}}}}}({L}_{{{{{{{{\rm{TeV}}}}}}}}})\propto {L}_{{{{{{{{\rm{TeV}}}}}}}}}^{-1.5}$. The possibility of λ being correlated to the spin-down power, i.e., $\lambda ={\lambda }_{0}{(\dot{E}/{\dot{E}}_{0})}^{\delta }$, was suggested by Abdalla et al.¹⁶ that found ${L}_{{{{{{{{\rm{TeV}}}}}}}}}=\lambda \,\dot{E}\propto {\dot{E}}^{1+\delta }$ with 1 + δ = 0.59 ± 0.21 by studying a sample of PWNe in the HPGS catalog. In this case, one obtains γ ≃ 1.2 in Eq. (9) that corresponds to a source luminosity function ${Y}_{{{{{{{{\rm{TeV}}}}}}}}}({L}_{{{{{{{{\rm{TeV}}}}}}}}})\propto {L}_{{{{{{{{\rm{TeV}}}}}}}}}^{-1.8}$. We consider this last scenario (α = 1.8) as our reference case, conforming to previous analyses on the subject^2,4. In order to discuss thoroughly the dependence of our results on the performed assumptions, we also consider, however, the alternative hypothesis α = 1.5.

It is finally useful to introduce the function ${{{{{{{{\mathcal{Y}}}}}}}}}_{{{{{{{{\rm{TeV}}}}}}}}}\left({F}_{{{{{{{{\rm{TeV}}}}}}}}}\right)$ that describes the source emissivity distribution. This is related to the luminosity function by the expression ${{{{{{{{\mathcal{Y}}}}}}}}}_{{{{{{{{\rm{TeV}}}}}}}}}\left({F}_{{{{{{{{\rm{TeV}}}}}}}}}\right)={E}_{{{{{{{{\rm{TeV}}}}}}}}}\,{Y}_{{{{{{{{\rm{TeV}}}}}}}}}({F}_{{{{{{{{\rm{TeV}}}}}}}}}{E}_{{{{{{{{\rm{TeV}}}}}}}}})$. By using Eq. (3), we see that the emissivity distribution is not modified, if the ratio ${F}_{{{{{{{{\rm{TeV}}}}}}}},\max }\equiv {L}_{{{{{{{{\rm{TeV}}}}}}}},\max }/{E}_{{{{{{{{\rm{TeV}}}}}}}}}$ is kept constant.

Consistency among HGPS and Fermi-LAT catalogs

The adopted range for the flux ratio parameter R_Φ can be further validated by comparing the predicted source flux distribution in the GeV domain with the 4FGL-DR2 catalog (Fig. 4). It should be remarked that, while PWNe provides the prominent contribution of the observed emission at TeV energies, they are instead a subdominant component in the GeV domain. The 4FGL-DR2 catalog includes 5788 sources which are mostly extragalactic objects³⁶. The total number of identified and/or associated Galactic sources is 486. The largest source class, including 271 objects, is given by pulsars that typically have soft emission spectra with cut-off at few GeV and are not expected to contribute to the population of TeV emitting sources potentially detectable by HGPS. In addition to pulsars, the 4FGL-DR2 catalog encompasses 18 PWNe, 43 SNRs, and 96 objects (labeled as SPP) of unknown nature but overlapping with known SNRs or PWNe. The magenta line in Fig. 4 corresponds to the cumulative number N(Φ_GeV) of PWNe with flux larger than Φ_GeV in the latitude range ∣b∣ ≤ 20.25^∘ while the black line also includes SPP sources. The SPP source class is not expected to fully correspond to the population considered in this work; it can be however regarded as an upper limit for theoretical calculations.

Fig. 4: Pulsar wind nebulae flux distribution in the GeV domain.

We report with shaded bands the cumulative number N(Φ_GeV) of sources with fluxes larger than Φ_GeV predicted in our model and in the latitude range ∣b∣ ≤ 20.25^∘. In particular, the red (blue) band is obtained by assuming α = 1.5 (α = 1.8) and by considering 250 ≤ R_Φ ≤ 1500. The magenta line represents the observed cumulative number of pulsar wind nebulae (PWNe) with fluxes larger than Φ_GeV in the 4FGL-DR2 catalog. The black line also includes SPP sources with potential association with supernova remnant or pulsar wind nebulae. The gray vertical line represents the Fermi-LAT sensitivity threshold for objects in the Galactic plane.

The two shaded bands in Fig. 4 show theoretical predictions of our population model for two different values of the power-law index of the luminosity function (α = 1.5 and 1.8). Namely, the red (blue) shaded band is obtained by assuming the best-fit values ${L}_{{{{{{{{\rm{TeV}}}}}}}},\max }=4.9\times 1{0}^{35}\,{{{{{{{\rm{erg}}}}}}}}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}$ (${L}_{{{{{{{{\rm{TeV}}}}}}}},\max }=6.8\times 1{0}^{35}\,{{{{{{{\rm{erg}}}}}}}}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}$) and τ = 1.8 × 10³â€‰y (τ = 0.5 × 10³â€‰y) for α = 1.5 (α = 1.8) given in Cataldo et al.²⁰ and by varying the flux ratio in the range 250 ≤ R_Φ ≤ 1500. The lower bound for the flux ratio (i.e., R_Φ = 250) is obtained by requiring that sources are not underpredicted (within statistical fluctuations) in the flux region where the catalog can be considered complete. More precisely, it corresponds to assuming 6 sources with flux larger than 5 × 10^−9 cm^−2 s^−1 to be compared with 9 PWNe in the 4FGL-DR2 catalog. The upper bound (i.e., R_Φ = 1500) is instead obtained by requiring that very bright sources are not overpredicted and corresponds to assuming 3 sources with flux larger than 5 × 10^−8 cm^−2 s^−1 to be compared with no observed PWNe + SPP sources in the 4FGL-DR2 catalog.

In general, we see that a reasonable agreement exists with theoretical expectations, supporting the phenomenological description adopted in this paper. We also note that the performed comparison provides by itself a proof that the average spectral index β_GeV of PWNe at GeV energies should be smaller than the value β_TeV = [2.3–2.4] observed in the TeV domain. Indeed, if we assume that source spectrum is described by undistorted power-law with spectral index β_TeV ~ [2.3–2.4], we obtain ${R}_{{{\Phi }}}=1{0}^{3({\beta }_{{{{{{{{\rm{TeV}}}}}}}}}-1)} \sim 1{0}^{4}$. Considering that bright sources in the HGPS catalog have fluxes Φ_TeV ~ 10^−11 cm^−2 s^−1, we should expect an ensemble of sources with fluxes Φ_GeV ~ 10^−7 cm^−2 s^−1 in the GeV domain. This is not observed by Fermi-LAT, indicating that TeV galactic sources typically have a spectral break and harder emission spectrum below ~1 TeV. Coherently with this conclusion, most of the PWNe in the 4FGL-DR2 catalog have a spectral indexes ≤ 2 at GeV energies.

Total luminosity and flux

The total luminosity produced by the considered population in the TeV domain is given as a function of ${L}_{{{{{{{{\rm{TeV}}}}}}}},\max }$ and τ by:

$${L}_{{{{{{{{\rm{TeV}}}}}}}}}^{{{{{{{{\rm{MW}}}}}}}}}=\frac{{{{{{{{\mathcal{N}}}}}}}}\,{L}_{{{{{{{{\rm{TeV}}}}}}}},\max }}{\left(2-\alpha \right)}\left[1-{{{\Delta }}}^{\alpha -2}\right]$$

(12)

where ${{{{{{{\mathcal{N}}}}}}}}=R\,\tau \,(\alpha -1)$ and ${{\Delta }}\equiv {L}_{{{{{{{{\rm{TeV}}}}}}}},\max }/{L}_{{{{{{{{\rm{TeV}}}}}}}},\min }$. Unless otherwise specified, we quote results obtained for Δ → ∞ that can be easily recalculated by using the above equation, if other values are considered. The flux in the TeV domain produced at Earth by all sources included in a fixed observational window (OW) can be expressed as:

$${{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}^{{{{{{{{\rm{tot}}}}}}}}}=\xi \,\frac{{L}_{{{{{{{{\rm{TeV}}}}}}}}}^{{{{{{{{\rm{MW}}}}}}}}}}{4\pi {E}_{{{{{{{{\rm{TeV}}}}}}}}}}\,\langle {r}^{-2}\rangle$$

(13)

where the parameter ξ, which is defined as

$$\xi \equiv {\int}_{{{{{{{{\rm{OW}}}}}}}}}{d}^{3}r\,\rho ({{{{{{{\bf{r}}}}}}}}),$$

(14)

represents the fraction of sources of the considered population which are included in the OW, while the quantity

$$\langle {r}^{-2}\rangle \equiv \frac{1}{\xi }{\int}_{{{{{{{{\rm{OW}}}}}}}}}{d}^{3}r\,\rho ({{{{{{{\bf{r}}}}}}}})\,{r}^{-2}$$

(15)

is the average value of their inverse square distance. By considering Eqs. (12) and (13), we see that the total flux produced by the considered population in the TeV domain is directly determined by the maximal emissivity ${F}_{{{{{{{{\rm{TeV}}}}}}}},\max }={L}_{{{{{{{{\rm{TeV}}}}}}}},\max }/{E}_{{{{{{{{\rm{TeV}}}}}}}}}$.

The total flux produced in the GeV domain can be calculated as a function of the parameter R_Φ (for fixed values of ${L}_{{{{{{{{\rm{TeV}}}}}}}},\max }$ and τ) and it is given by Eq. (4).

Source flux distributions

The source flux distribution in the TeV domain dN/dΦ_TeV can be calculated as a function of ${L}_{{{{{{{{\rm{TeV}}}}}}}},\max }$ and τ by using:

$$\frac{dN}{d{{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}}=\int dr\,4\pi {r}^{4}{E}_{{{{{{{{\rm{TeV}}}}}}}}}\,{Y}_{{{{{{{{\rm{TeV}}}}}}}}}(4\pi {r}^{2}{E}_{{{{{{{{\rm{TeV}}}}}}}}}{{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}})\,\overline{\rho }(r),$$

(16)

where $\overline{\rho }(r)\equiv {\int}_{{{{{{{{\rm{OW}}}}}}}}}d{{\Omega }}\,\rho (r,{{{{{{{\bf{n}}}}}}}})$. The above expression can be recasted in terms of the source emissivity distribution as:

$$\frac{dN}{d{{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}}=\int dr\,4\pi {r}^{4}\,{{{{{{{{\mathcal{Y}}}}}}}}}_{{{{{{{{\rm{TeV}}}}}}}}}(4\pi {r}^{2}{{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}})\,\overline{\rho }(r),$$

(17)

By considering that the function ${{{{{{{{\mathcal{Y}}}}}}}}}_{{{{{{{{\rm{TeV}}}}}}}}}({F}_{{{{{{{{\rm{TeV}}}}}}}}})$ only depends on the parameter ${F}_{{{{{{{{\rm{TeV}}}}}}}},\max }$, we can understand the effects of a variation of β_TeV in our analysis. Indeed, a modification of β_TeV reflects into a variation of the photon average energy E_TeV. This can be reabsorbed by a shift of ${L}_{{{{{{{{\rm{TeV}}}}}}}},\max }$ in such a way that the ratio ${F}_{{{{{{{{\rm{TeV}}}}}}}},\max }={L}_{{{{{{{{\rm{TeV}}}}}}}},\max }/{E}_{{{{{{{{\rm{TeV}}}}}}}}}$ is kept constant, with no effects on the predicted source distribution dN/dΦ_TeV, on the cumulative flux ${{{\Phi }}}_{{{{{{{{\rm{TeV}}}}}}}}}^{{{{{{{{\rm{tot}}}}}}}}}$ in the TeV domain and on the quality of the fit to HGPS catalog.

Finally, the source flux distribution dN/dΦ_GeV in the GeV domain is connected to that in TeV domain by the R_Φ parameter, according to Eq. (5). Note that the distribution dN/dΦ_GeV is predicted independently on the assumed values of the spectral parameters E₀ and β_TeV. Faint sources that produce a flux at Earth below the Fermi-LAT observation threshold ${{{\Phi }}}_{{{{{{{{\rm{GeV}}}}}}}}}^{{{{{{{{\rm{th}}}}}}}}}$ are not resolved and contribute to the large-scale diffuse emission from the Galaxy. This contribution is evaluated by using Eq. (6).