The isotropic $\gamma$-ray emission above 100 GeV: where do very high energy $\gamma$ rays come from?

Astrophysical sources of very high energy (VHE; $>100$ GeV) $\gamma$ rays are rare, since GeV and TeV photons can be only emitted in extreme circumstances involving interactions of relativistic particles with local radiation and magnetic fields. In the context of the Fermi Large Area Telescope (LAT), only a few sources are known to be VHE emitters, where the largest fraction belongs to the rarest class of active galactic nuclei: the blazars. In this work, we explore Fermi-LAT data for energies $>100$ GeV and Galactic latitudes $b>|50^{\circ}|$ in order to probe the origin of the extragalactic isotropic $\gamma$-ray emission. Since the production of such VHE photons requires very specific astrophysical conditions, we would expect that the majority of the VHE photons from the isotropic $\gamma$-ray emission originate from blazars or other extreme objects like star-forming galaxies, $\gamma$-ray bursts, and radio galaxies, and that the detection of a single VHE photon at the adopted Galactic latitudes would be enough to unambiguously trace the presence of such a counterpart. Our results suggest that blazars are, by far, the dominant class of source above 100 GeV, although they account for only $22.8^{+4.5}_{-4.1}\%$ of the extragalactic VHE photons. The remaining $77^{+4.1}_{-4.5}\%$ of the VHE photons still have an unknown origin.


INTRODUCTION
Emitting photons with energies greater than 100 GeV is challenging. In nature, such photons are typically produced under extreme conditions surrounding compact astrophysical sources (Rieger et al. 2013), such as pulsars or black holes, where charged particles can be accelerated up to relativistic energies and hence emit very high energy (VHE, > 100 GeV) γ rays via inverse Compton (IC) scattering and/or curvature radiation (Sturrock 1971;Caraveo 2014;Blandford et al. 2019). Above ∼ 50 GeV, only a few hundred γ-ray sources are significantly detected (Wakely & Horan 2008;Ackermann et al. 2016) with the Fermi Large Area Telescope (LAT) and ground-based observatories, like the High Energy Stereoscopic System (Funk et al. 2004), the Very Energetic Radiation Imaging Telescope Array System (Weekes et al. 2002), and the Major Atmospheric Gamma Imaging Cherenkov Telescopes (Baixeras et al. 2003).
As listed in TeVCat 1 (Wakely & Horan 2008) and the Second Fermi -LAT Catalog of High-Energy Sources 2 (2FHL; Ackermann et al. 2016), the extragalactic VHE sky is dominated by blazars, which are radio-loud active galactic nuclei with a relativistic jet directed toward the observer (Urry & Padovani 1995;Blandford et al. 2019). Blazars of 105 • is applied, and we apply the standard good time interval filter DATA QUAL > 0 and the recommended instrument science configuration LAT CONFIG == 1.
For the main analysis, no exposure maps, livetime cubes or source maps are computed. We aim to analyze each event of our dataset individually. To estimate the Galactic contribution at such high latitudes and energies, we compare the predicted number of photons from four distinct Galactic models, namely gll iem v07, gll iem v06, gll iem v05 and gal 2yearp7v6 v0 in the energy range 100-500 GeV for both ROIs. This energy range has been chosen because the Galactic models prior to gll iem v07 do not reach TeV energies. For each model, we compute the total number of VHE photons by directly integrating them over the time, energy, and solid angle intervals adopted here. These Galactic models are compositions of several multiwavelength templates mapping the gas and dust in the Milky Way and differ in several aspects, as in the angular resolution of the gas distributions, and the gas and dust density/velocity profiles over the Galaxy 4 . We select the upper and lower values for the total number of photons estimated with these Galactic models as the systematic error band (in the 100-500 GeV energy range) and consider the predicted number of photons from gll iem v07 to be the most accurate value among all the models, since it is the latest model released to date. When we extrapolate these values for energies > 500 GeV (based on the energy distribution of photons in our dataset) and consider the Poissonian statistical error on the total number of photons, we get a total of 315 ± 39 (sys) ±18 (stat) and 239 ± 30 (sys) ±15 (stat) photons predicted by the Galactic model for the northern and southern ROIs, respectively.
Other major background contributors are the Fermi Bubbles. However, as shown in the next section, only the Southern bubble is a significant source of VHE photons on the analysis performed here and we remove it from our ROI. Our two ROIs are centered in the North and South Galactic poles and include 1165 and 794 VHE photons, respectively. This substantial difference in the number of photons is caused mainly by the uneven coverage of the sky performed by Fermi -LAT, which has a deeper exposure towards the Northern Galactic hemisphere (Ackermann et al. 2012b) and by the presence of Mkn 421 in the northern ROI.

METHODS
Several γ-ray blazars are relatively bright VHE sources. The VHE photons that they emit appear as clusters of events in the LAT data. In order to avoid to erroneously associate these clustered photons with multiple sources, we start our analysis by crossmatching the VHE photons with 4FGL-DR3 and masking those photon clusters coincident with 4FGL-DR3 sources, such that, to the aims of our analysis, they are represented by a single photon centered on Figure 2. Equivalent photons above 100 GeV in the northern and southern ROIs after masking for photon clusters coincident with 4FGL sources. Red dots represent equivalent photons coincident with the Southern Fermi Bubble and are excluded from the rest of the analysis. In principle, we assume that each blue dot can be associated to a blazar. the position of the γ-ray counterpart listed in 4FGL-DR3. From now on, we refer to our sample of photons after masking as "equivalent photons", since some of them may represent clusters of photons. In Figure 1 we show the effect of applying the mask to Mkn 421, where all VHE photons lying within a radius of 0.5 • from the γ-ray source center are substituted by a single equivalent photon. Adopting this strategy, we aim to be very conservative in excluding any clustering of VHE photons related to 4FGL-DR3: as can be seen in Figure 1, for strong VHE sources, like Mkn 421, a small percentage (∼ 3%) of VHE photons can be spread over a radius of 0.3 • ∼ 0.4 • . The total number of equivalent photons in the northern ROI drops to 743 after masking 58 clusters, while in the Southern ROI it drops to 656 after masking 41 clusters.
The spatial distribution of equivalent photons after masking photon clusters coincident with 4FGL-DR3 sources is shown in Figure 2, where the panels are centered at the North and South Galactic poles. Differently from the Northern Fermi Bubble, which barely reaches b = 50 • in Galactic latitude and does not contaminate our ROI, the Southern Fermi Bubble extends down to b = −54 • and is a significant source of background for our analysis. We therefore discard all equivalent photons spatially coincident with the bubble, as shown by the red dots in Figure 2. The chosen edges are b > −54 • , in Galactic latitudes, and −22 • < l < 16 • in Galactic longitude, in accordance with the bubble edges available in Sarkar (2019). The number of discarded equivalent photons is 44, reducing the total number of equivalent photons in the Southern ROI to 612. In Figure 3 we summarize all cuts applied to Fermi -LAT data for both ROIs.
After masking, we assume that each VHE equivalent photon is a potential γ-ray source and crossmatch them with several catalogs of astrophysical counterparts. Our main catalog is comprised of 5648 blazars collected from Roma-BZCat (Massaro et al. 2015), 4FGL-DR3, and identified in optical spectroscopic campaigns (Massaro et al. 2013;Paggi et al. 2014;Crespo et al. 2016;Marchesini et al. 2019;Peña-Herazo et al. 2019). Most of the blazars listed in this catalog have a confirmed optical spectrum (nearly 100% in the case of BL Lacs) or are confirmed γ-ray emitters. The impact of adopting other catalogs of γ-ray sources is discussed in §4.2. To test if these blazars indeed contribute to the VHE isotropic emission, we repeat the matching process for 5000 lists of mock VHE equivalent photons, generated by displacing the original positions of the VHE equivalent photons by a random value between 0 • and 5 • in a random direction of the sky (see §4). In order to define the optimal association radius, we first perform the crossmatches between the VHE equivalent photons and the sources listed in our main catalog considering a quite large constant association radius r = 0.5 • , and then compare the distribution of angular separation between the real and mock equivalent photon lists. From Figure  4, we define the optimal association radius to be r assoc ≈ 0.15 • ± 0.03 • , beyond which the associations are consistent with noise. This procedure has been successfully adopted in several works searching for the counterparts of γ-ray sources in radio, infrared and optical bands (more details can be found in Massaro et al. 2013Massaro et al. , 2014bGiroletti et al. 2016). We assume r assoc to be constant over the adopted energy range, which is reasonable, given that the Fermi -LAT containment radius for photons above 100 GeV is nearly constant (see section 2 in Abdo et al. 2011). If instead we use other catalogs to define the optimal association radius (as those discussed in §4.2), we find very similar values for r assoc , all around 0.15 • and consistent within the σ r error region.

RESULTS
The results for both ROIs are shown in Figure 5. For the northern ROI (hereinafter represented in blue), we find that 114 +1 −6 equivalent photons (after applying the masks previously discussed and considering the 1 σ error in r assoc ) are associated with 114 +1 −6 blazars. The expected number of random associations derived from our mock lists of VHE equivalent photons is 10 +5 −3 . For the southern ROI (hereinafter represented in orange), we have 69 +4 −2 equivalent photons associated with 69 +4 −2 blazars, with an expected number of random matches of 7 +2 −3 . In this figure, we fit the distributions of mock matches with a Poissonian function and the results indicate that blazars indeed account for a significant fraction of the VHE equivalent photons detected by Fermi -LAT, with an association level of confidence of 32.4σ in the North and 24.3σ in the South. However, most of the VHE equivalent photons still lack an astrophysical counterpart listed in the catalogs adopted in this work.
The basic statistics of our findings are shown in Table 1. We observe that 15.3 +1.5 −2.0 % of the VHE equivalent photons in the northern ROI are associated with blazars, while for the southern ROI this fraction is 11.3 +1.8 −1.5 %. Furthermore, only 10.5 +1.0 −1.1 % of all blazars scattered over both ROIs, as listed in our main catalog, have a corresponding VHE γ-ray counterpart detected with Fermi -LAT, most of which (∼ 70%) are BL Lac objects, although this class of source only makes 28% of our main catalog. These results tell us that 86.5 +1.6 −1.4 % of Fermi -LAT VHE equivalent photons at high Galactic latitudes are unlikely to be originated in blazars. Moreover, as the fraction of associated blazars is smaller than the fraction of associated equivalent photons, it is unlikely that the high number of non-associations is due to incompleteness in our blazar catalog (more details in §5).
As discussed in §2, the expected number of Galactic VHE photons for both ROIs is 315 ± 39 sys ± 18 stat in the North and 239 ± 30 sys ± 15 stat in the South, which correspond to 42 ± 8% total and 39 ± 7% total of all VHE equivalent photons . Distribution of angular separation between VHE equivalent photons and the blazars in our main catalog (orange line). The black line represents the average of the same distribution computed for 5000 mock equivalent photon lists and the gray shadow is the standard deviation of this distribution. We define the association radius as rassoc = 0.146 • , which is the median value for the intersection between the orange line and each one of the 5000 angular separation distributions for the mock lists of equivalent photons. The error in rassoc is taken as the standard deviation of these intersection points and has a value of 0.030 • . Throughout this work, we use the rounded value of the association radius rassoc = 0.15 • ± 0.03 • . Inset: zoom in the region around rassoc.  Table 1. Results from crossmatches. Columns represent the northern and southern ROIs individually and summed. The errors are estimated from the uncertainty in rassoc (i.e. computing the number of matches for rassoc, rassoc + σr, and rassoc − σr) and Galactic models. For the uncertainties in the fractions, we also take into account the binomial formula (also known as Wald method). We see that only 22.8 +4.5 −4.1 % of the extragalactic VHE equivalent photons have a blazar counterpart.

North
available in both ROIs, respectively. This means that after taking the Galaxy and the blazars into account, roughly 50% of the observed VHE equivalent photons still have unknown origin (assuming that none of the equivalent photons coincident with blazars are spurious or have Galactic origin). This is equivalent to saying that only 22.8 +4.5 −4.1 % of the extragalactic photons have a blazar counterpart and that the remaining 77 +4.1 −4.5 % of the extragalactic photons have unknown origin. A possible origin for a fraction of these photons, however, can be due to spurious signals induced by cosmic rays at the detector level (see §5). We fit each distribution with a Poissonian function, which ensures that blazars significantly contribute to the VHE sky at the 32.4σ level in the North and 24.3σ in the South (combined significance of 40.3σ).

Clustering of VHE equivalent photons
There are basically no clusters of VHE equivalent photons within an angular separation of 2 × r assoc . The exceptions are 8 pairs of equivalent photons in the southern and 15 pairs in the northern ROIs, and none of them match with the blazars in our catalog. These values are compatible with what is expected by chance (15 ± 5 for the north and 10 ± 4 for the south) when we crossmatch the lists of VHE equivalent photons from both ROIs with 5000 random distributions of VHE photons scattered over the same areas of the sky and with the same size of our original sample of equivalent VHE photons. In Figure 6, we show the distribution of angular separation between each pair of equivalent photons in our ROIs and their closest neighbors. We see that most VHE equivalent photons are truly isolated in the sky, being typically more than 1 • away from each other.

Repeating the analysis for other catalogs
We repeat the crossmatching process with other catalogs of candidate γ-ray emitters, again considering the whole sample of VHE equivalent photons and r assoc = 0.15 • ± 0.03 • as the association radius (see §3). For FRICAT and FRIICAT, two catalogs of Fanaroff-Riley (FR) radio galaxies (Capetti et al. 2017a,b), we find only 3 matches with FR I galaxies in the Northern ROI. Similarly, if we adopt a catalog containing only the non-blazars in 4FGL-DR3, we find only 8 matches in the North (three radio galaxies, one active galactic nucleus and three unidentified γ-ray sources, aka UGSs) and 7 in the South (one radio galaxy, one starburst galaxy and five UGSs).
When we use the combined sample of WIBRaLS and KDEBLLACS (hereinafter WISECATS), which are two catalogs of γ-ray blazar candidates (D'Abrusco et al. 2019; de Menezes et al. 2019) with 15121 sources in total (2468 sources in the Northern and 1929 in the Southern Galactic ROIs, i.e., far more sources than in our main catalog), the total number of matches is only 142 +15 −18 (87 +7 −11 equivalent photons in the North and 55 +8 −7 equivalent photons in the South), which is significantly smaller than what we have with our main catalog (see Figure 5). Out of these matches, 36 +2 −3 associations are not common to our main catalog. By taking these matches into account, the fraction of associated equivalent photons shown in Table 1 increases from 15.3 +1.5 −2.0 % to 18.0 +1.5 −1.9 % in the North and from 11.3 +1.8 −1.5 % to 13.9 +1.7 −1.3 % in the South (after subtracting the Galactic photons, the fraction of associated equivalent photons becomes 27.3 +4.8 −4.3 % over both ROIs, against the 22.8 +4.5 −4.1 % value reported in Table 1). The results for WISECATS are shown in Figure 7 together with the distribution of matches obtained from 5000 mock lists of equivalent photons. In order to properly compute these noise distributions, we first select WISECATS blazar candidates with b > |50 • | and, out of those, we randomly select 959 candidates in the North and 730 in the South, so that they match the number of blazars that we have for both ROIs in the main catalog (see Table 1). We then generate 10 lists of mock VHE equivalent photons (as discussed in §3) and crossmatch them with the randomly selected blazar candidates. After repeating this process 500 times, we obtain the final noise distributions with a total of 5000 lists of mock VHE equivalent photons. We therefore conclude that WISECATS, although presenting a higher completeness, do not present better results than our main catalog in spotting VHE blazar candidates. We also test other catalogs of possible γ-ray emitters, like the Catalogue of Extreme & High Synchrotron Peaked Blazars (3HSP; Chang et al. 2019) and the Candidate Gamma-Ray Blazar Survey Source Catalog (CGRaBS;Healey et al. 2008). For 3HSP we have a total of 129 +4 −6 matches for both ROIs, but only 11 ± 1 of them are not common to our main catalog, while for CGRaBS we find 19 ± 1 matches in the North and 10 ± 1 in the South, all of them being common to our main catalog. For the Chandra ACIS Survey for X-Ray AGN in Nearby Galaxies (CHANSNGCAT; She et al. 2017), we find only 6 ± 1 matches in the North and 3 +0 −1 in the South, while for the Australia Telescope National  von Kienlin et al. 2020), this time assuming a larger and fixed association radius of r assoc = 0.5 • instead of r assoc = 0.15 • ± 0.03 • . The arrival times of the VHE photons, however, are not consistent with the trigger times of the γ-ray bursts, meaning that a correlation is unlikely. Among these alternative catalogs, 3HSP and WISECATS lead the significance on the number of sources associated with VHE equivalent photons. We summarize the matches obtained with these alternative catalogs and the combined significance of both ROIs for each one of them in Table 2 . We see that the number of associations is smaller (87 +7 −11 in the North and 55 +8 −7 in the South) than those shown in Figure 5. We therefore conclude that WISECATS do not perform better than our main catalog in spotting VHE blazar candidates.

Catalog
Matches  Table 2. List of alternative catalogs and their total number of matches with the VHE equivalent photons considering rassoc = 0.15 • ± 0.03 • . The last column gives the combined significance of the matches for both ROIs. For the catalogs tagged with "*", the adopted association radius is fixed at rassoc = 0.5 • .
magnitudes available in these catalogs, but they do not show any peculiar property over the general populations of VERONCAT and MILLIQUAS sources. We also tested if the VHE equivalent photons match with a list of 35 active galactic nuclei known to host ultrafast outflows (see this list in Fermi-LAT Collaboration 2021), but we found only a couple of matches, leading to an association significance below 1.5σ. In a recent work by Roth et al. (2021), the authors found that most (or possibly all) of the isotropic γ-ray emission could be explained by star-forming galaxies. The authors also argue that for energies ∼ 1 TeV, the major contributors to the isotropic γ-ray emission are nearby (z ∼ 0.1) star-forming galaxies. To test this scenario, we collected 64 nearby (z max ∼ 0.06) star-forming galaxies from Ackermann et al. (2012c), which have unambiguous ongoing star formation, selected based on the presence of dense molecular clouds spotted by the HCN survey (Gao & Solomon 2004). Our results show that only a couple of these star-forming galaxies match with the VHE equivalent photons, giving an association significance of ∼ 1σ, and suggesting that nearby (z ≤ 0.06) star-forming galaxies (as selected here) may not significantly contribute to the isotropic γ-ray emission at energies > 100 GeV (assuming that the γ-ray emission of star-forming galaxies is correlated with their infrared emission collected in the HCN survey).

DISCUSSION AND CONCLUSIONS
For energies > 100 GeV and Galactic latitudes b > |50 • |, we would expect the extragalactic isotropic γ-ray emission to be dominated by blazars or other potential VHE sources, like star-forming and radio galaxies. It is unreasonable to assume that these VHE equivalent photons originate from Galactic VHE objects, like pulsar wind nebulae and supernova remnants, as out of 62 such objects listed in 4FGL-DR3, only two have b ∼ 32 • (both in the Large Magellanic Cloud), and the rest has b < |10 • | (Abdollahi et al. 2020(Abdollahi et al. , 2022. Above |50 • | in Galactic latitude, the majority of VHE γ-ray emitters must be extragalactic. In fact, in the third data release of 4FGL, only 1.4% of the γ-ray sources with b > |50 • | are associated with Galactic counterparts (i.e. 18 pulsars and one binary system), while 79% of them are associated with blazars.
After subtracting the expected number of Galactic VHE photons from our sample, we are left with ∼ 800 equivalent photons (North + South). In the most optimistic scenario, where all blazar associations are real, only 22.8 +4.5 −4.1 % of these equivalent photons can be associated with blazars, leaving 75% of the isotropic γ-ray emission unexplained. Given the incompleteness of our main catalog (described in §3), the real fraction of associations can be higher, however, if incompleteness is really the problem, in §4.2 we should observe more matches with WISECATS than with our main catalog, as this catalog has a much higher completeness for BL Lac objects, which are the main candidates for VHE emission (Massaro et al. 2013). Our results therefore tend to disfavor unknown blazars as the sources of the majority of VHE equivalent photons from the extragalactic isotropic γ-ray emission. Instead, we could speculate other origins for them, like for instance the scattering of the EBL by ultra-high energy cosmic rays or the scattering of soft photons in the Galactic halo that are not appropriately accounted for in the gll iem v07 model. The electromagnetic cascades resulting from such interactions could in principle create a truly diffuse isotropic γ-ray emission (Ackermann et al. 2015). Other possibilities for this emission could rise from the interaction of cosmic rays with the Solar radiation field (Orlando & Strong 2008) and Solar System debris (Moskalenko & Porter 2009), which are expected to contribute to the isotropic emission at some still unknown level. Furthermore, at the detector level, Fermi -LAT can eventually misinterpret a cosmic-ray induced event with a γ-ray, then introducing a spurious signal to the real data. We reduce the impact of this systematic problem by adopting the SOURCEVETO event class (see §2) in our analysis, which is one of the most stringent Fermi -LAT event reconstruction classes in discriminating cosmic rays from γ rays, even though it is not 100% effective. In the P8R3 SOURCEVETO Fermi -LAT data adopted in this work, the intensity of cosmic-ray induced events above 100 GeV is roughly constant, with a value of ∼ 10 −6 MeV cm −2 s −1 sr −1 (Fermi -LAT Collaboration, 2022, in preparation), which can account for ∼ 1% of the events at 100 GeV, ∼ 3% at 300 GeV, and less than 10% at 500 GeV. Since ∼ 90% of the photons in our sample (before applying the masks) have energies below 300 GeV, the final level of contamination by cosmic rays is << 10%. These levels of contamination do not significantly affect the fraction of blazar associations shown in Table  1 since, except for a handful of blazars, all matches in our analysis (before masking) have at least one VHE photon below 300 GeV, where the contamination is minimal.
Previous works on the isotropic γ-ray emission (see e.g. Harding & Abazajian 2012; Singal et al. 2012;Ajello et al. 2015) also found that blazars cannot account for all of the extragalactic isotropic γ-ray emission, although none of them used a method similar to the method adopted here. In Ackermann et al. (2015), the authors found that roughly half of the total extragalactic γ-ray emission at 100 GeV come from blazars which, in our work, would be equivalent to measuring the total number of photons in both ROIs (before masking) that are originated by blazars. For a matter of comparison, in our work the fraction of photons (before masking) coincident with blazars is 43.1%, in reasonable agreement with Ackermann et al. (2015), given the significant differences in our analyses. Furthermore, (Roth et al. 2021) recently found that a significant fraction of the isotropic γ-ray emission at ∼ 1 TeV can be attributed to the emission of nearby (z ≈ 0.1) star-forming galaxies using a theoretical approach. In this scenario, the VHE emission of star-forming galaxies is generated by cosmic rays accelerated in supernova remnants interacting with the interstellar medium. Our investigation of an observed list of low redshift (z ≤ 0.06) star-forming galaxies described in §4.2, disfavors the hypothesis that star-forming galaxies are the origin/counterparts of our list of VHE equivalent photons, at least in the very local Universe (z ≤ 0.06).
In this work we used Fermi -LAT observations to investigate whether the > 100 GeV extragalactic isotropic γ-ray emission above b = |50 • | can be explained by blazars. Our results suggest that blazars are not enough to explain this emission and are summarized below.
1. Only 15.3 +1.5 −2.0 % of all equivalent photons in the northern and 11.3 +1.8 −1.5 % in the southern ROIs are associated with blazars. After subtracting the Galactic photons from our sample, the fraction of extragalactic equivalent photons associated with blazars is 22.8 +4.5 −4.1 %. If additionally to the main catalog we also consider WISECATS, this fraction slightly increases to 27.3 +4.8 −4.3 %. This result suggests that the detection of a single VHE γ ray at b > |50 • | does not unambiguously lead to the detection of a blazar counterpart.
2. Based on the adopted catalogs, we found that 75% of the extragalactic isotropic γ-ray emission above 100 GeV and b > |50 • | has no clear origin.
3. Almost 70% of the matches are with BL Lac objects, which are sources known to typically present harder γ-ray spectra if compared to other blazars like Flat Spectrum Radio Quasars.
4. The tests we performed with 5000 mock catalogs ensure that the association of blazars with VHE equivalent photons is not by chance and that blazars do have a contribution, although small, in the VHE sky. It is unlikely that the observed small fraction of associations is due to the incompleteness of our catalog. The fact that we have more blazars than photons in both ROIs and also that we performed the same analysis with alternative, and more complete catalogs, rule out this possibility.
With the advent of the CTA (Actis et al. Acharya et al. 2013) in the upcoming years, new possibilities for investigating VHE sources are expected, helping in the quest of unveiling the origin of the extragalactic isotropic γ-ray emission.