Unresolved sources naturally contribute to PeV $\gamma$-ray diffuse emission observed by Tibet AS$\gamma$

The Tibet AS$\gamma$ experiment provided the first measurement of the total diffuse gamma-ray emission from the Galactic disk in the sub-PeV energy range. Based on the analysis of TeV sources included in the HGPS catalog, we predict the expected contribution of unresolved sources in the two angular windows of the Galactic plane observed by Tibet AS$\gamma$. We show that the sum of this additional diffuse component due to unresolved sources and the truly diffuse emission, due to cosmic ray interaction with the interstellar medium, well saturates the Tibet data, without the need to introduce a progressive hardening of the cosmic-ray spectrum toward the Galactic center.


INTRODUCTION
The Tibet ASγ collaboration has recently obtained a measurement of the Galactic diffuse γ-ray emission, providing the first evidence that this component extends to sub-PeV energies (Amenomori et al. 2021). The origin of this high-energy diffuse emission is extensively debated in the very recent literature, see e.g. Liu & Wang (2021); Koldobskiy et al. (2021); Fang & Murase (2021); Neronov et al. (2021).
An essential step for the interpretation of Tibet ASγ data is the evaluation of the cumulative flux produced by sources that are too faint to be individually resolved. These give rise to a large-scale diffuse emission, superimposed to that produced by cosmic ray (CR) interactions with the interstellar medium, that can affect the space and energy distribution of the observed signal. Population studies of TeV Galactic sources have shown that unresolved contribution is not negligible in the TeV domain, see e.g. Cataldo et al. (2020); Steppa & Egberts (2020). Recently Tibet ASγ (Amenomori et al. 2019) and HAWC (Abeysekara et al. 2020) have shown that several Galactic sources produce γ−rays above ∼ 50 TeV. Moreover, LHAASO-KM2A reports the detection of more than 530 photons at energies above 100 TeV and up to 1.4 PeV from 12 ultra-high-energy Corresponding author: Giulia Pagliaroli giulia.pagliaroli@gssi.it γ-ray sources (Cao et al. 2021). It is thus natural to expect that unresolved source contribution may be relevant also in the sub-PeV energy range.
In this paper, we took advantage of the observations provided by H.E.S.S. Galactic Plane Survey (HGPS) in the 1 − 100 TeV energy domain to discuss the implications of unresolved pulsar-powered sources for the interpretation of Tibet ASγ results. Our approach is based on standard assumptions for the spatial and intrinsic luminosity source distributions (also employed e.g. Steppa & Egberts (2020); Strong (2007); Pothast et al. (2018)). These prescriptions are physically well-motivated for a population of sources powered by pulsar activity, like Pulsar Wind Nebulae (PWNe) or TeV halos which could represent the dominant component of the TeV sky, as recently suggested by Sudoh et al. (2019); Abdalla et al. (2018a). We show that the inclusion of unresolved contribution from these sources, adding up to standard predictions for CR diffuse flux (Lipari & Vernetto 2018), naturally explains the Tibet ASγ data, without the need to introduce a progressive hardening of the CR spectrum toward the Galactic center.

SOURCE POPULATION
It is naturally expected that sources provide a relevant contribution to high energy γ-ray Galactic emission, see e.g. Ahlers & Murase (2014) for an estimate of the total emission expected from different classes of Galactic sources. However, the key point is to quantify the fraction of sources that are not resolved by experiments, since these generate a cumulative diffuse emission that adds up to the one produced by CR interactions. The goal of this paper is to predict the unresolved pulsarpowered sources contribution to galactic γ-ray emission in TeV and sub-PeV energy domain. In order to estimate this quantity, we consider the information provided by HGPS in the energy range 1 TeV ≤ E γ ≤ 100 TeV.
We assume that average source emission spectrum can be described by a power-law with exponential cutoff, i.e.
The index β is given by the average spectral index of sources included in the HGPS catalog, i.e. β = 2.3. The cutoff energy is chosen as E cut = 500 TeV. This value is not yet constrained and it corresponds to assuming that TeV source spectrum observed by H.E.S.S. can be extrapolated with a moderate suppression in the sub-PeV region. In this respect, it is interesting to note that the Tibet ASγ experiment has recently shown that Crab Nebula emits γ−rays in the sub-PeV region with 5.6σ statistical significance (Amenomori et al. 2019). Moreover, HAWC experiment has reported evidence of several Galactic sources emitting above ∼ 50 TeV (Abeysekara et al. 2020) whose γray flux can be described as leptonic emission from electrons/positrons injected with power-law spectrum and exponential cutoff around 1 PeV (Sudoh et al. 2021). We remark that the adopted value of E cut moderately affects flux predictions for the low energy data points given by Tibet ASγ that are relatively close to the range probed by H.E.S.S.. At larger energies (E γ ≥ 1 PeV), the source emission spectrum is limited by photons absorption in the interstellar radiation field (mainly due to Cosmic Microwave Background radiation) that suppresses the flux produced by distant sources. We take this into account, as it is described in Vernetto & Lipari (2016). The sources space and luminosity distribution is described by: where r indicates the distance from the Galactic Center. The function ρ(r) is assumed to be proportional to the pulsar distribution parameterized by Lorimer et al. (2006) and to scale as exp (− |z| /H) with H = 0.2 kpc, along the direction z perpendicular to the Galactic plane. It is conventionally normalized to one when integrated in the entire Galaxy. The function Y TeV (L TeV ) gives the source intrinsic luminosity distribution in the TeV energy domain 1 . It is parameterized as a power-law: in the luminosity range L TeV, Min ≤ L TeV ≤ L TeV, Max . This distribution is naturally obtained for a population of fading sources, such as PWNe or TeV Halos, created at a constant rate R and having 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 birth rate of PWNe or TeV Halos is similar to that of SN explosions in our Galaxy, i.e. R R SN = 0.019 yr −1 (Diehl et al. 2006). Since γ−ray emission is powered by pulsar activity, the TeV-luminosity L TeV is a fraction of the spin-down powerĖ, i.e. we assume L TeV = λĖ where λ ≤ 1. If pulsar energy loss is dominated by magnetic dipole radiation (braking index n = 3) and the efficiency of TeV emission does not depend on time (λ ∼ const), the exponent in Eq. (3) is γ = 2, that corresponds to a source luminosity function 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 2 . In this case, one obtains γ 1.2 in Eq.

THE UNRESOLVED SOURCE CONTRIBUTION
In a recent paper (Cataldo et al. 2020), we have shown that H.E.S.S. observations can be used to efficiently constrain the TeV sources population. In particular, we considered HGPS sources producing a photon flux above 1 TeV larger than 10% of the CRAB flux. Above this threshold, the HGPS catalog can be considered complete (Abdalla et al. 2018b) and it includes 32 sources. We removed from the analysis sources firmly associated with SNRs (i.e. Vela Junior, RCW 86, RX J1713.7-3946) and we treated the residual 29 sources as pulsarpowered. This assumption is justified being this subset composed by 8 firmly identified PWNe, 2 composite objects, showing evidence of both shell and nebular emission, and 19 unidentified sources. In the unidentified sub-sample, 12 sources have been considered as candidate PWNe in further studies on the basis of new data and/or phenomenological considerations (Sudoh et al. (2021), Wakely & Horan (2008), Abdalla et al. (2018a), Giacinti et al. (2020)). The average spectral index of the considered sample of 29 HGPS sources is 2.34 supporting our assumption β = 2.3.
Namely, by fitting the flux, latitude and longitude distribution of bright sources in the HGPS catalog, we have obtained L TeV, Max = 4.9 +3.0 −2.1 × 10 35 erg s −1 and τ = 1.8 +1.5 −0.6 × 10 3 yr for α = 1.5 (L TeV, Max = 6.8 +6.7 −3.0 × 10 35 erg s −1 and τ = 0.5 +0.4 −0.2 × 10 3 yr for α = 1.8). These results are also valid for extended TeV sources, provided that they have dimensions that do not exceed ∼ 40 pc. Moreover, our results are consistent with those obtained by a completely independent analysis of HGPS catalog performed by Steppa & Egberts (2020). In this paper, we take advantage of the results of Cataldo et al. (2020) to estimate the unresolved flux produced by the considered population in the TeV and sub-PeV energy domain. For definiteness, we take as a reference the case α = 1.5. Slightly larger fluxes are obtained for α = 1.8.
In order to evaluate the cumulative contribution to diffuse γ−ray signal of sources which are too faint to be individually detected, we introduce a flux detection threshold Φ th TeV whose value is estimated by considering the performances of H.E.S.S. detector. H.E.S.S. is able to resolve point-like sources if they produce an integrated flux Φ TeV in the [1 − 100] TeV energy domain that is larger than 0.01Φ CRAB , where Φ CRAB = 2.26 × 10 −11 cm −2 s −1 is the flux produced by CRAB in the same energy range. Extended sources can however escape detection even if their flux exceeds this value, as it is e.g. understood by looking the sensitivity curve presented in Fig. 13 of Abdalla et al. (2018b). The HGPS catalog can be considered complete for objects producing fluxes larger than 0.1Φ CRAB (with the exception of sources having an angular extension larger than ∼ 1 • which cannot be observed by H.E.S.S). Taking this into account, we calculate the cumulative emission of unresolved sources as: where ϕ(E γ ) is the average source emission spectrum 3 : the quantity dN/dΦ TeV is the source flux distribution in a given region of the sky (see Eq. A.6 of Cataldo et al. (2020)). The function η(E γ , Φ th TeV ) is the average survival probability of photons with energy E γ emitted by unresolved sources. It includes the gamma absorption due to gamma-gamma interactions on CMB photons following the numerical approach described in Lipari & Vernetto (2018). We vary the flux detection threshold in the range Φ th TeV = [0.01 − 0.1]Φ CRAB . Note that the diffuse Galactic γ-ray flux at sub-PeV energies measured by Tibet ASγ is obtained by subtracting/masking the contribution of sources which are included in the TeVCAT catalog (Wakely & Horan 2008). This implies that sources should be faint not only at sub-PeV energies (but also at TeV energies) to escape detection. As a consequence, the above approach which is based on the detection capabilities of experiments operating in the TeV domain is also adequate to investigate unresolved source contribution in the sub-PeV energy range.

RESULTS
The flux produced by faint sources which are not individually detected adds up the CR diffuse emission, shaping the radial and spectral behaviour of the total diffuse γ-ray flux observed by different experiments. In Fig. 1, we shows the theoretical predictions for the total diffuse γ−ray flux (green band) as a function of energy in the regions of the Sky probed by Tibet ASγ experiment (Amenomori et al. 2021). The upper panel (a) refers to the region 25 • < l < 100 • , while the lower panel (b) shows the region 50 • < l < 200 • ; both of them correspond to the latitude range |b| < 5 • . The unresolved sources contribution is obtained as described in the previous section and the thickness of the darker green band corresponds to the uncertainty in flux detection threshold Φ th TeV . Namely, the upper and lower green lines are obtained by assuming Φ th TeV = 0.1 Φ CRAB and Φ th TeV = 0.01 Φ CRAB , respectively. The light green band also includes the uncertainty on this prediction due to the correlated variations of the source population parameters τ and L TeV, Max within their 1σ uncertainty. The truly diffuse emission, produced by CR interactions with the interstellar gas, is shown by grey solid lines in Fig.1 and corresponds to the "space-independent" model of Lipari & Vernetto (2018). Red data points show the diffuse flux measured by Tibet ASγ. These are obtained after subtracting events within 0.5 • from known TeV sources included in the TeVCAT catalog (Wakely & Horan 2008). The error bars show 1σ statistical errors. Finally, we also display the CR diffuse flux corresponding to the "space-dependent" model of 50°< l < 200°, |b| < 5°F igure 1. Differential energy spectra of diffuse γ−rays from the Galactic plane in two different angular regions. Red data points are the measurements provided by Tibet (Amenomori et al. 2021). Blue data points in the upper panel are provided by Argo-YBJ (Bartoli et al. 2015), while blue triangles in the lower panel are upper limits by the CASA-MIA experiment (Chantell et al. 1997). Solid and dashed curves display the predicted energy spectra by the space-independent and space-dependent models by Lipari & Vernetto (2018), respectively. The green shaded band represents the total diffuse γ−ray emission obtained by adding the unresolved source contribution estimated in this paper to the γ-ray truly diffuse emission from space-independent model by Lipari & Vernetto (2018). Lipari & Vernetto (2018) (gray dashed lines) to permit comparison with our predictions. This is obtained by assuming that the CR spectrum in the inner Galaxy is harder than at the Sun position, as it seems to be suggested by Fermi-LAT data (Acero et al. 2016;Yang et al. 2016;Pothast et al. 2018) (see Vecchiotti et al. (2022), Peron et al. (2021) Pagliaroli et al. (2016) similar to that adopted by Lipari & Vernetto (2018). Fig. 1. First, we see that unresolved source contribution is not negligible for E γ ≥ 1 TeV and becomes progressively more relevant as energy increases as it was also pointed out in Linden & Buckman (2018). This is a natural consequence of the fact that sources are expected to have, on average, harder spectrum (a part from cutoff effects) than CR diffuse emission. At the energy E γ 150 TeV, corresponding to the first data point of Tibet ASγ, the cumulative flux produced by faint sources is estimated to be 49%−154% (25%−79%) of the truly diffuse signal in the region 25 • < l < 100 • (50 • < l < 200 • ).

Several important conclusions can be obtained from
It is useful to investigate the stability of our results concerning the possibility that some of the 29 sources considered in our analysis turn out to be not pulsarpowered. For this purpose, we repeat the fit of the flux, latitude and longitude distribution of HGPS sources taking into account only the 22 sources with a clear or potential association to PWNe. The new best-fit values with this reduced sample are L TeV, Max = 6.2 × 10 35 erg s −1 and τ = 1.1 × 10 3 yr. The unresolved flux due to the PWNe population obviously decreases, however, the prediction still falls inside the 1σ statistical uncertainty band reported in light green in Fig. 1.
We remark that our prediction is only marginally dependent on the high-energy extrapolation of the source spectra (e.g. the adopted value of the cutoff energy) since the lowest energy data point of Tibet ASγ experiment nearly overlaps with the energy range probed by H.E.S.S.. This is evident in Fig.2 where green lines represent the total diffuse γ-ray emission for intermediate sensitivity threshold Φ th TeV ∼ 0.05 Φ CRAB , obtained by assuming different values for the spectral energy cut-off from 1 PeV to 100 TeV. In particular, if we assume a lower energy cutoff E cut = 300 TeV (E cut = 100 TeV), the unresolved source flux at 150 TeV decreases by ∼ 17% (∼ 68%) with respect to the reference case E cut = 500 TeV. We can also consider the effects of source spectral index variations (not shown in Fig. 2 to avoid overcrowding it). If we take β = 2.4 (β = 2.5 TeV), the unresolved source flux at 150 TeV decreases by ∼ 32% (∼ 54%) with respect to the reference assumption β = 2.3. In all these cases, our predictions for the total diffuse γ-ray emission are still consistent with the first two Tibet data points.
It is interesting to investigate the typical age of sources that give relevant contribution to unresolved signal. In Fig.3 we show the relative contribution to unresolved emission as a function of log 10 (t/1kyr) where t is the age of the source, i.e. the time passed since pulsar formation. Thick lines refer to the sky region |b| < 5 • and 25 • < 50°< l < 200°, |b| < 5°F igure 2. We highlight the role of the assumed energy cut for the sources spectra. The green lines show the differential energy spectra of diffuse γ−rays from the Galactic plane in two different angular regions for an intermediate sensitivity threshold of Φ th TeV ∼ 0.05 ΦCRAB and assuming different values for the spectral energy cut-off from 1 PeV to 100 TeV. l < 100 • while dashed ones refer to the most lateral region between 50 • < l < 200 • . With blue lines we show the case of detection threshold Φ th TeV = 0.01Φ CRAB , corresponding to the lower bound of the green band for the total diffuse flux reported in Fig.3. In this case the dominant contribution is provided by PWNe with an age ranging between t ∼ (22 − 33) kyr, depending on the sky region considered. Red lines correspond to the upper bound of the green band for the total diffuse flux reported in Fig.1 with Φ th TeV = 0.1Φ CRAB . In this case a younger population provides a dominant contribution peaking between t ∼ (7 − 11) kyr, depending on the sky region considered. The unresolved flux contribution due to sources older than ∼ 100 kyr (likely TeV Halos) is expected to be at most the 20%.
It is important to remark that our calculations naturally reproduce the Tibet ASγ results both in the low and high longitude observation window, corroborating the conclusion that unresolved PWNe provide a relevant contribution in this energy range. A test of our calcu-  Figure 3. The distribution of sources contributing to unresolved signal as a function of log 10 (t/1 kyr). Thick lines refer to the sky region |b| < 5 • and 25 • < l < 100 • while dashed ones refer to the region between 50 • < l < 200 • . The blue (red) lines are obtained for Φ th TeV = 0.01ΦCRAB (Φ th TeV = 0.1ΦCRAB).
lations could be provided by additional observations in the TeV and PeV energy domain. Future experiments with improved sensitivity will be e.g. able to resolve more sources in the regions of interest. As an example, CTA, with a sensitivity according to Sudoh et al. (2019), should be able to resolve about 280 (140) pulsar-powered sources in the whole Galaxy, if the typical source size is 10 pc (40 pc). We also note that LHAASO-KM2A experiment recently presented preliminary determinations of the diffuse γ−ray signal in the sub-PeV energy domain (Zhao et al. 2021). The LHAASO-KM2A diffuse flux measurement is obtained by masking sources included in TeV-CAT catalog (Wakely & Horan 2008) with a relatively large mask radius 5 . LHAASO-KM2A preliminary data for the longitude window 25 • ≤ l ≤ 100 • are lower than Tibet ASγ results and are close to predictions for the "space-dependent" CR-diffuse emission of Lipari & Vernetto (2018). They cannot be however directly compared with theoretical determinations of the diffuse signal because the masking procedure cuts out a large part of the galactic plane, as it is seen by Fig. 2 of Zhao et al. (2021), where most of the CR-diffuse signal and unresolved emission is produced. They thus naturally represent a lower limit to the total diffuse emission.
Finally, by comparing our predictions in panel a) with the truly-diffuse flux given by the "space-dependent model" of Lipari & Vernetto (2018) (grey dashed lines), we see that the unresolved sources can mimic, at rela-5 The adopted mask radius is R = 2 σ 2 ext + p.s.f. 2 , where σext represents the extension of the source while p.s.f. is LHAASO-KM2A point-spread-function, and it is typically comparable to or larger than 1 • (Zhao et al. 2021). tively low galactic longitudes, the effects produced by CR-spectral hardening in the inner Galaxy (and viceversa). This is in agreement with what we have found in Vecchiotti et al. (2022). By looking at panel b), we note, however, that the previous statement is not valid in the longitude range 50 • ≤ l ≤ 200 • where unresolved sources, differently from CR-spectral hardening, produce an enhancement of the diffuse emission and provide a good description of the Tibet ASγ data.
The different behaviour has a direct explanation. The diffuse emission at l ≥ 60 • is mainly produced, for geometrical reasons, at galactocentric distances comparable to or larger than the distance of the Sun from the Galactic center. As a consequence, it is not sensitive to variations of the CR distribution in the inner Galaxy. On the contrary, it can be increased by unresolved sources, which are distributed in the whole Galaxy, and thus produce a non vanishing contribution even at large latitudes. We remark that that this point is important because it provides the possibility to distinguish between the two effects in present and future experiments. At the moment, the inclusion of unresolved PWNe contribution produces a better description of the Tibet ASγ data than CR spectral hardening.