Searching for TeV Candidates in 4LAC High-synchrotron-peaked Frequency BL Lac Objects

Most extragalactic sources detected in the $\gamma$-ray band belong to the blazar category(Abdollahi et al., 2020). Blazars are an important subclass of active galactic nuclei (AGNs) and are characterized by their strong and rapid variability and high levels of brightness (e.g., Blandford & Rees 1978; Urry & Padovani 1995). The spectral energy distributions (SEDs) of blazars are dominated by two components, which are illustrated by a double-bump spectral shape in $\rm log\nu$-$\rm log\nu F\nu$ space. The origin of the low energy bump, seen from the radio band to the ultraviolet or soft X-ray band, is attributed to the synchrotron emission of a relativistic electrons in the jet. Either leptonic models (e.g., Dermer et al. 1992; Maraschi et al. 1992; Dermer & Schlickeiser 1993; Bloom & Marscher 1996; Zheng & Yang 2016) or hadronic models (e.g., Aharonian 2000; Mücke & Protheroe 2001; Mücke et al. 2003; Zheng & Kang 2013) can be used to reproduce the high energy emission of blazars. According to the presence or absence of broad emission lines in their optical spectra, blazars are divided into flat spectrum radio quasars (FSRQs) and BL Lac objects. The equivalent widths (EWs) of FSRQ optical spectra emission lines in the comoving frame are greater than 5$\rm\AA$, while the EWs of BL Lac objects are less than 5$\rm\AA$ (Stickel et al., 1991). The peak frequency, $\nu_{\rm syn}$, of the low energy bump (synchrotron bump) can also be used to classify blazars, as follows: low-synchrotron-peak blazars (LSP; $10^{14}\leq\nu_{\rm syn}\leq 10^{15}\rm{Hz}$); intermediate-synchrotron-peak blazars (ISP; $10^{14}\leq\nu_{\rm syn}\leq 10^{15}\rm{Hz}$); high-synchrotron-peak blazars (HSP; $\nu_{\rm syn}\textgreater 10^{15}\rm{Hz}$) (Abdo et al., 2010); and extreme HSP blazars (EHSP, $\nu_{\rm syn}\textgreater 10^{17}\rm{Hz}$) (Arsioli et al., 2018). Multiwavelength observations show that the origin of the nonthermal emission of blazars extends from the radio band to the $\gamma$-ray band. Some can even extend to the very high energy band (VHE, $E\textgreater 0.1\ \rm{TeV}$, e.g., de Oña Wilhelmi 2008).

VHE $\gamma$-ray astronomy research focuses on the high energy emissions from extreme objects and their physical mechanisms. Research regarding TeV sources has promoted the development of both VHE $\gamma$-ray astronomy and neutrino astronomy (Arsioli et al., 2015; Chang et al., 2017, 2019; Chiaro et al., 2019). However, this research has long been hindered by limited observations across the TeV energy band and the limited number of TeV sources. Currently, the observation of TeV sources mainly depends on imaging atmospheric Cerenkov telescopes (IACTs), such as WHIPPLE, MAGIC, H.E.S.S., and VERITAS, or extensive air shower (EAS) arrays, such as Tibet AS-$\gamma$ and ARGO-YBJ. Traditional ground-based detectors suffer from a small field of view (FOV), limited sensitivity, and short exposure times. They are also affected by weather and the Earth’s magnetic field. In addition, although a number of sources emit TeV $\gamma$-rays, it is difficult for TeV photons to reach the Earth because of the absorption caused by photons from extragalactic background light (EBL) or cosmic microwave background (CMB), with which they interact. This is especially the case for distant sources (Hauser & Dwek, 2001; Dom´ınguez et al., 2011; Franceschini & Rodighiero, 2018). Taking these various issues together, TeV sources are difficult to detect and identify. The online VHE $\gamma$-Ray Source Catalog, TeVCat¹¹1TeVCat http://tevcat.uchicago.edu/ (Horan & Wakely, 2008) reports 232 TeV sources so far. There are 77 TeV blazars, including 66 BL Lacs, which dominate the extragalactic TeV sources. The next generation of TeV detectors, i.e., the Cherenkov Telescope Array (CTA) and the Large High Altitude Air Shower Observatory (LHAASO), are expected to perform significantly better. This makes it worthwhile to search for a collection of high-confidence TeV candidates for follow-up observations.

The Fermi Large Area Telescope (Fermi-LAT) has observed the entire sky at GeV energy levels for more than 12 years (Atwood et al., 2009). Making use of these Fermi GeV $\gamma$-ray observations provides an alternative approach to explore for TeV candidates. On the basis of the first seven years of observations, the Third Catalog of Hard Fermi-LAT Sources (3FHL, Ajello et al. 2017) reported 1556 sources detected in the energy range from 10 GeV to 2 TeV. From these sources, TeV candidates can be selected based on their GeV fluxes and $\gamma$-ray spectral indices. When focusing on TeV blazars, HSP blazars, especially HSP BL Lac objects (HBLs), can serve as a good data set from which to select TeV candidates. Several HSP blazar catalogs, such as the second Wide-field Infrared Survey Explorer (WISE) High Synchrotron Peak Catalog (2WHSP, Chang et al. 2017) and the catalog of extreme and high-synchrotron peak blazars (3HSP; Chang et al. 2019) have been utilized in the search for blazars that are likely to be detected by their TeV energies. Using artificial neural network (ANN) machine-learning (ML) algorithms, Chiaro et al. (2019) searched for HBL candidates in the collections of unclassified blazars and unassociated sources observed with Fermi-LAT, relying on the GeV data obtained from the 3FGL. They then analyzed the Fermi spectra of the HBL candidates and used an EBL model to identify TeV candidates. In other bands outside of the GeV $\gamma$-ray band, TeV blazars also show unique characteristics. Massaro et al. (2013) have suggested that TeV BL Lac objects have a higher X-ray to infrared flux ratio and a lower WISE magnitudes and that these characteristics could be used to find TeV candidates. In an alternative approach, Costamante (2020) combined X-ray, $\gamma$-ray, and infrared data to analyze the clustering of TeV BL Lac objects in the entire BL Lac objects and managed to identify some TeV candidates.

Recently, the Fermi-LAT collaboration released the fourth Fermi-LAT Gamma-ray source catalog (4FGL-DR1, eight years of exposure; Abdollahi et al. 2020), a fourth catalog of AGNs detected by the Fermi-LAT (4LAC, Ajello et al. 2020), and a second release of 4FGL data (4FGL-DR2, 10 years of exposure; Ballet et al. 2020). 4FGL contains more gamma-ray sources than the previous source catalog, has higher energy upper-limits, and has longer exposure times (Acero et al., 2015; Abdollahi et al., 2020). In addition, the reports provide wide cross-matching with radio observations (VLBI Radio Fundamental Catalog, e.g., RFC; see Petrov et al. 2019) and optical observations (Gaia Data Release 2, e.g., Gaia-DR2, see Gaia Collaboration 2018) . These advances facilitate the exploration for candidates for TeV observations. In this paper, we present a new approach to search for TeV candidates in the 4LAC HBLs, which incorporates two steps: 1) radio, optical, and GeV ${\gamma}$-ray data are combined and several ML algorithms are employed to search for possible TeV candidates (PTCs) in the 4LAC HBLs; 2) the ${\gamma}$-ray spectra of the PTCs are analyzed by using the longer exposure times and higher energy upperlimits of the new Fermi data than reported in the published catalogs. The EBL model, FOV, and sensitivity of CTA and LHAASO are then used to evaluate the likelihood of a PTC being detected.

The remainder of this paper is organized as follows. In Section 2, we describe how ML was used to select PTCs from the 4LAC HBLs. In Section 3, we briefly review the Fermi spectral analysis of the PTCs and the EBL correction. In Section 4, we compare the corrected spectra with the sensitivities of the CTA and LHAASO. A discussion of the results and our conclusions are given in Section 5. Throughout the paper, we assume a Hubble constant $H_{0}=75\rm\ km\ s^{-1}\ Mpc^{-1}$, the matter energy density $\Omega_{M}=0.27$, the radiation energy density $\Omega_{r}=0$, and the dimensionless cosmological constant $\Omega_{\Lambda}=0.73$. We set the random seed as “120” in the algorithms with random processes.

ML techniques are popular in the field of astronomical data mining and data analysis (Ball & Brunner, 2010; Mirabal et al., 2012; Chiaro et al., 2016; Saz Parkinson et al., 2016; Salvetti et al., 2017; Baron, 2019; Kang et al., 2019a, b; Arsioli & Dedin, 2020; Fraga et al., 2021; Zhu et al., 2021). According to whether the classification of a sample is given, ML can be mainly divided into supervised machine-learning (SML) and unsupervised machine-learning (USML) techniques (Baron, 2019). SML approaches are usually used for classification and regression. USML approaches focus on clustering and dimensionality reduction by searching for potential relationships among the given samples. Numerous algorithms are applied to SML classification, including logistic regression (LR), decision trees, random forest(RF), support vector machines (SVM), NNs, and Bayesian networks (see, e.g., Feigelson & Babu 2012; Kabacoff 2015 etc.). In SML classification, the data set contains a certain number of objects, with each object containing several features and a label (Ballet et al., 2020). Features are usually observations that characterize the physical properties of the object, while the labels mark the classes. A classification algorithm model learns the relation between the features of known objects and their labels. The model can then be employed to evaluate a potential class of objects without clear classifications based on their features. The simple application of an SML classifier involves several steps, including data set preparation, model training, model optimization, model testing, and finally prediction.

Generally, for the binary classification of an unknown object between class A and class B, SML classifiers compute likelihood possibilities, $L_{\rm A}$ and $L_{\rm B}$, rather than give a classification directly. $L_{\rm A}$ or $L_{\rm B}$ correspond to the probability that the object belongs to A or B, respectively. Their relation can be expressed as $L_{\rm A}=1-L_{\rm B}$ (Chiaro et al., 2019). A binary “soft” classification can be achieved by selecting an $L_{\rm th}$ classification threshold in the range of $0\sim 1$. In standard SML binary classification, the boundary between class A and class B is clear and, after training, the classifier is always applied to an independent unknown data set. In the context of this paper, the two types of objects are HBLs detected in the TeV energy band (labeled as TeV) and HBLs which are not detected in the TeV energy band (labeled as non-TeV). However, if an object is labeled as “non-TeV”, this does not mean that it does not emit TeV $\gamma$-rays; it may not have been detected simply because of observational limitations. So, there will be TeV sources in the non-TeV data set. With soft SML classifiers, it is possible to train the model on these kinds of data sets to establish a probability space. In the probability space, we can search for TeV candidates among the non-TeVs that bear some resemblance to a TeV source. This makes it possible to train, optimize, and test an SML model in the same way as the standard process, so as to obtain a result with a higher level of confidence. The SML model can then be employed to compute the likelihood probabilities for each source to create a probability space, after which, the PTCs can be selected.

When evaluating a classifier, there are different metrics for performance evaluation that can be used, such as accuracy, recall, precision, and the receiver operating characteristic (ROC) curve (Baron, 2019). Accuracy is widely used, but simple classification accuracy is usually dominated by large samples from imbalanced data sets (Zhu et al., 2021). In this study, we used the balanced-accuracy (i.e., the average accuracy of positive and negative samples) rather than simple accuracy to evaluate the classifier performance.

Combining the results of multiple SML classification algorithms can reduce misclassifications (Zhu et al., 2021). In our approach, we combined two SML methods, SVM (Vapnik, 1999) and LR (e.g., Cox 1958; Walker & Duncan 1967). These are popular SML classifiers used in astronomy. Scikit-learn, originally known as sklearn (Pedregosa et al., 2011), is a complete ML integration Python package that makes it easy to construct various ML algorithm models. We created the SVM and LR models using the Scikit-learn package in a Python environment (version 3.8.5). The SML classification flowchart used in the context is shown in the Table 1.

The newest release of the Fermi $\gamma$-ray source catalog, 4FGL-DR2, contains GeV $\gamma$-ray spectral data for a 10 year period (2008-2018) in seven energy bins in the energy range of 50 MeV-500 GeV. By turning to more than 12 years of $\gamma$-ray observations provided by the Fermi Science Support Center (FSSC), we were able to analyze the $\gamma$-ray spectra of the PTCs obtained using the approach described above by utilizing more Fermi data than provided in the published catalogs. As the Fermi data were updated to P8R3 on 2008 November 12, we used the Fermi P8R3 data from 2008 October 1 to 2020 October 1 (mission elapsed time, MET, from 244548001 to 623253605). The photon events in the energy range from 100 MeV to 1 TeV were selected using the default data quality and a $90^{\circ}$ zenith angle. We used the corresponding instrument response functions for P8R3 SOURCE V3, the galactic interstellar emission model, gll$\_$iem$\_$v07 (i.e., gll$\_$iem$\_$v07.fits⁴⁴4https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/4fgl/Galactic_Diffuse_Emission_Model_for_the_4FGL_Catalog_Analysis.pdf), and the new isotropic spectral template (iso$\_$P8R3$\_$SOURCE$\_$V3$\_$v1.txt). Together with the Fermi science tool, Fermitools (version v11r5p3), the open-source Python package, Fermipy (Wood et al., 2017), was used to calculate the SEDs of the PTCs.

Each spectra was divided into 12 energy bins. For each energy bin, we provided the energy flux, and there was a 1$\sigma$ error when $TS(TestStatistic)\geq 9$ ($\textgreater 3\sigma$), with the energy flux upper limit having a confidence level of 95$\%$ when $0\textless TS\leq 9$. We removed the energy bin when $TS\leq 0$. Aside from the data points, Fermitools provided the best-fitting lines for characterizing the evolution of the spectrum when analyzing the Fermi spectra. The 12 year average $\gamma$-ray spectra of the 24 PTCs are shown in Figure 4. Six sources are displayed as a LogParabola (LP) spectrum in the $\gamma$-ray band and the other 19 PTCs have PL-type SEDs.

TeV $\gamma$-rays crossing interstellar space are attenuated by $\gamma+\gamma\to e^{+}+e^{-}$ through their interaction with EBL and CMB photons in wavelengths in the region of $0.1\sim 1000\ \mu\rm m$ (Gould & Schréder, 1967; Dwek & Krennrich, 2005). Depending on the redshift, EBL absorption may only have a strong effect on the flux above a few tens of GeV. However, the best-fitting lines were mainly dominated by Fermi’s low energy region, because there are higher photon counts in the low-energy region. The currently published Fermi source catalog has a detection energy upper limit of 1-3 TeV, although Fermi-LAT can detect higher energy photons. An increase in the detection energy band leads to a decrease in the effective detection area, an increase in the systemic uncertainty, and an insufficient high energy photon count. So, the Fermi spectra provide a good indication of the shape of the intrinsic source spectrum and can be extrapolated in the absorbed band.

We extended the Fermi best-fitting lines to the Fermi spectra to an energy of 100 TeV (see the black lines in Figure 4) assuming the spectra had an EBL optical optical depth, i.e., $\tau(E_{\gamma},z)$, of 0. Franceschini et al. (2008) and Franceschini & Rodighiero (2017) have provided an EBL optical depth table that takes into account the contribution of EBL photons. Using the EBL model from Franceschini et al. (2008) and Franceschini & Rodighiero (2017), we first calculated the EBL optical depth for all 24 TPCs from 20 GeV to 100 TeV, then corrected the Fermi spectra. The corrected best-fitting Fermi lines are shown in blue in Figure 4.

The upcoming CTA and LHAASO will form the next generation of TeV detectors. They will be characterized by an extremely high sensitivities and large FOVs. The CTA is a new generation IACT, which contains two arrays. The North array is located at $28^{\circ}.7$, while the South array is at $-24^{\circ}.7$. The sensitivity of the southern array is slightly higher than that of the northern one. Drawing upon these two arrays, the CTA can achieve all sky observations in the energy range of 20 GeV-10 TeV. The sensitivity of the CTA is affected by the zenith angle ($z_{\rm max}$), location, and the Earth’s magnetic field. Using the sensitivities provided by the CTA online calculator ⁵⁵5The sensitivities of CTA are available at https://www.cta-observatory.org/science/cta-performance/, for each PTC at a decl. of $b$, we obtained the CTA sensitivity curves for $\rm 5hrs\ 1yr^{-1}$ and $\rm 50hrs\ 1yr^{-1}$ as follows: [$b\geq 88^{\circ}.7$] $-$ sensitivity of the northern array at $z_{\rm max}=60^{\circ}$; [$58^{\circ}.7\leq b\textless 88^{\circ}.7$] $-$ sensitivity of the northern array at $z_{\rm max}=40^{\circ}$; [$2^{\circ}.7\leq b\textless 58^{\circ}.7$] $-$ sensitivity of the northern array at $z_{\rm max}=20^{\circ}$; [$-54^{\circ}.7\leq b\textless 2^{\circ}.7$] $-$ sensitivity of the southern array at $z_{\rm max}=20^{\circ}$; [$-84^{\circ}.7\leq b\textless-54^{\circ}.7$] $-$ sensitivity of the southern array at $z_{\rm max}=40^{\circ}$; [$b\leq-84^{\circ}.7$] $-$ sensitivity of the southern array at $z_{\rm max}=60^{\circ}$. LHAASO, which is located in Daocheng, Sichuan province, China, aims to detect cosmic rays and $\gamma$-rays with an energy higher than 30 TeV using three detectors: KM2A, WCDA, and WFCTA. LHAASO can naturally survey half of the all sky from decl. $-20^{\circ}$ to $80^{\circ}$ in all weather when the zenith angle is set at $50^{\circ}$ (Cao et al., 2019). For the PTCs in the LHAASO’s FOV, we plotted the sensitivity curve for one year of operation (Bai et al., 2019), and the sensitivity curve is shown as a black dotted line in Figure 4. The CTA sensitivity curves for $\rm 50\ hrs\ 1yr^{-1}$ and $\rm 5\ hrs\ 1yr^{-1}$ are shown by the red and green dotted lines, respectively.

We compared the results of the EBL-corrected PTCs with the CTA and LHAASO detection sensitivitied. There are 16 PTCs in the FOV of LHAASO, four of which are likely to be detected by LHAASO observations in light of the corrected Fermi spectrum. Out of the 24 PTCs, half are located in the northern sky area of the CTA northern array’s FOV. The energy spectra of all PTCs are above the sensitivity of the CTA’s $\rm 50\ hrs\ 1yr^{-1}$ observations, while only 13 PTCs are above the sensitivity of the CTA $\rm 5\ hrs\ 1yr^{-1}$ observations. Detailed information regarding the 24 PTCs is listed in Table 4. The distribution of 24 PTCs in the FOV of CTA and LHAASO under the Galactic coordinate system is shown in Figure 5.

In this study, we split the process of searching for TeV candidates in the 4LAC HBLs into two steps. First, we used SVM and RF algorithms to search for PTCs in the 4LAC HBLs by combining radio, optical, and GeV $\gamma$-ray data. This search revealed 24 PTCs that were above the minimum confidence standard in the SML probability space. We then analyzed PTC $\gamma$-ray spectra costructed from 12 years of Fermi observations and corrected them using an EBL model. Taking into account the sensitivity of the next generation CTA and LHAASO, we suggested four candidates for LHAASO observations and 24 candidates for CTA observations.

Current TeV detectors (e.g., IACTs and EAS arrays) are hindered by their limited sensitivities and FOVs, strong background interference, and their susceptibility to bad weather. The sensitivities of CTA and LHAASO are excepted to be approximately an order of magnitude higher than those of current detectors (Bernlöhr et al., 2013; Cui et al., 2014). CTA is excepted to yield good performance in the soft TeV band ($E\ \textless\ 10\ \rm TeV$), while LHAASO will focus on higher energy phenomena. The observation energy range of CTA and LHAASO may well therefore turn out to be complementary, enabling high detection sensitivity throughout the whole TeV energy range.

The samples used in the present work were limited to HBLs. However, there are also TeV FSRQs, TeV IBLs, and TeV LBLs in the overall population of TeV blazars. Emission models of blazars suggest that the peak frequencies of the two bumps in blazar SEDs are correlated (Abdo et al., 2010). A higher peak frequency in a synchrotron spectrum usually means there is a higher peak frequency for the high energy bump. This is why TeV blazars are dominated by HBLs. For example, 43 4LAC TeV HBLs account for 68.3$\%$ of the 4LAC TeV blazars, while the 283 Fermi HBLs account for approximately just 10$\%$. It can be seen in Figure 6 that the distributions of the TeVs and non-TeVs in terms of the peak frequency of synchrotron emissions are quite different. If sources displaying intermediate or low synchrotron peaks are accommodated within studies, a sample selection bias is inevitable (e.g., Richards et al. 2012; Luo et al. 2020; Zhu et al. 2021). On the other hand, in the 4LAC, the BL Lac objects with peaked frequency of synchrotron emission $\nu_{\rm syn}\textgreater 10^{15}\rm{Hz}$ are always recognized as HBLs, it is hard to clearly distinguish EHBLs from HBLs. Arsioli et al. (2020) proposed that the search for TeV blazars may benefit from considering HSP and EHSP as a whole. It is consistent with our sample selection criteria.

The best-fitting lines after EBL correction for most PTCs were a better match for the Fermi data points than the uncorrected fits. This also confirms that the part of the BL spectra of blazars that break at the GeV level can be attributed to EBL attenuation (Dwek & Krennrich, 2005). However, several sources did not match the best-fitting lines so well. For example, the spectra of 4FGL J0630.9-2406 (z=1.238) and 4FGL J1243.2+3627 (z=1.065) showed an obvious cutoff after the EBL correction was made that is not seen in the data points. This means that the result of the EBL correction-based TeV fluxes limitation for these two sources may be unacceptable. The excess of TeV $\gamma$-rays in high-redshift blazars remains an open issue. Several theories have been put forward to explain the spectra of hard $\gamma$-rays, such as axion-like particles (Simet et al., 2008; Sánchez-Conde et al., 2009), or Lorentz invariance violation (Protheroe & Meyer, 2000), but, to date, there is no definitive conclusion. In addition, the spectra of PTCs 4FGL J1120.8+4212 (z=0.124) and 4FGL J2323.8+4210 (z=0.059) has earlier breaks than the best-fitting lines in the GeV energy band. This manifested as a mismatch between the data points and the fitted line. It is possible that these resulted from the intrinsic nature of the spectra rather than EBL absorption. Liu & Bai (2006) and Liu et al. (2008) have indicated that the GeV break of blazars may result from the absorption in the board-line region. That thee break points of the PTC spectra we uncovered below 100 GeV provided evidence that the GeV breaks of a blazars do not have any single EBL origin.

The $\gamma$-ray spectra for some of the PTCs, such as 4FGL J0325.5-5635, 4FGL J1548.8-2250, and 4FGL J2042.1+2427, suggested the presence of a hard-soft-hard trend, which may correspond to a third component in addition to synchrotron radiation and inverse Compton scattering. The origin of the TeV peak remains another open issue. Exploring more HBLs with a TeV peak and investigating emission models, such as a secondary radiation model of the interactions between cosmic rays and galaxy background light in neighboring space, is planned for consideration in our future work.

We thank the anonymous referee for their very constructive and helpful comments and suggestions, which greatly helped us to improve our paper. We particularly thank Dr. C. Y. Yang from the Yunnan Observatory who provided us with many helpful comments and suggestions. We also thank the Fermi collaboration and Gaia collaboration for their data support. This work is partially supported by the National Natural Science Foundation of China (Grant Nos. 11763005 and 11873043) and the Science and Technology Foundation of Guizhou Province (QKHJC[2019]1290). This work is also supported by the Graduate Research Foundation of Yunnan Normal University. We would like to express our gratitude to EditSprings (https://www.editsprings.com/) for the expert linguistic services provided.