The contribution of Galactic TeV pulsar wind nebulae to Fermi-LAT diffuse emission

The large-scale diffuse $\gamma-$ray flux observed by Fermi-LAT in the 1-100 GeV energy range, parameterized as $\propto E^{-\Gamma}$, has a spectral index $\Gamma$ that depends on the distance from the Galactic center. This feature, if attributed to the diffuse emission produced by cosmic rays (CR) interactions with the interstellar gas, can be interpreted as the evidence of a progressive CR 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 HESS 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, viz. $\sim 30\%$ of the total diffuse $\gamma$-ray emission in the inner Galaxy. This additional component naturally accounts for a large part of the spectral index variation observed by Fermi-LAT, weakening the evidence of CR 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 (Gabici et al. 2019). The diffuse γ-ray emission, produced by 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 Pothast et al. (2018). This feature can be considered as the indirect evidence of a progressive CR spectral hardening towards the Galactic center Yang et al. (2016); Acero et al. (2016).This conclusion, however, challenges standard implementations of the CR diffusion paradigm, in which uniform diffusion throughout the Galaxy is as-sumed, and would require a more complex description of CR transport (Recchia et al. 2016;Cerri et al. 2017). It is thus extremely important to consider any possible alternative explanations of Fermi-LAT results Nava et al. (2017).
An essential step for the observational identification of CR diffuse emission, is the evaluation of the cumulative flux produced by sources which 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 to that produced by CR interactions. To investigate the role of this additional component recent works Acero et al. (2016); Pothast et al. (2018) 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 (IACT), like H.E.S.S. (Aharonian et al. 2006a), MAGIC (Aleksić et al. 2016) and VERITAS (Weekes et al. 2002), and air shower arrays, such as Argo-YBJ (Bartoli et al. 2013), Milagro (Atkins et al. 2004) and HAWC (Abeysekara et al. 2016(Abeysekara et al. , 2017(Abeysekara et al. , 2020, 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) (Abdalla et al. 2018a) or TeV halos (Linden & Buckman 2018;Sudoh et al. 2019;Giacinti et al. 2020), whose properties can be effectively constrained by observations at TeV energies (Cataldo et al. 2020;Steppa & Egberts 2020). 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 strongly 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 powerlaw with different spectral indexes β GeV and β TeV in the GeV and TeV energy domain and with a transition energy E 0 = [0.1 − 1.0] TeV located between the ranges probed by Fermi-LAT and H.E.S.S.. At high energies (E ≥ E 0 ), we allow the source spectral index to move inside the range β TeV = [1.9 − 2.5] measured by H.E.S.S. (Abdalla et al. 2018b) for identified PWNe, see Appendix A. The index β GeV is instead determined by requiring realistic values for the parameter R Φ , defined as the ratio The properties of the considered source population can be constrained by observation in the TeV energy domain. Following a previous work Cataldo et al. (2020), 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 Lorimer et al. (2006). The source density along the direction perpendicular to the Galactic plane is assumed to scale as exp (− |z| /H) 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 TeV, Min ≤ L TeV ≤ L TeV, Max Strong (2007). 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 Methods for details.
Previous analyses on the subject (Acero et al. 2016;Pothast et al. 2018) 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   . 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 (Pothast et al. 2018), the error bars represent the statistical error. The red bands represent the predicted contribution of unresolved TeV pulsar wind neubulae for α = 1.8, E0 = 0.8 TeV, βTeV = 2.4 and RΦ ranging between 250 − 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 − 1500. The gray lines show a power-law with an index of 2.7 for comparison. 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 spindown power as suggested by H.E.S.S. data Abdalla et al. (2018a). 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 corecollapse SN explosions (this analysis is only sensitive to the product Rτ . A smaller PWN formation rate can be balanced by a higher value of τ (and viceversa) with no consequences for the present discussion), i.e. R = 0.019 yr −1 ), one obtains L TeV, Max = 6.8 × 10 35 erg s −1 (L TeV, Max = 4.9 × 10 35 erg s −1 ) and τ = 0.5 × 10 3 y (τ = 1.8 × 10 3 y) for α = 1.8 (α = 1.5) (Cataldo et al. 2020). The fit to HGPS sources permits us to constrain the cumulative flux Φ tot TeV produced by PWNe population in the TeV domain with ∼ 30% statistical accuracy, being it equal to Φ tot TeV = 5.9 +1.8 −1.5 × 10 −10 cm −2 s −1 for α = 1.8. 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. Φ tot TeV = 3.8 +1.2 −1.1 × 10 −10 cm −2 s −1 for α = 1.5 that corresponds to ∼ 35% reduction with respect to the reference case α = 1.8 Cataldo et al. (2020).
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 Φ tot TeV , 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.
2.2. 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 Φ th GeV is the Fermi-LAT detection threshold. For objects contained in the Galactic plane, this is estimated as Φ th GeV = 10 −9 cm −2 s −1 Acero et al. (2015) by looking at the turnover of the observed source number as a function of the photon flux above 1 GeV (see their Fig. 24, panel (a)). By considering that the flux distribution scales as dN/dΦ We remark that the total and unresolved fluxes, Φ tot GeV and Φ NR GeV , only depend on R Φ and are independent on the assumed values for the spectral parameters E 0 and β TeV .
In the last line of Tab. 1, we give the flux Φ NR GeV 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 Φ diff GeV detected by Fermi-LAT (see second column in Tab.1) in the 1 − 100 GeV energy range and determined in Pothast et al. Pothast et al. (2018) 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 (Acero et al. 2016;Pothast et al. 2018) that have been performed under this hypothesis. Results for α = 1.5 are smaller and are reported in the last two columns of Tab. 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. Pothast et al. (2018). 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 Tab. 1, the unresolved contribution becomes more relevant in the central rings, due 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 Pothast et al. (2018). 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 parameters 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 = 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 Tab. 2 and displayed in Fig. 3, also take into account effects of possible variations of E 0 and β TeV .
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. Pothast et al. (2018) with a single power-law ∝ E −Γ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 Tab.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 (Acero et al. 2016;Yang et al. 2016) performed by using a similar approach. One can get a visual perception of the situation by comparing the green dashed lines with the grey 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 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 0 and β TeV on the unresolved PWNe emission are shown in the Supplementary Material, see Fig.6. 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 Sect. Method (see Eq. 8). As a result of this, PWNe unresolved emission in the inner Galaxy is either relatively large or has an 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 ∝ E −Γ 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 are larger and closer to the value measured at the Sun position with respect to those obtained in previous analyses (Pothast et al. 2018;Acero et al. 2016;Yang et al. 2016) that do not take unresolved sources into account, see Tab. 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 0 and β TeV and provide a conservative estimate of the total systematical error for CR diffuse emission, as it discussed in the next paragraph.

Spectral index
In order to estimate the uncertainties in our approach, we repeat our analysis for different combinations of the spectral parameters (R Φ , E 0 , β TeV ). The final results of our analysis are given in Tab. 2. Here, the first errors describe systematic uncertainties, evaluated as the maximal variations of Γ BF that are obtained when (R Φ , E 0 , β TeV ) are simultaneously varied in the 3-dim parameter space defined by the ranges R Φ = [250 − 1500], E 0 = [0.1 − 1.0] TeV and β TeV = [1.9 − 2.5]. 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 0 = 1.0 TeV (E 0 = 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 Tab. 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 (Abdalla et al. 2018b).
In addition to the reference case α = 1.8 (third column) that is obtained by using the luminosity distribution index considered by previous analyses (Acero et al. 2016;Pothast et al. 2018), 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 Tab. 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 (Kelner et al. 2006). Black points show the results of our work that are compared to those given by Pothast et al. Pothast et al. (2018)

(red points) and Acero et al. Acero et al. (2016) (orange points).
We also show with grey points the results obtained by Peron et al. Peron et al. (2021) 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. Pothast et al. (2018) and Peron et al. Peron et al. (2021). 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 (Strong et al. 2004) (dashed blue line. The theoretical CR distribution is shown e.g. in Fig. 8 of Acero et al. (2016) where the specific GALPROP configuration is also given. It basically coincides with the solution of 3D isotropic diffusion equation with uniform diffusion coefficient, stationary CR injection and infinite smearing radius, as it is e.g. shown in Fig. 1 of Cataldo et al. (2019).). The inclusion of unresolved PWNe strongly 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 account 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 PWNe. This phenomenological description could be improved in the future adopting a more refined PWN spectrum as the one in Fiori et al. (Fiori et al. 2022) 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 behaviour 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 modulate 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 behaviour of the total diffuse emission as an indirect evidence for CR spectral hardening toward the Galactic center (Acero et al. 2016;Yang et al. 2016;Pothast et al. 2018). 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 strongly 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. Cataldo et al. (2019), between total γ-ray emission measured by H.E.S.S., Milagro, Argo and HAWC and that obtained by implementing CR spectral hardening.

METHOD
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: 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 behaviour at GeV and TeV energies. We take this into account by parameterizing it with a broken powerlaw with different spectral indexes β GeV and β TeV in the GeV and TeV energy domain and with a transition energy E 0 located between the ranges probed by Fermi-LAT and H.E.S.S.. Even if our approach is completely phenomenological, the postulated spectral behaviour 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 Sudoh et al. (2019). In this scenario, the observed gamma-ray emission is produced by IC scattering of HE electron and positrons on background photons (CMB, 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 (Bucciantini et al. 2011;Torres et al. 2014;Sudoh et al. 2021). In the assumption of a broken powerlaw for the gamma ray spectrum, the flux ratio R Φ , the energy break E 0 and the two spectral indexes β GeV and β TeV are not independent and are related by the following expression: where ϵ inf GeV ≡ (1.0 GeV/E 0 ) and ϵ sup GeV ≡ (100 GeV/E 0 ) (ϵ inf TeV ≡ (1.0 TeV/E 0 ) and ϵ sup TeV ≡ (100 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 behaviour at high energies is known (i.e. β TeV and E 0 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 0 , β TeV ) in the 3-dim parameter space defined by the ranges R Φ = [250 − 1500], E 0 = [0.1 − 1.0] TeV and β TeV = [1.9 − 2.5]. The spectral index β GeV of GeV emission is determined as a function of (R Φ , E 0 , β 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 0 , β 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 Tab. 1 of the Supplementary Mate-rial). 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 0 , 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) 2 ). 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 ≃ 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 (Cataldo et al. 2020). 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: 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 (Gaensler & Slane 2006) or TeV Halos (Linden & Buckman 2018), 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 ≃ R SN = 0.019 yr −1 Diehl et al. (2006). If gammaray emission is powered by pulsar activity, the TeVluminosity can be connected to the pulsar spin-down power, i.e.: L TeV = λĖ (10) where λ ≤ 1 and:Ė =Ė 0 1 + t τ sd −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 (λ ∼ const), the exponent in Eq. 9 is γ = 2, that corresponds to a source luminosity function Y TeV (L TeV ) ∝ L −1.5 TeV . The possibility of λ being correlated to the spin-down power, i.e. λ = λ 0 (Ė/Ė 0 ) δ , was suggested by Abdalla et al. Abdalla et al. (2018a) that found L TeV = λĖ ∝Ė 1+δ 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 TeV (L TeV ) ∝ L −1.8 TeV . We consider this last scenario (α = 1.8) as our reference case, conforming to previous analyses on the subject (Acero et al. 2016;Pothast et al. 2018). 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 Y TeV (F TeV ) that describes the source emissivity distribution. This is related to the luminosity function by the expression Y TeV (F TeV ) = E TeV Y TeV (F TeV E TeV ). By using Eq. 3, we see that the emissivity distribution is not modified, if the ratio F TeV, Max ≡ L TeV, Max /E 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 (see Fig. 4). It should be remarked that, while PWNe provide 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 (Abdollahi et al. 2020). 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 (labelled 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.
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 TeV, Max = 4.9 · 10 35 erg s −1 (L TeV, Max = 6.8 · 10 35 erg s −1 ) and τ = 1.8·10 3 y (τ = 0.5·10 3 y) for α = 1.5 (α = 1.8) given in Cataldo et al. Cataldo et al. (2020) 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 Φ = 10 3 (βTeV−1) ∼ 10 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 TeV, Max and τ by: where N = R τ (α − 1) and ∆ ≡ L TeV, Max /L 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: where the parameter ξ, which is defined as represents the fraction of sources of the considered population which are included in the OW, while the quantity 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 TeV, Max = L TeV, Max /E TeV . The total flux produced in the GeV domain can be calculated as a function of the parameter R Φ (for fixed values of L 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 TeV, Max and τ by using: where ρ(r) ≡ OW dΩ ρ(r, n). The above expression can be recasted in terms of the source emissivity distribution as: By considering that the function Y TeV (F TeV ) only depends on the parameter F 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 TeV, Max in such a way that the ratio F TeV, Max = L TeV, Max /E TeV is kept constant, with no effects on the predicted source distribution dN/dΦ TeV , on the cumulative flux Φ tot TeV 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 0 and β TeV . Faint sources that produce a flux at Earth below the Fermi-LAT observation threshold Φ th GeV are not resolved and contribute to the large scale diffuse emission from the Galaxy. This contribution is evaluated by using Eq. 6.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

ACKNOWLEDGEMENT
We thanks Daniele Gaggero and Mart Pothast for providing us the Fermi-LAT data points of the total gamma ray diffuse emission. The work of GP and FLV is partially supported by the research grant number 2017W4HA7S "NAT-NET: Neutrino and Astroparticle Theory Network" under the program PRIN 2017 funded by the Italian Ministero dell'Istruzione, dell'Universita' e della Ricerca (MIUR).

AUTHOR CONTRIBUTIONS
All the Authors equally contributed to this work.

COMPETING INTERESTS
The authors declare no competing interests. Table 1. The cumulative gamma fluxes due to unresolved PWNe. The cumulative flux (Φ NR GeV ) of unresolved TeV PWNe in the GeV domain for α = 1.8 and α = 1.5 and for the two extreme values of RΦ allowed in our analysis. In brackets, we give the percentage of unresolved sources emission with respect to the total diffuse γ-ray flux (Φ diff GeV ) measured by Fermi-LAT in each galactocentric ring and in the latitude window |b| < 20.25 • .
Ring ( Figure 5. Comparison between observational and theoretical cumulative number of sources in the TeV energy range: The cumulative distribution of source in HGPS is shown by a black thick line. The blue (red) dashed line represent theoretical predictions from our population model for the best fit values of the maximal luminosity LTeV, Max = 6.8 · 10 35 erg s −1 (LTeV, Max = 5.0 · 10 35 erg s −1 ) and spin-down timescale τ = 0.5 · 10 3 y (τ = 1.7 · 10 3 y) with luminosity index α = 1.8 (α = 1.5).
note that the unresolved PWNe emission at 50 GeV can be increased (decreased) by a factor ∼ 3 (∼ 2) with respect to the reference case when (R Φ , E 0 , β TeV ) are simultaneously varied in the 3-dim parameter space defined by the ranges R Φ = [250 − 1500], E 0 = [0.1 − 1.0] TeV and β TeV = [1.9 − 2.5]. For completeness, we also show with a dashed line the predicted emission that was used to obtain the central values for Γ BF quoted in Tab. 2. This is 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 (Abdalla et al. 2018b).