The contribution of Galactic TeV pulsar wind nebulae to Fermi-LAT diﬀuse emission

,


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 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 the standard CR diffusion paradigm, in which uniform diffusion throughout the Galaxy is assumed, and would require non-standard or anomalous CR transport mech-anisms, see e.g. (Recchia et al. 2016;Cerri et al. 2017). It is thus extremely important to consider any possible alternative explanations of Fermi-LAT results (e.g. 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 Acero et al. (2016) and 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 HESS (Aharonian et al. 2006a), ... (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, see e.g. (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.

MAGIC
In this paper, we took advantage of the constraints provided by HESS 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
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 (HESS). 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 HESS. At high energies (E ≥ E 0 ), we take the average spectrum observed by HESS (Abdalla et al. 2018b) as a reference, i.e. we assume that all sources have β TeV = [2.3 − 2.4]. 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 Sect. 4, we ob-tain a consistent description of HGPS and 4FGL-DR2 Fermi-LAT catalogs for R Φ = [250 − 1500] that corresponds to β GeV = [1.18 − 2.0] (see Eq. 7). The assumed source spectrum can be further validated by considering the average observational properties of PWNe observed by Fermi-LAT and HESS in the GeV/TeV domain, see Section 4 for details. Moreover, the corresponding spectral shapes are consistent with theoretical predictions of γ-ray emission from PWNe, as e.g. discussed by Bucciantini et al. (2011);Torres et al. (2014). The properties of the considered source population can be constrained by observation in the TeV energy domain. Following 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 parameterized by 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 , see e.g. Strong (2007). This distribution 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. As working hypotheses, two different values for the luminosity index, α = 1.5 and α = 1.8, are considered. 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 1 , i.e. R = 0.019 yr −1 ), one obtains 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) (Cataldo et al. 2020). These values have been obtained Table 1. We provide the cumulative flux of resolved (Φ R GeV ) and unresolved (Φ NR GeV ) TeV PWNe in the GeV domain for α = 1.8 and for the two extreme values of RΦ allowed in our analysis. It is shown in bracket the percentage of unresolved sources with respect to the total diffuse γ-ray emission measured by Fermi-LAT in each galactocentric ring and in the latitude window |b| < 20.25 • .
The determination of L TeV, Max and τ by HGPS data allows us to constrain the properties of PWNe population in the TeV domain. As an example, the total flux Φ tot TeV produced at Earth by the considered sources can be calculated as a function of L TeV, Max and τ by using Eq. (12), as discussed in Section 4. The total flux produced at Earth by TeV PWNe in the GeV domain depends also on the parameter R Φ and it is given by: A fraction of this flux is emitted by sources which are too faint to 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 dN/dΦ GeV is the source flux distribution in the GeV domain while Φ 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 in 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)).
In the last line of Tab.1, we give the cumulative flux produced by TeV PWNe that are not resolved (resolved) by Fermi-LAT for the two extreme values R Φ = 250 and 1500. These fluxes are compared with the total large scale diffuse emission Φ diff GeV detected by Fermi-LAT (see second column in Tab.1) in the 1 − 100 GeV energy range and determined by Pothast et al. (2018) by using 9.3 years of Fermi-LAT Pass 8 data 2 . We see that unresolved emission by PWNe corresponds to a fraction ∼ 3% (for R Φ = 250) and ∼ 11% (for R Φ = 1500) of the total large scale diffuse emission. The above predictions 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 reported as additional material in the last section.
In order to probe the radial dependence of the PWNe contribution, we repeat our calculations by considering the Galactocentric rings adopted by 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 10 50 100 5. × 10 -8 10 50 100 1. × 10 -9 5. × 10 -9 1. × 10 -8 5. × 10 -8 16.5-50 kpc Γ BF =2.78 Figure 1. Black data points show the total diffuse γ-ray emission measured by Fermi-LAT in each galactocentric ring Pothast et al. (2018). The red lines represent the predicted contribution of unresolved TeV PWNe for α = 1.8, E0 = 0.3 TeV, βTeV = 2.3 and RΦ in the range 250 − 1500. Green lines show the diffuse CR emission inferred by fitting the data with (solid) and without (dashed) including the PWNe contribution. Blue lines represent the total gamma fluxes predicted as a function of the energy for α = 1.8. The gray lines show a power-law with an index of 2.7 for comparison. 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 for R Φ = 250 (R Φ = 1500). 3 This clearly shows that this component is not negligible and cannot be ignored in the interpretation of Fermi-LAT diffuse emission data.
The effect of the unresolved TeV PWNe population on the determination of CR diffuse emission is discussed in Fig.1 where for illustrative purposes, we fix the transition energy E 0 = 0.3 TeV and the spectral index in the TeV domain β TeV = 2.3 and we only consider the effects of variations of the flux ratio R Φ . Black data points represent the total diffuse γ−ray flux observed by Fermi-LAT in each galactocentric ring given by Pothast et al. (2018) in 25 log-spaced energy bins between 0.34 − 228.65 GeV and in the latitude window |b| < 20.25 • . These data have been fitted with a single power-law ∝ E −Γ1 by Pothast et al. (2018), 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 two extreme cases R Φ = 250 and 1500. The red shaded area can be considered as the systematical uncertainty associated to the parameter R Φ . 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. 4 (see Eq. 7). 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 im-portant piece of information cannot be neglected and it is included for the first time in our work.
If we take unresolved PWNe emission into account, the evidence for CR spectral hardening in the inner Galaxy is 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 reported in the various panels of Fig.1 and are mildly dependent on the assumed value of R Φ . We see that they 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.
This statement does not critically depend on the assumed source spectrum, as can be verified by repeating our analysis for different combinations of the spectral parameters (R Φ , E 0 , β TeV ). The final results of our work are given in Tab. 2 and include systematic errors, evaluated as the maximal variation of Γ BF that is obtained when (R Φ , E 0 , β TeV ) are simultaneously varied in the 3-dim parameter space defined by the ranges The results of our work are compared with those given by other analyses in Fig. 2 where we show the γ-ray emissivity per H atom at 2 GeV (upper panel) which is a proxy of the CR spatial distribution in the Galaxy, and the CR proton spectral index (lower panel), 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. (2018) (red points) and Acero et al. (2016) (orange points). We also show with grey points the results obtained by Peron et al. (2021) by studying γ−ray emission in the direction of giant molecular clouds. The displayed error bars (for the black points) are obtained by summing in quadrature statistical and systematical uncertainties. We see that the emissivity calculated in this work is in good agreement with that obtained by Pothast et al. (2018) and 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. Notably, the three data-sets agree quite well with theoretical expectations for the CR dis-... tribution from GALPROP code (Vladimirov et al. 2011) (dashed blue line). 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. Table 2. Spectral indexes of the CR diffuse emission obtained by fitting the Fermi-LAT data with (ΓBF ) and without (Γ1) TeV PWNe unresolved contribution. The first error associated to ΓBF represents the systematic uncertainty (due to variations of RΦ, E0 and βTeV) while the second is the statistical one. The indexes Γ1 coincide with those obtained by Pothast et al. (2018). The TeV Galactic sky is dominated by a population of bright young PWNe whose properties are constrained by present HESS 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. 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.

Ring
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. (2019), between total γ-ray emission measured by HESS, Milagro, Argo and HAWC and that obtained by implementing CR spectral hardening.
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 (HESS). 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.
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 power-law 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 HESS. 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, suggested e.g. by Sudoh et al. (2019), 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 IC scattering of HE electron and positrons on background photons (CMB, starlight, infrared). In the Thompson regime, this naturally produces hard gammaray 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, see e.g. (Bucciantini et al. 2011;Torres et al. 2014;Sudoh et al. 2021). In the assumption of a broken power-law 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 through 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 do-mains. 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 = [2.3 − 2.4]. The spectral index β GeV of GeV emission is determined as a function of (R Φ , E 0 , β TeV ) by inverting Eq.(7) and it spans in the range β GeV = 1.18 − 2.00. 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. 3 of the Supplementary Material). 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 HESS 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.3). This is obtained by considering the best-fit spectra of each source, as given by Fermi-LAT and HESS 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 cumulative spectrum is well reproduced by a broken power-law with an energy break at about E 0 0.8 TeV, spectral indexes β TeV = 2.4 and β GeV = 1.89 and a R Φ 770 (red line in Fig.3). 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 ...
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 as recently measured by Diehl et al. (2006). If gamma-ray emission is powered by pulsar activity, the TeV-luminosity can be connected to the pulsar spin-down power, i.e.: where λ ≤ 1 and:Ė 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. (8) 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. (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. (8) that corresponds to a source luminosity function Y TeV (L TeV ) ∝ L −1.8 TeV . These two scenarios (α = 1.5 and α = 1.8) are considered as working hypotheses in this paper.
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 con- tribute 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. (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 4 . 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. The total flux produced by the considered population 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 Unresolved source contribution -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 can be evaluated as: where dN/dΦ GeV is the source flux distribution expected in GeV domain. This is connected to the source flux distribution in TeV domain by the R Φ parameter, according to: The latter can be calculated as a function of L TeV, Max and τ by using: where ρ(r) ≡ OW dΩ ρ(r, n). As it is discussed in Cataldo et al. (2020), the flux distribution scales as dN/dΦ TeV ∝ Φ −α TeV for Φ TeV → 0 and dN/dΦ TeV ∝ Φ −5/2 TeV for Φ TeV → ∞. This implies that Φ NR GeV ∝ R α−1 Φ , if L TeV, Max and τ are fixed, as it can be understood by using the low-flux scaling of dN/dΦ TeV .

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). Table 2. Spectral indexes of the CR diffuse emission for α = 1.5 obtained by fitting the Fermi-LAT data with (ΓBF ) and without (Γ1) TeV PWNe unresolved contribution. The two errors associated to our best fit values are the systematic one (due to variation of RΦ and E0) on the left and the statistical one on the right. The indexes Γ1 coincide with those obtained by Pothast et al. (2018).
In Fig. 5, we finally compare theoretical predictions from our population model with the cumulative distribution of HGPS sources (black solid line). The theoretical distribution for α = 1.8 (α = 1.5) is shown by a blue dashed line (red dashed line) and it is calculated by assuming that the maximal PWNe luminosity and spin-down timescale are L TeV, Max = 6.8 · 10 35 erg s −1 (L TeV, Max = 5.0 · 10 35 erg s −1 ) and τ = 0.5 · 10 3 y (τ = 1.7 · 10 3 y), respectively. These values have been obtained by performing an unbinned likelihood fit of 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 ).  Figure 5. 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).