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

The Large Area Telescope (LAT) on board the Fermi Gamma ray Space Telescope has revealed various types of astrophysical objects as γ-ray emitters. The most numerous of them are blazars, i.e., radio-loud quasars with powerful relativistic jets pointed toward the observer, followed by narrow-line Seyfert 1 galaxies, pulsars, and many others (see 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 γ-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 γ-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 γ-ray band. They are, similar to the Fermi-LAT detected objects, crucial to study their contribution to the diffuse γ-ray background emission. Focusing on blazars, a characterization of the γ-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; McEnery et al. 2019) and the Cherenkov Telescope Array (CTA; Actis et al. 2011). In particular, AMEGO (with an energy range 200 keV to ≳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 ∼0.02–300 TeV band, will observe the most efficient particle accelerator jets from BL Lac objects, a subclass 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 was successfully applied earlier to Energetic Gamma-Ray Experiment Telescope data to search for the γ-ray signal from, e.g., clusters 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 γ-ray emission from undetected populations. Huber et al. (2012) proposed a method which coadds the Fermi-LAT count maps and performs a maximum likelihood analysis on the combined data to derive the γ-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 γ-ray properties of the Fermi-LAT undetected population. The developed technique is sensitive to extract the faint γ-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 γ-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σ level, unless specified.

The input for the stacking pipeline is the list of sky positions, i.e., R.A. and decl., for the objects under consideration. The tool works in the following two steps.

2.1. Preprocessing

2.2. Stacking

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 subsample of BL Lac objects that have the synchrotron peak⁵ located at very high frequencies (${\nu }_{\mathrm{Syn}}^{\mathrm{peak}}\,\geqslant {10}^{17}\,\mathrm{Hz}$) indicating they 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 (VHEs; >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 γ-ray undetected extreme blazars and determine whether they are GeV γ-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 }_{\mathrm{Syn}}^{\mathrm{peak}}\,\geqslant {10}^{15}\,\mathrm{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).⁶ 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 ∼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 point spread function (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 (see Acero et al. 2016). We use a squared ROI of 10° × 10° and adopt a zenith angle cut (z_max) of 105°. With these settings, the preprocessing led to a significant γ-ray detection of 165 extreme blazars.⁷ The TS distribution of the remaining 172 γ-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., χ² distribution with 2 degrees of freedom, at low TS, and shows excess above TS ≳ 10 that can be attributed to the jet emission.

Zoom In Zoom Out Reset image size

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 (∼32σ significance for 2 degrees of freedom). The associated peak γ-ray flux and photon index are ${F}_{10-1000\mathrm{GeV}}\,={6.51}_{-0.35}^{+0.36}\times {10}^{-12}$ ph cm⁻² s⁻¹ and ${{\rm{\Gamma }}}_{10-1000\mathrm{GeV}}={2.08}_{-0.06}^{+0.07}$. These results demonstrate the capabilities of the stacking technique to extract the faint γ-ray signal from a Fermi-LAT undetected population.

The strong γ-ray signal has also allowed us to perform the stacking in smaller energy bins and generate a stacked γ-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 γ-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 time.⁸

Zoom In Zoom Out Reset image size

The γ-ray spectrum of undetected sources remains hard up to ∼100 GeV $\left(\tfrac{{dN}}{{dE}}\propto {E}^{-2}\right)$ 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 ∼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 γ-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 spectrum.⁹ A similar result holds for the γ-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 γ-ray signal from a population that is ∼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 γ-ray undetected extreme blazar population, on the other hand, remains below the detection limit of CTA both due to their low level of γ-ray emission and the EBL absorption (see also Franceschini et al. 2019). This result does not rule out the 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 γ-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 plane¹⁰ ($| {\text{}}b| \gt 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, γ-ray undetected extreme blazars are responsible for at least ∼10% of the total EGB at 100 GeV.

4.1. Background Stacking

Since we combine the γ-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 γ-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 (χ² 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 γ-ray emission observed from extreme blazars.

Zoom In Zoom Out Reset image size

4.2. Stacking of Simulated Sources

We have described a stacking technique to extract the γ-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 γ-ray detection of the population. The tool is capable of extracting about an order of magnitude fainter γ-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 (∼10%) of the total EGB originates from the γ-ray undetected extreme blazars.

The simplicity of the developed technique offers a possibility of applying it to various astrophysical problems other than just to determine the γ-ray detection/nondetection of a source population. One such example could be the known correlation between the IR and γ-ray luminosities (L_IR and L_γ, respectively) in star-forming galaxies (e.g., Ackermann et al. 2012). Since L_IR is well known, the stacking technique can be used to explore the strength of the correlation and associated parameters that maximize the TS profile for a given set of γ-ray spectral parameters (and thus L_γ). This is fully explored in a companion paper where we apply the developed pipeline to a sample of star-forming galaxies (M. Ajello et al. 2019, 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 (see Dubus 2013). This makes it a versatile tool for γ-ray astronomy.

We are thankful to the referee for 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. M.A. and V.S.P. acknowledge funding under NASA contract 80NSSC18K1718.

Software: fermiPy (Wood et al. 2017).