Identifying TeV Source Candidates among Fermi-LAT Unclassified Blazars

Blazars, some of the most powerful active galactic nuclei (AGNs), have a relativistic jet pointing toward the observer (e.g., Abdo et al. 2010; Massaro et al. 2015), and show rapid variability and high optical and radio polarization. Such objects are the most numerous class of extragalactic sources detected by TeV telescopes, the most sensitive of which are the Imaging Atmospheric Cerenkov Telescopes (IACTs) such as MAGIC,⁵ H.E.S.S.,⁶ VERITAS,⁷ and the upcoming Cerenkov Telescope Array (CTA).⁸ Despite their high sensitivity, however, observations by current IACTs are limited by their fairly small fields of view (∼few degrees), weather conditions, the need for relatively dark night skies, and by a high background that requires fairly long observations. A source with a flux of ∼1% of the Crab Nebula flux requires around 50 hr of observation time for a detection at 5σ. IACTs typically take data for only about 1200 hr per year (De Naurois & Rolland 2009). Such constraints provide a strong incentive to identify likely targets for IACT observations.

Blazar spectral energy distributions (SEDs) show two broad peaks in a νf_ν representation. The low-energy peak is attributed to synchrotron radiation, and the high-energy peak is usually thought to be produced by inverse Compton radiation (IC; e.g., Sikora et al. 1994). Based on the position of the synchrotron peak (${\nu }_{\mathrm{peak}}^{S}$) in the SED, blazars are divided into three subclasses: low-synchrotron-peaked (LSP, with ${\nu }_{\mathrm{peak}}^{S}\leqslant {10}^{14}\,$ Hz), intermediate-synchrotron-peaked (ISP, with ${10}^{14}\,\mathrm{Hz}\,\lt {\nu }_{\mathrm{peak}}^{S}\leqslant {10}^{15}\,$ Hz) and high-synchrotron-peaked (HSP, with ${\nu }_{\mathrm{peak}}^{S}\gt {10}^{15}\,$ Hz) (Abdo et al. 2010). HSPs, primarily BL Lac objects, represent the most numerous class of extragalactic TeV-energy sources. The TeVCat⁹ is an online, interactive catalog for very-high-energy (VHE energies, E > 100 GeV) γ-ray astronomy (Horan & Wakeley 2008). This catalog reports 223 TeV sources as of this writing. Among the 61 objects associated with blazars, 51 of them are HSPs and only 10 are LSP/ISP flat-spectrum radio quasars (FSRQs).¹⁰

All-sky observations with the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi; Atwood et al. 2009) at GeV energies offer opportunities to find such targets. An example is the Third Catalog of Hard Fermi-LAT Sources (3FHL: Ajello et al. 2017), which reports the locations and spectra of sources significantly detected in the 10 GeV–2 TeV energy range during the first 7 yr of the Fermi mission. From the 3FHL it is possible to select TeV candidates by γ-ray flux and photon index.

An alternative approach to searching for TeV candidates is to find objects belonging to a class of sources likely to be seen at TeV energies. In the case of blazars, this can be done by identifying those objects whose synchrotron emission peaks at high frequencies. An example of this approach is the second Wide-field Infrared Survey Explorer (WISE) High Synchrotron Peak Catalog (2WHSP; Chang et al. 2017), which is a list of HSP candidates based on multifrequency analysis of γ-ray source candidates with $| b| \gt 10^\circ$.

Here, we present a third approach to search for TeV HSP source candidates. This method includes two steps: (1) we use γ-ray variability information to search out potential HSPs among the unclassified Fermi-LAT sources and (2) we analyze γ-ray spectra of these sources using more Fermi-LAT data than used in published catalogs. The starting point is the third Fermi-LAT all-sky catalog of sources detected at energies between 100 MeV and 300 GeV (3FGL: Acero et al. 2015). The 3FGL catalog lists 3033 γ-ray sources, of which 1745 are AGNs, mostly BL Lacs and FSRQs, and includes γ-ray source locations, energy spectra, variability information on monthly timescales, and likely associations with objects seen at other wavelengths. In this catalog, 573 sources are listed as blazars of uncertain type (BCUs) and 1010 objects lack a plausible counterpart at other wavelengths (Unassociated Catalog Sources; UCSs).¹¹

Although BCUs and UCSs often lack optical spectra and sufficient information for a rigorous classification, statistical methods such as the Artificial Neural Network (ANN) algorithm can potentially provide classifications of these sources (e.g., Lefaucheur & Pita 2007; Chiaro et al. 2016; Saz Parkinson et al. 2016; Salvetti et al. 2017). In particular, Saz Parkinson et al. (2016) found 559 of the UCS sources have characteristics similar to those of AGN. These UCS_agns are combined with the original 573 BCUs from the 3FGL catalog to provide the targets for our search for HSP/TeV candidates.

The paper is organized as follows. In Section 2 we present the machine-learning method used in this study; in Section 3 we describe the selection of HSP candidates among the uncertain 3FGL objects; and in Section 4 we discuss the results of a dedicated Fermi-LAT analysis of the sources found analyzing 104 months of data. In Section 5 we examine the detectability of the targets by the present generation of IACTs and CTA. In Section 6 we summarize the conclusions of this study.

The starting point for selecting HSP candidates is the ANN method previously applied to Fermi-LAT sources to distinguish FSRQ-like sources from those with BL Lac characteristics (Chiaro et al. 2016; Salvetti et al. 2017). The key idea is that the γ-ray flares of BL Lacs tend to be smaller and less frequent than those of FSRQs. The input is the empirical cumulative distribution function (ECDF) of monthly γ-ray flux values for each source, taken from the 3FGL. We included in the algorithm values corresponding to increments of 10% from the 10th to the 100th percentile of the ECDF. The algorithm computes a likelihood value arranged to have two possibilities: class A or class B, with a likelihood (L) assigned to each analyzed source so that the likelihoods to belong to one or the other of the two classes are related by L_A = 1 − L_B. In this way, the greater the value of L_A, the greater the likelihood that the source is a class A candidate. In this case, the likelihood applies to the source having HSP characteristics, L_HSP. This approach uses the two-layer-perceptron ANN technique (Gish 1990; Bishop 1995), which is probably the most widely used architecture for practical applications of neural networks. Data enter the neural network through nodes in the input layer. The information travels across the links and is processed in the nodes through an activation function. Each node in the output layer returns the likelihood of a source to be a specific class. We applied the algorithm to the three synchrotron peak subclasses as classified in the Third Catalog of AGNs detected by the Fermi-LAT (3LAC: Ackermann et al. 2015) to train it to distinguish each source class.

Repeating the analysis of Chiaro et al. (2016) as a cross-check on that analysis, we considered 289 HSPs and 824 non-HSP objects classified in the third Fermi-LAT AGN catalog (3LAC; Ackermann et al. 2015). Maintaining the same ratio as in the catalog, that is, one third HSPs and two thirds non-HSP sources, we randomly mixed the sample and divided it into three subsamples: training, validation, and testing. The training sample is used to optimize the network. The validation sample is used to avoid over-fitting. The testing sample is independent both of the training and validation ones and was used to monitor the accuracy of the network. Although the random sampling resulted in a different training set from the previous analysis, the results were the same: for L_HSP > 0.8, 75% of the sources have characteristics of HSPs, while for L_HSP > 0.89, we expect 90% of the sources to be HSP-like.

We then applied the ANN HSP algorithm to the 573 BCUs and the 559 UCS_agn of our sample. The resulting likelihood distributions plotted in Figure 1 show, as expected, a peak at L_HSP = 0.0 due to the non-HSP populations (ISP and LSP), which are much more numerous than HSPs at energies covered by the Fermi-LAT. As for the original analysis (Chiaro et al. 2016), the lack of a peak at L_HSP = 1.0 indicates that the ANN network was not able to separate HSPs cleanly, but for the purpose of selecting candidates for additional analysis we are primarily interested in finding a high fraction of the sources with the desired characteristics. Requiring the L_HSP > 0.8 value, we identified 48 BCU and 32 UCS_agn as HSP candidates. Applying the higher threshold value L_HSP > 0.89, we identified 11 BCUs and 5 UCS_agn as very high confidence (VHC) HSP candidates. Tables 1 and 2 show the full lists of candidates. The names of the VHC candidates are shown in bold.

Zoom In Zoom Out Reset image size

We compared the results of our analysis with the list of blazars presented in the 2WHSP catalog (Chang et al. 2017), as both methods attempt to identify HSP objects. Of the 11 VHC candidates from the BCU list, six were also identified by 2WHSP: 3FGL J0506.9−5435, 3FGL J0921.0−2258, 3FGL J1155.4−3417, 3FGL J1711.6+8846, 3FGL J1714.1−2029, and 3FGL J1944.1−4523. None of the five VHC candidates from the UCS list were found in the 2WHSP catalog. Similar fractions were found for the other parts of the lists, indicating that the two methods are complementary in finding candidate HSPs.

Because the 4FGL catalog will soon be available with a larger sample of sources to study, we chose to focus our spectral analysis for this new method of selecting HSP candidates on the 16 VHC sources. We analyzed 104 months of Fermi-LAT Pass 8 (Atwood et al. 2013) data from 2008 August 4 to 2017 April 4, selecting γ-ray events in the energy range E = [0.1, 1000] GeV, passing standard data quality selection criteria and zenith angle cuts for AGN (e.g., Meyer et al. 2019), to find the γ-ray properties of our HSP candidates. We considered events belonging to the Pass 8 SOURCE event class and used the corresponding instrument response functions P8R2_SOURCE_V6, as we were interested in point-source detection. We used the interstellar emission model (IEM) released with Pass 8 data (Acero et al. 2016) (i.e., gll_iem_v06.fits). This is the model recommended for use with Pass 8 analyses. We also included the standard template for the isotropic emission (iso_P8R2_SOURCE_V6_v06.txt).¹²

We developed an analysis pipeline using FermiPy, a Python package that automates analyses with the Fermi Science Tools (Wood et al. 2017).¹³ FermiPy includes tools that can (1) generate simulations of the γ-ray sky, (2) detect point sources, and (3) calculate the characteristics of their SEDs.

The likelihood analysis works on a square region of interest (ROI). We used a 16° × 16° ROI centered on the sources of our sample. We analyzed each ROI separately. In each ROI, we binned the data with a pixel size of 008 and 8 energy bins per decade. Our background model included the IEM, isotropic template and sources from the preliminary 8 yr list, FL8Y,¹⁴ except for the source being analyzed. We allowed the normalization and slope of the IEM and the normalization of the isotropic template to vary. We first relocalized the source of interest, and then we searched for new point sources with Test Statistic TS > 25, defined as twice the difference between the log-likelihood for the null hypothesis (no source) and the hypothesis of a source at the location. New sources were then added to the analysis. We then analyzed each of our HSP candidates, assuming a power-law spectrum after determining that all the VHC candidates used that spectral form in the 3FGL catalog (Acero et al. 2015). For the individual energy bins, we plotted as upper limits those points with a $\mathrm{TS}\lt 9$.

In Tables 1 and 2, we report the following parameters for the HSP candidates: the best fit and 1σ error of the position, the TS of detection, the photon index found for a power-law SED shape, the flux, and the value of L_HSP. The photon index parameter, if less than 2, can be a relevant indicator for an IC peak at TeV energies and therefore quite useful in selecting IACT and CTA candidates. The mean and rms of the photon indexes of HSPs are 1.87 ± 0.20 while for LSPs and ISPs these are 2.21 ± 0.18, 2.07 ± 0.20 respectively (Ackermann et al. 2015).

We have repeated the analysis with a log-parabola intrinsic spectrum instead of a power law. Only for three sources (3FGL J0921.0–2258, 3FGL J0153.4+7114, and 3FGL J0506.9–5435) did we find that a log-parabola is preferred at the 4σ level, i.e., TS_curv > 16, where TS_curv is twice the difference of the log-likelihood values of the best fit with a log-parabola and a power law, respectively. These results are largely consistent with the spectral parameters listed in the preliminary 4FGL catalog.¹⁵ The 4FGL sources corresponding to 3FGL J0921.0−2258 and 3FGL J0153.4+7114 have power-law spectra instead of the curved spectra we found, although the evidence for curvature is of low significance in our analysis. If we had used the power-law fit, the extrapolations would have obviously been higher.

In this section, we compare the extrapolated fluxes of the Fermi-LAT spectra of our VHC sources with the sensitivity of present IACTs and the future CTA. In the case that the log-parabola fit to the Fermi-LAT data is preferred with TS_curv > 16, we use the log-parabola parameterization for the extrapolation and the best-fit power law otherwise. To evaluate whether the VHC HSP candidates can realistically be observed with IACTs or CTA, we must take into account the interaction of γ-rays with photons of the extragalactic background light (EBL). The EBL spans the wavelength regime from ultraviolet to far-infrared wavelengths and mainly consists of the integrated starlight emitted over the history of the universe and starlight absorbed and reemitted by dust in galaxies (Hauser & Dwek 2001; Kashlinsky 2005). During the propagation of γ rays through the EBL, the electron–positron pairs produced via γγ → e⁺ e⁻leads to an attenuation of the initial γ-ray flux (Nikisov 1962; Gould & Schreder 1967; Dwek & Krennrich 2013). To properly evaluate the absorption effect of the EBL, it is necessary to know the redshift of the analyzed γ-ray source. Because the redshifts of the selected VHC HSP candidates are unknown, we assume redshifts of z = 0 and z = 0.5 for the calculation. These are typical values of observed γ-ray BL Lacs.¹⁶ We use these z values while acknowledging that the blazar redshift range is a very open and longstanding debate. Some authors argue that the BL Lacs without a redshift are likely much more distant than those with a measured one (e.g., Padovani et al. 2012), so it could be possible that some objects fall beyond the 0.5 redshift.

Using the EBL model of Dominguez et al. (2011), we extrapolate the best-fit spectra obtained with the Fermi-LAT analysis in Section 4 up to 10 TeV with the assumed redshift values. The results are shown in Figures 2 and 3 where we compare the extrapolated fluxes with the CTA sensitivity (Acharyya et al. 2019) for 50 hr (5 hr) of observations as a solid (dashed) gray line. The sensitivities for currently operating IACT arrays H.E.S.S., MAGIC, and VERITAS are shown as a red band (Holler et al. 2015; Aleksć et al. 2016).¹⁷ The CTA sensitivity curves are available for the northern and southern array and for zenith angles (Z) of 20° and 40°. We choose the sensitivity curve depending on the source decl. δ assuming CTA site locations at 287 northern latitude and 247 southern latitude. For δ > 587, we choose the northern array with Z = 40°, for 587 ≤ δ < 27 the northern array with Z = 20°, for 2.7 ≥ Z ≥ − 547 the southern array with Z = 20°, and the southern array with Z = 40° for δ < −547.¹⁸

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Table 3, we report an estimate of the maximum redshift values of the VHC TeV candidates so that the sources are still detectable at 5σ for 50 hr observations with current IACTs as well as 5 and 50 hr CTA observations. For these estimates, we assume that the spectra are detectable if the extrapolation is above the sensitivity curves. We also report how these values change if the best-fit spectral parameters are changed by their 1σ statistical uncertainties. An actual detection will depend on the spectral shape and the exact exposure time, therefore these numbers should be regarded as a rough estimate. If no value is given, the source might not be significantly detected with the assumed observation within the assumed redshift range.

Table 3.  Estimates of the Maximum Redshift Values of the VHC Sources for Which the Sources Are Still Detectable at 5σ in 50 hr of Current IACT and 50 (5)hr of CTA Observations, Respectively

Source Name	Tobs = 50 hr, Current IACTs	Tobs = 5 hr, CTA	Tobs = 50 hr, CTA
BCU HCTeV candidates
3FGL J0047.9+5447	${0.18}_{-}^{\gt 0.50}$	${0.15}_{-}^{\gt 0.50}$	>0.50−0.38
3FGL J1155.4−3417	${0.26}_{-0.17}^{\gt 0.50}$	${0.35}_{-0.24}^{\gt 0.50}$	>0.50
3FGL J1434.6+6640	${0.15}_{-0.13}^{+0.27}$	${0.14}_{-}^{+0.33}$	>0.50−0.31
3FGL J0921.0−2258	⋯	⋯	⋯
3FGL J0648.1+1606	${0.03}_{-}^{+0.19}$	${0.01}_{-}^{+0.19}$	${0.29}_{-0.26}^{\gt 0.50}$
3FGL J1711.6+8846	${0.04}_{-}^{+0.22}$	${0.02}_{-}^{+0.24}$	${0.32}_{-0.29}^{\gt 0.50}$
3FGL J1714.1−2029	${0.38}_{-0.28}^{\gt 0.50}$	${0.47}_{-0.35}^{\gt 0.50}$	>0.50
3FGL J1910.8+2855	${0.18}_{-0.15}^{+0.25}$	${0.16}_{-0.15}^{+0.28}$	>0.50−0.17
3FGL J0153.4+7114	⋯	⋯	${0.43}_{-}^{\gt 0.50}$
3FGL J0506.9−5435	${0.05}_{-}^{\gt 0.50}$	${0.25}_{-}^{\gt 0.50}$	>0.50—
3FGL J1944.1−4523	${0.27}_{-0.21}^{\gt 0.50}$	${0.34}_{-0.26}^{\gt 0.50}$	>0.50−0.07
UCSagn HCTeV candidates
3FGL J1549.9−3044	${0.13}_{-}^{+0.27}$	${0.15}_{-0.13}^{\gt 0.50}$	>0.50−0.29
3FGL J2142.6−2029	${0.05}_{-}^{+0.28}$	${0.07}_{-}^{+0.35}$	${0.41}_{-0.35}^{\gt 0.50}$
3FGL J2321.6−1619	${0.23}_{-0.15}^{+0.20}$	${0.30}_{-0.20}^{\gt 0.50}$	>0.50
3FGL J2145.5+1007	${0.04}_{-}^{+0.33}$	${0.02}_{-}^{+0.35}$	${0.26}_{-}^{\gt 0.50}$
3FGL J2300.0+4053	${0.15}_{-}^{+0.32}$	${0.13}_{-}^{+0.35}$	>0.50−0.37

Note. If no value is given, the source will not be significantly detected within the assumed observation time. Sub- and superscript numbers give the change in the redshift values if the Fermi-LAT spectrum is extrapolated within 1σ uncertainties. See the main text for further details and caveats.

Download table as:  ASCIITypeset image

Any of the VHC sources that lie at low redshift should be detectable with CTA, with the exception of 3FGL J0921.0−2258, which shows strong intrinsic curvature. Similarly, most sources should already be detectable by the present-generation IACTs if their redshifts are close to 0. For a redshift of z = 0.5, most sources appear to be good candidates for detection with 50 hr of CTA observations, primarily in the 100 GeV energy range. Only 3FGL J1714.1−2029 seems possible as a candidate for current IACTs if its distance is at the upper end of the range considered, due to the EBL attenuation.

Motivated to expand the sample of TeV sources, we applied a machine-learning algorithm to variability parameters of Fermi-LAT blazar-like sources without firm identifications combined with new analysis of the LAT data for these sources. Follow-up work will require additional multiwavelength studies, including finding redshifts for most of the candidates and targeted observation by IACTs. We also recognize that this search is necessarily incomplete because of the difficulty to distinguish the blazar subclasses by the γ-ray flux properties only. As already pointed out by Ackermann et al. (2015), the γ-ray sources with unknown properties are generally fainter than the well-defined classes. The fainter sources offer less of the flaring information needed for the machine-learning method, so there may be HSP blazars among the sample of 3FGL sources rejected in the first step of our method. The level of incompleteness is difficult to quantify due to the very similar values of the synchrotron peaks of the three blazar subclasses. Nevertheless, some VHC HSP candidates, also thanks to the analysis of Fermi-LAT data, are convincing as TeV candidates and should be promising targets for currently operating IACTs, especially if the sources are located below redshifts of ∼0.1. If the sources are further away, CTA should be capable of significantly detecting them.

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT, as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

M.D.M. acknowledges support by the NASA Fermi Guest Investigator Program 2014 through the Fermi multi-year Large Program No. 81303 (PI: E. Charles) and by the NASA Fermi Guest Investigator Program 2016 through the Fermi one-year Program No. 91245 (PI: M. Di Mauro).

This paper has gone through internal review by the CTA Consortium. This research has made use of the CTA instrument response functions provided by the CTA Consortium and Observatory, see http://www.cta-observatory.org/science/cta-performance/ (version prod3b-v1) for more details.

Facilities: Fermi - Fermi Gamma-Ray Space Telescope (formerly GLAST), H.E.S.S. - , MAGIC - , VERITAS. -