{\it Fermi}-LAT Stacking Analysis Technique: An Application to Extreme Blazars and Prospects for their CTA Detection

TL;DR

This paper introduces a likelihood profile stacking method using Fermi-LAT data to analyze undetected extreme blazars, revealing their potential TeV emission and contribution to the extragalactic gamma-ray background, with implications for CTA observations.

Contribution

The paper develops a novel stacking technique for Fermi-LAT data to study gamma-ray undetected sources, providing new insights into their spectra and detectability with CTA.

Findings

Detected a significant cumulative gamma-ray signal from undetected extreme blazars.

Found that their intrinsic spectra likely extend above 1 TeV without softening.

Estimated that ~10% of the extragalactic gamma-ray background at 100 GeV originates from these sources.

Abstract

We present a likelihood profile stacking technique based on the {\it Fermi}-Large Area Telescope (LAT) data to explore the $\gamma$-ray characteristics of {\it Fermi}-LAT undetected astrophysical populations. The pipeline is applied to a sample of $\gamma$-ray unresolved extreme blazars, i.e., sources with the highest synchrotron peak frequencies ($\nu_{\rm Syn}^{\rm peak}\geqslant 10^{17}$ Hz), and we report a cumulative $\gamma$-ray detection with more than 32$\sigma$ confidence for 2 degrees of freedom. Comparing the generated stacked $\gamma$-ray spectrum with the sensitivity limits of the upcoming Cherenkov Telescope Array (CTA), we find that the {\it Fermi}-LAT undetected population of such extreme blazars, on average, may remain well below the CTA detection threshold due to their faintness and extragalactic background light (EBL) absorption. However, $\gamma$-ray detected blazars belonging to the same class are promising candidates for CTA observations. The EBL corrected stacked spectra of these sources do not show any softening up to 1 TeV. This finding suggests the inverse Compton peak of extreme blazars to lie above 1 TeV, thus indicating a hard intrinsic TeV spectrum. Our analysis also predicts that at 100 GeV, at least $\sim$10\% of the diffuse extragalactic $\gamma$-ray background originates from the $\gamma$-ray undetected extreme blazars. These results highlight the effectiveness of the developed stacking technique to explore the uncharted territory of $\gamma$-ray undetected astrophysical objects.

Full text

methods: data analysis — gamma rays: general — galaxies: active — galaxies: jets — BL Lacertae objects: general

††slugcomment: ApJ Letters accepted††software: fermiPy (Wood et al., 2017)

1 Introduction

The Large Area Telescope (LAT) onboard the Fermi Gamma ray Space Telescope has revealed various types of astrophysical objects as $\gamma$-ray emitters. The most numerous of them are blazars, i.e., radio-loud quasars with powerful relativistic jets pointed towards the observer, followed by narrow-line Seyfert 1 galaxies, pulsars, and many others (cf. The Fermi-LAT collaboration, 2019b). Concerning blazars, the Fermi-LAT observations have allowed us to explore various unsolved problems related to jet physics and/or $\gamma$-ray astronomy in general. A few examples are: the cosmic evolution of blazar jets (Ajello et al., 2012, 2014), connection of the central engine (i.e., central black hole and the accretion disk) with the relativistic jet (e.g., Petropoulou & Dermer, 2016; Paliya et al., 2017), the measurement of the Extragalactic Background Light (EBL, Domínguez & Ajello, 2015; Abdollahi et al., 2018), and the contribution of $\gamma$-ray blazars to the Extragalactic Gamma-ray Background (EGB, Ackermann et al., 2015; Ajello et al., 2015).

There are probably astrophysical source populations that are yet to be detected in the $\gamma$-ray band. They are, similar to the Fermi-LAT detected objects, crucial to study their contribution to the diffuse $\gamma$-ray background emission. Focusing on blazars, a characterization of the $\gamma$-ray properties is also pivotal to determine their detectability with the next generation high-energy missions, e.g., the All-sky Medium Energy Gamma-ray Observatory (AMEGO, Moiseev & Amego Team, 2017) and the Cherenkov Telescope Array (CTA, Actis et al., 2011). In particular, AMEGO (with an energy range 200 keV to $\gtrsim$10 GeV) is expected to discover some of the most powerful blazars, especially at high-redshifts ($z>3$, Paliya et al., 2019). On the other hand, CTA, which will operate in the $\sim$0.02$-$300 TeV band, will observe the most efficient particle accelerator jets from BL Lac objects, a sub-class of blazars with no or weak optical emission lines (Stickel et al., 1991), along with other types of sources, e.g. flat-spectrum radio quasars. Until then, the high-energy properties of these peculiar objects can be explored using Fermi-LAT observations.

A useful methodology to explore the characteristics of any astrophysical populations, especially undetected ones, is the stacking technique. This has been successfully applied earlier to Energetic Gamma-Ray Experiment Telescope data to search for the $\gamma$-ray signal from, e.g., cluster of galaxies (Reimer et al., 2003) and infrared (IR) galaxies (Cillis et al., 2005). Also in the Fermi-LAT era, various stacking algorithms have been developed to search for $\gamma$-ray emission from undetected populations. Huber et al. (2012) proposed a method which co-adds the Fermi-LAT count maps and performs a maximum likelihood analysis on the combined data to derive the $\gamma$-ray parameters. The Fermi-LAT data analysis software, Fermitools provides a package, Composite2, which makes use of summed log-likelihood functions to combine the likelihood fitting of multiple sources at once. This tool ties together the spectral parameters of interest for all sources under consideration before performing the fit and has been used in some studies (e.g., Ackermann et al., 2011, 2014).

Here we present a new approach based on the stacking of the individual source likelihood profiles to explore the $\gamma$-ray properties of the Fermi-LAT undetected population. The developed technique is sensitive to extract the faint $\gamma$-ray signal, is quick, and is flexible to be used for any kind (binned or unbinned) of Fermi-LAT data analysis. It can also be applied to $\gamma$-ray detected sources to estimate the average properties of the population. Unlike Composite 2, the tool can be used to independently generate the likelihood profiles of as many sources as one needs before combining them, at the expense of parallel processing computational resources. We describe the steps of the stacking technique in Section 2 and present the results of its application to a sample of extreme blazars in Section 3. The results associated with the validation of the stacking technique are presented in Section 4. In Section 5, we discuss and summarize our findings. All the quoted uncertainties are estimated at 1$\sigma$ level, unless specified.

2 The Pipeline

The input for the stacking pipeline is the list of sky positions, i.e., right ascension and declination, for the objects under consideration. The tool works in the following two steps:

2.1 Pre-processing

First, a standard likelihood analysis is performed on the Fermi-LAT data covering a given time period, energy range, and a selected region of interest (ROI) using fermiPy (Wood et al., 2017). To model the $\gamma$-ray sky, we consider the recently released fourth catalog of Fermi-LAT detected sources (4FGL, The Fermi-LAT collaboration, 2019a). We use the latest interstellar emission model and the standard template for the isotropic emission111https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html. All 4FGL sources lying within ROI size+$R$ are considered in the likelihood fitting, where size of the ROI and additional radius $R$ depend on the minimum energy chosen for the analysis. Spectral parameters associated with the sources lying within the ROI are allowed to vary during the likelihood fitting and are kept fixed to the 4FGL values if they lie outside the ROI. After a first round of the optimization, we scan the ROI to search for unmodeled $\gamma$-ray objects by generating test statistic (TS) maps. The maximum likelihood TS is defined as TS = $2\log(\mathcal{L}_{1}-\mathcal{L}_{0}$), where $\mathcal{L}_{0}$ and $\mathcal{L}_{1}$ denote the likelihood values without (i.e., null hypothesis) and with (alternative hypothesis) a point source at the position of interest, respectively (Mattox et al., 1996). Since unmodeled sources can have hard or soft spectra, TS maps are generated for various photon indices, e.g., 1.5, 2, 2.5. When an excess emission (TS$>$25) is found, it is added to the sky model following a power law spectral model and a second set of TS maps is generated. Once all excesses above the background are identified and inserted to the model, a final likelihood fit is performed to optimize all the free parameters in the ROI and derive the spectral parameters for the source of interest. This exercise enables us to segregate the whole sample in the $\gamma$-ray detected and undetected objects.

2.2 Stacking

The pipeline considers the source lists generated in the previous step and proceeds as follows:

Assuming that the average spectral behavior of the source population is well represented with a power law model, we generate a grid of photon flux ($10^{-15}-10^{-8}$ ${\rm\,ph\,\,cm^{-2}\,s^{-1}}$ in 50 logarithmic steps) and photon index (1.5$-$3.5 in the interval of 0.1). Note that the lower limit of the photon flux grid should be small enough to represent the absence of any $\gamma$-ray source so that the likelihood value computed at this photon flux refers to the null hypothesis of the $\gamma$-ray detection. 2. 2.

The photon flux and index at a given grid point are then used as the source model parameters and the fitting is performed to determine the log-likelihood value. This step is repeated at every grid point, thus effectively generating a likelihood profile. To speed-up the process, we freeze the spectral parameters of all other modeled sources to the optimized values derived in the previous step, except the diffuse background models which are allowed to vary. By subtracting the log-likelihood value at the lowest flux (representing the null hypothesis) from the generated profile, we compute the TS or detection significance profile for a given source. This step is repeated for all $\gamma$-ray undetected sources in the sample and a set of TS profiles is created. 3. 3.

Since the log-likelihood is additive in nature, we add the generated TS profiles to produce a combined profile representing the whole sample. The TS peak position and 1$\sigma$ uncertainties in the associated spectral parameters are then estimated by fitting a spline function. The photon flux and index values associated with the TS peak represent the average spectral parameters of the whole population.

3 Application to Extreme Blazars

We apply the developed stacking technique to a sample of extreme blazars and use them to demonstrate the robustness of the pipeline by performing simulations and background checks.

Extreme blazars are a sub-sample of BL Lac objects that have the synchrotron peak222The spectral energy distribution (SED) of a blazar is characterized by a double hump structure. The low-energy peak, usually located in the sub-mm-to-X-ray bands, originates from synchrotron emission, whereas, the high-energy peak is usually explained by inverse Compton scattering off low energy photons by the relativistic electrons present in the jet (see, e.g., Böttcher, 2019, for a review). located at very high frequencies ($\nu_{\rm Syn}^{\rm peak}\geqslant 10^{17}$ Hz) indicating them to host some of the most efficient particle accelerator jets (e.g., Costamante et al., 2001). A high synchrotron peak frequency also suggests their inverse Compton peak to be located at very high energies (VHE, $>$100 GeV) making them bright TeV candidates and a promising source population for CTA detection. However, the same phenomenon causes the extreme blazars to be faint and less variable in the Fermi-LAT energy range (e.g., Tavecchio et al., 2011; Arsioli et al., 2018). Therefore, it is instructive to perform a stacking analysis of the $\gamma$-ray undetected extreme blazars and determine whether they are GeV $\gamma$-ray emitters as a whole. The derived spectral parameters can be used to explore the detectability of extreme blazars with CTA and also use them to probe the EBL and the extragalactic VHE background (Padovani et al., 1993; Ajello et al., 2015).

We select 337 extreme blazars with known redshift from a sample of 2011 high-synchrotron peaked (HSP, $\nu_{\rm Syn}^{\rm peak}\geqslant 10^{15}$ Hz, Abdo et al., 2010) BL Lac objects included in the third catalog of HSP blazars (3HSP, Chang et al. 2019, see also Chang et al. 2017)333http://www.ssdc.asi.it/3hsp/. The sample has a redshift range of 0.01$-$0.85 with a mean redshift of 0.35. We use the P8R3 LAT data acquired in $\sim$127 months (2008 August 4 to 2019 March 18) of Fermi operation and in the energy range of 10$-$1000 GeV for the analysis. The minimum energy is chosen as 10 GeV to allow a maximum overlap with the frequency band covered by CTA. Furthermore, the Fermi-LAT PSF considerably improves above 1 GeV (e.g., Atwood et al., 2013) thereby enabling a better source localization and suppressing of the diffuse background emission which is bright at MeV energies (cf. Acero et al., 2016). We use a squared ROI of $10^{\circ}\times 10^{\circ}$ and adopt a zenith angle cut ($z_{\rm max}$) of 105*∘*. With these settings, the Pre-processing led to a significant $\gamma$-ray detection of 165 extreme blazars444Among 165 $\gamma$-ray detected objects, 154 are also present in the 4FGL catalog.. The TS distribution of the remaining 172 $\gamma$-ray undetected extreme blazars is shown in the left panel of Figure 1. This is consistent with the null hypothesis or background fluctuations, i.e., $\chi^{2}$ distribution with 2 degrees of freedom, at low TS and shows excess above TS$\gtrsim$10 that can be attributed to the jet emission.

The stacking pipeline is then applied to the 172 extreme blazars and we show the results in the right panel of Figure 1. As can be seen, the stacked emission is well detected by the Fermi-LAT with a TS = 1062 ($\sim 32\sigma$ significance for 2 degrees of freedom). The associated peak $\gamma$-ray flux and photon index are $F_{\rm 10-1000~{}GeV}=6.51^{+0.36}_{-0.35}\times 10^{-12}$ ${\rm\,ph\,\,cm^{-2}\,s^{-1}}$ and $\Gamma_{\rm 10-1000~{}GeV}=2.08^{+0.07}_{-0.06}$. These results demonstrate the capabilities of the stacking technique to extract the faint $\gamma$-ray signal from a Fermi-LAT undetected population.

The strong $\gamma$-ray signal has also allowed us to perform the stacking in smaller energy bins and generate a stacked $\gamma$-ray spectrum. Using the best-fit photon flux and index derived in each of the seven energy bins, we make the bow-tie plot and show the cumulative spectrum in the left panel of Figure 2. In this diagram, we also show the stacked $\gamma$-ray spectrum of 165 Fermi-LAT detected extreme blazars generated using the same pipeline. For a comparison, we overplot sensitivity limits of CTA for 50 hr integration time555https://www.cta-observatory.org.

The $\gamma$-ray spectrum of undetected sources remains hard up to $\sim$100 GeV ($\frac{{\rm d}N}{{\rm d}E}\propto E^{-2}$) and declines after, most likely due to EBL absorption. The average spectrum of the Fermi-LAT detected sources also reveals a well defined peak at $\sim$100 GeV and steepens at higher energies. To explore the role of the EBL, we generate stacked spectra of both populations taking into account the EBL absorption (Domínguez et al., 2011). The results are plotted in Figure 2 (left panel) with cyan and gray colors for $\gamma$-ray undetected and detected extreme blazars, respectively. Since the EBL corrected stacked spectra do not exhibit any softening, the inverse Compton peak of the Fermi-LAT detected population lies above 1 TeV, thus indicating a very hard intrinsic TeV spectrum666We have also tested whether the shape of the stacked SED is dominated by a few bright sources. This is done by normalizing the SED of each blazar by its integrated 10$-$1000 GeV flux. We find that such a normalized stacked spectrum has a shape similar to that plotted in Figure 2 and therefore the results shown here truly represent the average behavior of the population.. A similar result holds for the $\gamma$-ray undetected sources though a strong claim cannot be made due to the flux upper limit in the highest energy bin.

The derived results also demonstrate the effectiveness of the developed stacking technique to extract the $\gamma$-ray signal from a population which is $\sim$an order of magnitude fainter than the Fermi-LAT detected one. By comparing this signal with the plotted sensitivity limits, we find that the probability of CTA detection is higher for objects that already have Fermi-LAT detections (Figure 2). The $\gamma$-ray undetected extreme blazar population, on the other hand, remains below the detection limit of CTA both due to their low level of $\gamma$-ray emission and the EBL absorption (see also Franceschini et al., 2019). This result does not rule out possibility of CTA detection for a few individual sources, especially during a TeV flaring state. However, considering the source population as a whole, we show that Fermi-LAT detected extreme blazars are better candidates for VHE observations.

With the knowledge of the average $\gamma$-ray spectrum, we can derive the contribution of the Fermi-LAT undetected extreme blazars to the total EGB. It is computed by assuming that the sources are uniformly distributed in the sky outside the Galactic plane777All extreme blazars studied in this work lie outside the Galactic plane. The assumption of the uniform sky distribution allows us to derive the lower limit to the EGB contribution by the astrophysical population under consideration. The accurate measurement can be done by considering the source count distribution or luminosity function (see, e.g., Ackermann et al., 2016) and requires a precise estimation of the solid angle covered in the sky. These aspects, however, are beyond the scope of this work. ($|{\rm b}|>10^{\circ}$). In the right panel of Figure 2, we show the derived results and compare them with the total EGB estimated in Ackermann et al. (2015). According to our analysis, $\gamma$-ray undetected extreme blazars are responsible for at least $\sim$10% of the total EGB at 100 GeV.

4 Validation of the Stacking Technique

4.1 Background Stacking

Since we combine the $\gamma$-ray undetected objects, a genuine question arises about the possibility of stacking the diffuse background emission rather than the radiation originated from actual point sources. Therefore, to test the robustness of the $\gamma$-ray signal reported above, we randomly select 172 empty sky positions not lying within 95% error radius of any 4FGL source and repeat the same procedure as adopted for extreme blazars. In the left panel of Figure 1, we show the TS histogram of these random positions. Comparing the distribution with the null hypothesis ($\chi^{2}$ curve), it can be concluded that the derived signal is fully compatible with the random background fluctuations. The corresponding stacking plot shown in the left panel of Figure 3 confirms the negligible contribution of the random background fluctuations to the stacked $\gamma$-ray emission observed from extreme blazars.

4.2 Stacking of Simulated Sources

The robustness of the stacking technique is also verified by performing simulations of $\sim$127 months of the Fermi-LAT data for 100 objects. The $\gamma$-ray spectra of these sources are assumed to follow a power law. The power law indices are randomly extracted from a Gaussian distribution peaking at 2.1 with a dispersion of 0.1. On the other hand, fluxes are selected from a log-normal distribution having the peak at 10*-11* ${\rm\,ph\,\,cm^{-2}\,s^{-1}}$ and a dispersion of 0.1 in log-space. We populate the background sky with 4FGL sources and Galactic and isotropic diffuse emissions and repeat the entire procedure as described above. The results of the stacking of the simulated objects are presented in the right panel of Figure 3. We find the best-fit photon flux and index values as $1.04^{+0.01}_{-0.01}\times 10^{-11}$ ${\rm\,ph\,\,cm^{-2}\,s^{-1}}$ and 2.11${}^{+0.07}_{-0.07}$, respectively at peak TS = 1216. Clearly, the agreement between the input and output spectral parameters confirms the effectiveness of the developed tool in measuring the average $\gamma$-ray properties of the unresolved population.

5 Discussion and Summary

We have described a stacking technique to extract the $\gamma$-ray signal from the Fermi-LAT undetected sources by combining their likelihood profiles. The pipeline is applied to a sample of 172 extreme blazars resulting in a strong $\gamma$-ray detection of the population. The tool is capable of extracting about an order of magnitude fainter $\gamma$-ray signal than that possible with the conventional point source Fermi-LAT data analysis. A crucial finding of the stacking analysis is that CTA may not be able to detect VHE emission from these objects mainly due to their faintness and EBL absorption. However, sources that already have been detected by Fermi-LAT should be primary targets for CTA observations. Another important finding is that the EBL-corrected stacked spectrum of extreme blazars exhibits no softening up to 1 TeV, thus indicating their inverse Compton peak to lie at $>$1 TeV. At 100 GeV, a significant fraction ($\sim$10%) of the total EGB originate from the $\gamma$-ray undetected extreme blazars.

The simplicity of the developed technique offers a possibility to apply it to various astrophysical problems other than just to determine the $\gamma$-ray detection/non-detection of a source population. One such example could be the known correlation between the IR and $\gamma$-ray luminosities ($L_{\rm IR}$ and $L_{\gamma}$, respectively) in star-forming galaxies (e.g., Ackermann et al., 2012). Since $L_{\rm IR}$ is well-known, the stacking technique can be used to explore the strength of the correlation and associated parameters that maximizes the TS profile for a given set of $\gamma$-ray spectral parameters (and thus $L_{\gamma}$). This is fully explored in a companion paper where we apply the developed pipeline to a sample of star-forming galaxies (Ajello et al., in preparation). Finally, the stacking technique is flexible to adopt any spectral models (other than power law used here), including EBL attenuated ones, and can be applied to any astrophysical populations, e.g., galaxy clusters and X-ray binaries (cf. Dubus, 2013). This makes it a versatile tool for $\gamma$-ray astronomy.

We are thankful to the referee for a constructive criticism. The Fermi-LAT Collaboration acknowledges support for LAT development, operation and data analysis from NASA and DOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN (Italy), MEXT, KEK, and JAXA (Japan), and the KA Wallenberg Foundation, the Swedish Research Council and the National Space Board (Sweden). Science analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged. This work performed in part under DOE Contract DE-AC02-76SF00515. .