Exploring TeV candidates of Fermi blazars through machine learning

Data preprocessing is an essential step in the dataset before starting the training model. By manipulating or deleting the data, preprocessing of data can change natural features into a more available representation for the downstream model. It contains the following steps:

Data Discretization

In some of the datasets, one feature is represented by multiple columns. For example, in 4FGL-DR2, such as $Flux\_band$, 7 columns represent integral photon flux in each spectral band. We divide $Flux\_band$ as $Flux\_band\ 1$, $Flux\ band\ 2$, … $Flux\_band\ 7$, respectively.

Data cleaning

There may be some useless and missing data, which will reduce model performance. We discard error features, string features, and most missing value features.

Data standardization

To make all features dimensionless, we use sklearn.preprocessing.StandardScaler to convert all features into the standard normal distribution.

Data partition

We randomly divide the dataset into a 4:1 training set and testing set using python package StratifiedKFold.split of sklearn.model_selection, and each training / testing set preserves the percentage of the dataset for each class. We repeat the division 5 times with 5 different random seeds: 0, 1, 2, 3, 4 to get 5 training sets and 5 testing sets. We mark the training / testing sets as training 1 / testing 1, training 2 / testing 2, …, training 5 / testing 5.

We chose the SML method based on the following considerations:

1. 1.

   We request a highly interpretable model that is not overly dependent on computer performance, which outputs an explicitly linear boundary to separate TeV blazars and non-TeV blazars in the feature space.

2. 2.

   Once we have the physical properties of a blazar, that are, the features, we want to get the conditional probability that the blazar is a TeV one under the features.

There are many highly-developed SML classification methods, such as: Logistic Regression (Korsós et al. 2021, LR) Support Vector Machines (Hearst 1998, SVM), artificial neural network (McCulloch & Pitts 1943, ANN), random forest (Tin Kam Ho, 1995), etc. In this work, we choose the LR model over other methods like SVM, ANN and the likes for two main reasons. Mathematically, we assumed that the binary labels of blazars follow a Bernoulli distribution, and, as the logistic function is the expectation of a Bernoulli distribution, we believe it is the most natural way to map the linear combination of physical properties of the source into a Bernoulli distribution. Algorithmically, SVM, RF, and ANN are more complicated than LR, resulting in higher calculation costs and an overall loss in the ability to interpret the model, moreover the output of SVM needs to go through the LR model to obtain a conditional probability. The same holds true for ANN, although other models may be used to obtain conditional probability. In addition, ANN and RF cannot guarantee a linear boundary.

To illustrate LR, we introduce the odds ratio: $\frac{p}{1-p}$, where $p$ represents the probability of the positive event labeled 1. Then we further define the logit function (log-odds), which can be written as ${\rm logit}\,(p)={\rm log}\,\frac{p}{1-p}$. We can express a linear relationship between features and the log-odds:

here $y$ is the label; the positive event is marked as ‘1’, while ‘0’ for the negative events. $p(y=1|\mathbf{x}$) is the conditional probability that an individual labeled as ‘1’, and $\mathbf{w}$ is the weight. $p(y=1|\mathbf{x}$) (or $logistic$) is the inverse function of the ${\rm logit}$ function called the ${\rm logistic}$ function (or sigmoid function), , which can be expressed as:

where $logit(p)$ means $logit\,(p(y=1|\mathbf{x}))$.

For the implementation, we used scikit-learn (sklearn, Pedregosa et al. 2011), an ML library in python, which provides many useful tools. In particular, we used a sklearn API (Buitinck et al., 2013): sklearn.linear_model.LogisticRegression which implements the LR model, to train the training set. Eventually, we used the LR to classify the dataset by evaluating the $logit$ value of each source in the dataset and calculating the conditional probability $logistic$ that a source emits TeV radiation. Note that the number of TeV sources in our sample is small. To increase the weight of TeV sources, we set the weight parameter class_weight to balance in sklearn.linear_model.LogisticRegression.

A confusion matrix is applied to visualize the model performance. True positive (TP) represents the number of positive events that are correctly classified, true negative (TN) represents the number of negative events that are correctly classified, false positive (FP) represents the number of misclassified positive events, false negative (FN) represents the number of misclassified negative events, as shown in Fig. 2.

The accuracy ($\frac{TP+TN}{FP+FN+TP+TN}$) is one of the most common indexes for evaluating model performance. However, we have an extreme imbalance of the two classes, which means non-TeV sources dominate the results. We also want the model to identify as many TeV sources as possible while reducing the number of non-TeVs misjudged as TeVs. To meet these two requirements, we introduce larger Area Under the Curve (Fawcett 2006, $AUC$) as the metric to evaluate LR performance, which is based on the receiver operating characteristic curve (Fawcett 2006, ROC) with the true positive rate ($TPR$) on the Y-axis and the 1- false positive rate ($FPR$) on the X-axis, where $TPR=\frac{TP}{TP+FN}$ and $FPR=\frac{FP}{FP+TN}$. Furthermore, we introduce Youden’s J statistic (or Youden’s index, Youden 1950), which is $TPR-FPR$. The optimal $p_{\rm thre}$ corresponds to the maximum Youden index, which means that when we fix TPR, FPR reaches a minimum, which indicates that our screening of TeV candidates is very stringent. So we consider a non-TeV blazar as TeV one if its $logistic$ beyond the optimal $p_{\rm thre}$.

Now we can train the models on the 5 training sets and further optimize the model performance in two ways:

Feature selection

It may degrade model performance if all features are used, which is called the curse of dimensionality (Taylor, 2019). Data becomes sparse in high-dimensional space, resulting in no more extended similarities between data, which hinders efficient organization and data processing.

The sklearn.feature_selection.SequentialFeatureSelector (or SFS for short) of python could filter useless features, where the basic idea is to sequentially remove features from the full feature space until the new feature subspace contains the desired number of features.

Hyper-parameters optimization

Some parameters in the model cannot be obtained from the training model, which is often called hyper-parameters; they are not the $\mathbf{w}$ in Formula (1). Instead, we need to specify them manually before the training model.

The sklearn.grid_search.GridSearchCV (or GS for short) is helpful in the selection of optimal hyper-parameters, which performs a brute-force search on a list of values that we specify for hyper-parameters, and picks out the optimal hyper-parameters. For sklearn.linear_model.LogisticRegression, we optimize 6 hyper-parameters ⁸⁸8https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html: penalty, C, solver, tol, max_iter, l1_ratio.

The optimization strategy is inspired by the k-fold cross-validation (Ball & Brunner, 2010), which randomly divides the training set into k folds without replacement. Then k - 1 folds are used for the training model, and the remaining one fold, namely the validation set, evaluates the model performance. This process was repeated k times.

Based on k-fold cross-validation, we do both feature selection and hyper-parameter optimization at the same time using the nested cross-validation (Raschka, 2015; Varma & Simon, 2006), which is characterized by a nested loop. The dataset in the outer loop is divided into a training set and a testing set. The training set in the inner loop undergoes k-fold cross-validation selecting optimal hyper-parameters, which is performed by GS. Next, the SFS in the outer loop picks out the optimal features,and the testing set evaluates the model generalization. We implement nested cross-validation on 5 partitions of the dataset and get 5 models with optimal features and hyper-parameters.

However, Bias-variance trade-off is one of the fundamental rules of machine learning. The bias measures the deviation between the predicted value and the actual value on the known dataset, and the variance measures the performance of the model on the unknown dataset. Overfitting is a common problem in machine learning, and it is mainly due to high variance, making the performance on the testing set worse than that on the training set. Therefore, for each of the 5 partitions, we calculate the gap that equals $AUC$ of the training set minus $AUC$ of the testing set and selects one optimal model corresponding to the smallest positive gap.

Now we could build linear boundaries to distinguish TeV emitting sources from non-TeV emitting sources. When the $logistic$ of a source is greater than a threshold value ($p_{\rm thre}$), the LR determines the source as a TeV source, otherwise as a non-TeV one. A $p_{\rm thre}$ corresponds a $logit$ threshold value: $logit_{thre}$ (See Formula (2)). So we could express the linear boundary as $logit^{\prime}=0$, where $logit^{\prime}=logit-logit_{thre}$.

To better understand the process of predicting labels for the unknown datasets, we define some transition concepts. According to section  2.1, we call the dataset that has only undergone data discretization as the initial set. In contrast, the dataset that has undergone complete preprocessing is called the learning set. We get the optimal model from the learning set, and the optimal features are also determined. In the initial set, a subset containing the optimal features is called a prediction set. The learning set requires that all features have no missing values, but the prediction set only needs the full optimal parameters, containing more samples that will help to find more TeV candidates. Consequently, we could calculate the $logistic$ for each blazar in the prediction set, and consider the non-TeV blazars with $logistic\geqslant 80\%$ as high-confidence TeV candidates.

We used four catalogs: 4FGL-DR2 / 4LAC-DR2, 3FHL, 3HSP, and 2BIGB to complete our tasks. 4FGL-DR2, based on 10-year Fermi-LAT data, contains over 5700 sources, 3511 of which are AGNs (3436 are blazars). In the observer frame, 4LAC-DR2 includes synchrotron peak frequency ($\log\nu_{\rm p}^{\rm s}$), the corresponding intensity ($\log\nu_{\rm p}^{\rm s}f_{\nu_{\rm p}}^{\rm s}$), and redshift ($z$), which are critical supplements to 4FGL-DR2; 3FHL (Ajello et al., 2017) is a 7-year-based catalog of 1556 VHE sources; 3HSP (Chang et al., 2019) includes 2013 HSP blazars containing a critical feature: the peak flux in units of the faintest HSP blazar peak flux (2.5$\times$ 10^-12 $\>\mathrm{erg\cdot cm^{-2}\cdot s^{-1}}$) that has been detected in the TeVCat(Wakely & Horan, 2008): $FOM$ (Figure of Merit); 2BIGB (Arsioli et al., 2020) is the result of $\gamma$-ray likelihood analysis for 1160 3HSP sources yielded photon flux from 500 MeV to 500 GeV: $\rm F_{0.5-500GeV}$ ($F^{\rm ph}$).

TeVCat is an online interactive catalog containing VHE sources. As of now, it contains 251 sources detected mainly by IACTs and EAS arrays in the TeV energy region, and at least 91 are extragalactic sources. Among those sources, 80 blazars are detected by Fermi-LAT, thus we assume that blazars emitting TeV radiation are usually detected by Fermi-LAT. We cross-matched the four catalogs with TeVCat, respectively. We marked labels of their blazars in common with TeVCat as ‘1’, while the others were‘ 0’. According to the SML workflow introduced in the previous section, we got four groups of initial, learning, and prediction sets:

4FGL-DR2 / 4LAC-DR2

We supplement 4FGL-DR2 with $\log\nu_{\rm p}^{\rm s}$, $\log\nu_{\rm p}f_{\nu_{\rm p}}$ and $z$ from 4LAC-DR2, obtaining a new catalog: 4FGL-DR2 / 4LAC-DR2. Initial set: 3436 blazars including 80 TeV ones, with 88 features. Learning set: 861 blazars including 71 TeV ones, with 32 features. 688 / 173 blazars for each training / testing set. Prediction set: 1459 blazars including 73 TeV ones, with 5 features.

3FHL

Initial set: 1207 blazars including 74 TeV ones, with 49 features. Learning set: 506 blazars including 63 TeV blazars, with 18 features. 404 / 102 blazars for each training / testing set. Prediction set: 540 blazars including 67 TeV ones, with 2 features.

3HSP

Initial set: 2013 blazars including 57 TeV ones, with 13 features. Learning set: 1411 blazars including 53 TeV ones, with 5 features. 1128 / 283 blazars for each training / testing set. Prediction set: 1771 blazars including 54 TeV ones, with 2 features.

2BIGB

Initial set: 1160 blazars including 57 TeV ones, with 31 features. Learning set: 1040 blazars including 54 TeV ones, with 7 features. 832 / 208 blazars for each training / testing set. Prediction set: 1040 blazars including 54 TeV ones, with 3 features.

Processing the SML method we ascertain the result meeting condition of $logistic\geqslant$ 80% which gives 40 high-confidence candidates out of 150 4FGL-DR2 / 4LAC-DR2 candidates (see also Fig. 11, Tab. 8). Among the 40 high-confidence candidates, 24 sources are in common with 3FHL candidates, 11 with 3HSP candidates, and 14 with 2BIGB candidates. We report the main results of SML as follows:

4FGL-DR2 / 4LAC-DR2

$$logit^{\prime}=4.808\Gamma+2.809V_{F}+3.889\log f^{ph}_{\rm 7}-3.34z+0.857\log\nu_{\rm p}^{\rm s}+20.244,$$

(3)

where $logit^{\prime}=0$ is ideal for $p_{\rm thre}=$40%. The optimal LR model is built on training 3, with $AUC$ being 97% on training 3 and 96% on testing 3. 150 TeV blazar candidates are obtained from 1459 blazars in the prediction set, with $AUC$: 98%, $FPR$: 11%, and $TPR$: 99%. Only 1 TeV blazar is mistaken as a non-TeV one.

3FHL

$$logit^{\prime}=0.116\log f^{\rm ph}_{\rm 2}-0.628z+1.028,$$

(4)

where $logit^{\prime}=0$ is ideal for $p_{\rm thre}=$52%, the optimal LR model is built on training 5, with $AUC$ being 89% on training 5 and 87% on testing 5. 126 TeV candidates are obtained from 540 blazars, with $AUC$: 88%, $FPR$: 27%, and $TPR$: 88%. Only 8 TeV blazars are mistaken as non-TeV ones.

3HSP

$$logit^{\prime}=0.306FOM-3.861z-0.173,$$

(5)

where $logit^{\prime}=0$ is ideal for $p_{\rm thre}=$61%, the optimal LR model is built on training 5, with $AUC$ being 99% on training 5 and 94% on testing 5, 40 TeV candidates are obtained from 1771 blazars, with $AUC$: 98%, $FPR$: 2%, and $TPR$: 91%, only 5 TeV blazars are mistaken as non-TeV ones.

2BIGB

$$logit^{\prime}=-0.016E_{\rm P}+0.248FOM-4.395z-0.122$$

(6)

where $logit^{\prime}=0$ is ideal for $p_{\rm thre}=$48%, the optimal LR model is built on training 4, with $AUC$ being 97% on training 4 and 97% on testing 4. 83 TeV candidates are obtained from 1040 blazars, with $AUC$: 97%, $FPR$: 8%, and $TPR$: 96%, only 2 TeV blazars are mistaken as non-TeV ones.

The features of the 4 learning sets are reported in Tab. 1, where Col. (1) indicates the catalog information; Col. (2) gives how many columns the features contain, and each column represents a feature; Col. (3) and Col. (4) show the names and units of the features from the FITS version of the catalogs; Col. (5) is the description of the feature. The performance of the LR model on 4 learning sets is also shown in Tab. LABEL:LR_4FGL to Tab. 5, where Col. (1) shows the datasets name; Col. (2) the $AUC$, where the top is for training sets while the bottom for testing sets; Col. (3) $AUC$ difference between training sets and testing sets; Col. (4) is the $p_{\rm thre}$; Col. (5) is optimal features; Col. (6) is optimal hyper-parameters, note that the bold corresponds to minimum overfitting datasets, which represents the data and parameters of the optimal LR model.

We also visualize part of the results in figures: Fig. 3 shows feature selection; Fig. 4 is $AUC$ of the training / testing sets; Fig. 5 illustrates ROC curve, $AUC$, and $p_{\rm thre}$ of the prediction sets; Fig. 6 visualize the linear boundary of 3FHL, 3HSP, and 2BIGB. Fig. 7 is the confuse matrix.

The tactic is to correct the data by considering EBL first and then fitting SEDs on the intrinsic data. The SED contains two bumps, a low-energy one in the infrared to soft X-ray energy range due to synchrotron emission and a high-energy one in the region between hard X-ray to $\gamma$-ray that is associated with IC emission. To do so, we use an online tool ⁹⁹9http://tools.asdc.asi.it/SED/ provided by the space science data center (SSDC) of the Italian Space Agency to collect data from radio to TeV from many missions and experiments together with catalogs and archival data (Aharonian et al., 2001; Giommi et al., 2002; Amenomori et al., 2003; Aharonian et al., 2003, 2009; Daniel et al., 2005; Schroedter et al., 2005; Acciari et al., 2008, 2011; Godambe et al., 2008; Chandra et al., 2010, 2012; H. E. S. S. Collaboration et al., 2010; HESS Collaboration et al., 2013; Aliu et al., 2011, 2015; Bartoli et al., 2011, 2012; Archambault et al., 2013, 2014; Arlen et al., 2013; Abramowski et al., 2013, 2015; Biteau & Williams, 2015; Sharma et al., 2015). When we fit the SEDs, we discard the data bins whose flux error is bigger than flux upper limits. Note that there are TeV photons with zero error of the TeV candidates. Those TeV band photons are not certified by IACTs or EAS array (See Fig. 6). Likewise, we also discard data in the low radio energy range (${\rm log}\nu(\>\mathrm{Hz})<9$), since there is an asymmetry in the synchrotron bump of some sources. In this paper, we do not take into account the simultaneous observation of the data. EBL photons interact with VHE photons via pair production, a process that induces an attenuation of the observed gamma-ray spectra from blazars at cosmological distances, which may have a large impact on SED in IC bump. Assuming the EBL optical depth ($\tau\ (E,z)$) to depend on the photon energy ($E$) and redshift ($z$), Saldana-Lopez et al. (2021) provided a 2D grid of $\tau\ (E,z)$ ¹⁰¹⁰10https://www.ucm.es/blazars/eblin the range of 0.001 TeV $\leq E\leq$ 100 TeV, and 0 $\leq z\leq$ 6. Subsequently, the two bumps of the observed SED were fitted separately with a log-parabola as follows (Fan et al., 2016):

where $|2c|$ is the spectral curvature, $\log\nu_{\rm p}$ is the logarithm of the peak frequency, $\log\nu_{\rm p}f_{\nu_{\rm p}}$ is the logarithm of peak flux, and $\log(\nu f_{\nu})$ is the flux density.

we successfully fit the synchrotron bump for 40 TeV candidates and the IC bump for 20 TeV candidates and return the errors of the fitting $|2c|$, $\log\nu_{\rm p}$ (in units of $\mathrm{Hz}$) and $\log\nu_{\rm p}f_{\nu_{\rm p}}$ (in units of $\mathrm{erg\ cm^{-2}\ s^{-1}}$). The fitting process executes on scipy.optimize.curve_fit of python. The results are also listed in Tab. 6, in which Col. (1) gives the 4FGL name, Col. (2) to Col. (4) are parameters and errors for fitting synchrotron bump, Col. (5) to Col. (7) are parameters and errors for fitting IC bump, while the right side of the slash is the data without EBL-absorbed. For each high-confidence TeV candidate, the SEDs fitting results are shown in Fig. 11.

Being different from the all-sky scanning characteristics of space detectors, ground detectors can only scan part of the sky due to location constraints. IACTs need a dark sky for observations, and the field of view (FoV) of the current IACTs is small (3${}^{\circ}\textendash$ 5 ^∘). Besides, IACTs duty-cycle is restricted by the need to observe only during clear-sky, moonless nights, and, in addition, further constraints stem from the zenith angle: while the zenith angle increases, the sensitivity worsens and the energy threshold increases. This results in limiting the portion of the sky available for observation. At the same time, CTAO, as the next generation of IACT, consists of two arrays located in the Northern (28^∘ 45^′ N) and the Southern (24 ^∘ 41^′ S) hemispheres so that the FoV covers most of the sky, but not all, since CTAO suffers the same limitations affecting the current generation IACTs observations at high Zenith Angles ($Z$). Conversely, the EAS arrays can work under all weather conditions and have large FoV ($\sim{\rm 2}sr$), which can continuously monitor a significant fraction of the sky every day.

We compile the sensitivity curves of IACTs and EAS arrays from van Eldik et al. (2015); Aleksić et al. (2016); Cao et al. (2019); Abeysekara et al. (2017b, 2017) and online database: CTAO¹¹¹¹11https://www.cta-observatory.org/science/ctao-performance/, H.E.S.S. ¹²¹²12chrome-extension://efaidnbmnnnibpcajpcglclefindmkaj/https://www.mpi-hd.mpg.de/hfm/HESS/pages/home/proposals/sc_sens.pdf, VERITAS¹³¹³13https://veritas.sao.arizona.edu/about-veritas/veritas-specifications.. The FoV of these sensitivity curves are as follows: CTAO north (north site and $\mathrm{0}^{\circ}\leq Z\leq\mathrm{20}^{\circ}$), CTAO south (south site and $\mathrm{0}^{\circ}\leq Z\leq\mathrm{20}^{\circ}$), H.E.S.S. zenith ($Z\sim\mathrm{0}^{\circ}$), H.E.S.S. medium ($\mathrm{12}^{\circ}\leq Z\leq\mathrm{22}^{\circ}$), MAGIC low ($\mathrm{0}^{\circ}\leq Z\leq\mathrm{30}^{\circ}$), MAGIC medium ($\mathrm{30}^{\circ}\leq Z\leq\mathrm{45}^{\circ}$), VERITAS ($\mathrm{0}^{\circ}\leq Z\leq\mathrm{20}^{\circ}$), LHAASO ($Z\leq\mathrm{40}^{\circ}$ or $\mathrm{-11}^{\circ}\leq Dec\leq\mathrm{69}^{\circ}$), and HAWC ($\mathrm{-20}^{\circ}\leq Dec\leq\mathrm{60}^{\circ}$); while the exposure times are: 507 days for HAWC, one year for LHAASO, 25 hours for H.E.S.S. near the zenith, and 50 hours for the rests.

Since EAS arrays can observe all sources in the FoV simultaneously according to the source declination, whereas IACTs can only track one source at a time and require darkness sky, we employ two different strategies to infer the visibility of a given source. A source could pass the observation window for EAS arrays if its declination is within the most observable sky. On the other hand, for IACTs, we evaluate observation windows by 2 tools: the H.E.S.S. online visibility tool ¹⁴¹⁴14https://www.mpi-hd.mpg.de/hfm/HESS/pages/home/visibility/ and a python package: astroplan. Input the celestial coordinates of the source and the observation information of IACTs to these two tools. We can get visible months during darkness at a specific elevation angle (1 - $Z$), as shown in Fig. 8 and Tab. 7. The 2 tools have slight differences in the specific observable months. The observable months of 40 high-confidence TeV blazar candidates for the 2 tools are also shown in Tab. 7, where Col. (1) gives the source name; Col. (2) and Col. (7) give the observable months for IACTs and EAS arrays; The numbers represent the observable months, i.e. ‘1’ for January, ‘2’ for February, and so on.

We calculate the $\tau\ (E,z)$ for 40 TeV candidates with $E$ in the range of 0.001TeV to 100 TeV, and get the EBL-absorbed IC bumps. By comparing the 20 IC bumps and the sensitivity curves of IACTs and EAS arrays in the energy range of 1 TeV $\leq E\leq$ 100 TeV, we obtain 7 out of 40 high-confidence TeV blazar candidates that satisfy the two detectable criteria. The source number can be detected by the TeV facilities are (see Tab. 8 for detail): 6 CTAO (north or south) targets, 2 H.E.S.S targets, 5 MAGIC targets, and 1 VERITAS targets, and EAS arrays: 1 LHAASO targets, 0 HAWC targets. The SML and detectability results of the 40 high-confidence TeV blazar candidates are shown in Tab. 8, where Col. (1) gives the 4FGL name; Col. (2) and Col. (3) give the galactic coordinates (in degrees); Col. (4) synchrotron-peak frequency ($\log\nu^{s}_{p}$); Col. (5) redshift (z); Col. (6) flux density EBL-absorbed of sources at $1\>\mathrm{TeV}$ ($f_{1\>\mathrm{TeV}}$) which can be computed using Formula (7) the parameters of the IC bump; Col. (7) to Col. (10) $logistic$ in 4 catalogs, which is the likelihood of a source belonging to TeV sources given its optimal features Col. (11) the SED class in 4FGL-DR2 / 4LAC-DR2, where ‘bll’ stands for BL Lac, ‘fsrq’ for FSRQ, and ‘bcu’ for blazar candidate of uncertain type; Col. (12) to Col. (14) candidates compared with those of 3FHL, 3HSP, and 2BIGB, a candidate in common is marked as ‘Y’; Col. (15) a candidate comparing with those of other literature, ‘C & G’ represents the candidate is the same with Costamante & Ghisellini (2002), while ‘M’ for Massaro et al. (2013), ‘F’ for Foffano et al. (2019) ‘Chi’ for Chiaro et al. (2019); Col. (16) IC bump of candidates compare with the sensitivity of IACTs and particle detectors arrays in the range of $\mathrm{1TeV}\leq E\leq\mathrm{100TeV}$. ‘CN’, ‘CS’, ‘HZ’, ‘HL’, ‘ML’, ‘MM’, ‘V’, ‘L,’, and ‘H’, stand for CTAO north, CTAO south, H.E.S.S zenith, H.E.S.S low, MAGIC low, MAGIC medium, VERITAS, LHAASO, and HAWC, respectively.

We also want to know how far away the source becomes undetectable in the TeV energy band, i.e. the redshift upper limit of the source beyond which it will become undetectable. For this purpose, in the range of 0 $\leq z\leq$ 6, we compare the EBL-absorbed IC bumps for 20 blazars and the sensitivity curves of detectors, at $E=$ 1 TeV, 10TeV, and 100 TeV, respectively. Redshifts upper limits for the 20 blazars for each IACT and EAS array are listed in Tab. 9, where Col. (1) gives the 4FGL name; Col. (2) to Col. (10) are the redshift upper limits of the 20 blazars for each, where the values in each column are redshift upper limits at $E=$1 TeV, $E=$10 TeV, and $E=$100 TeV, respectively. The value ‘0’ means the IC bump can not be detected at the specific $E$; Col. (11) redshifts; Col. (12) is the detectability for Cherenkov detectors. One can see the VHE photons could survive from EBL easier than those with higher redshift. This can also be seen from Formula (3) $\sim$ (6), the negative coefficient of $z$ leads $logit$ to be inversely proportional to $z$. A larger redshift will result in a smaller logistic, indicating a source to be less likely to be a TeV candidate.

The SFS identified five features in 4FGL-DR2 / 4LAC-DR2 that are most critical to distinguish TeVs from non-TeVs: flux density, FOM, synchrotron peak frequency, redshift, spectral index, and variability, which are compatible with Lin & Fan (2016). Further considering X-ray emission and the $\log\nu_{\rm p}^{\rm s}$ shifts in behavior during a flare should be further analyzed in future work, sinces blazars that are not generally detected at TeV energies could be detected during flaring states. Therefore, it is possible that all blazars can emit the TeV photons. The blazars with TeV emission detected are only because their flux density in the TeV energy band is relatively high in specific periods, such as the flaring period. Note that 4FGL-DR2 includes features such as $Variability\_Index$, which mirrored the difference between the flux fitted in each time interval and the average flux over the entire catalog interval. Furthermore, $Frac\_Variability$ is the fractional variability computed from the fluxes in each year. The $logit^{\prime}$ of 4FGL-DR2 / 4LAC-DR2 contains the $Frac\_Variability$ making sure the LR is sensitive to the sources with synchrotron peak shifting during flaring states.

According to the classification proposed by Abdo et al. (2010) or Fan et al. (2016), sources with $\log\nu_{\rm p}^{\rm s}>$ 15 (or $\log\nu_{\rm p}^{\rm s}>$ 15.3) are HSPs. In TeVCat, if we considered only HSP/HSP BL Lacs, we would have ignored more than 11% of TeV sources. We presented here an alternative way, where we have not applied many initial cuts by filtering on SED class or other properties. To compare the performance of two criteria: HSPs ($\log\nu_{\rm p}^{\rm s}\geqslant 15.3$ Hz (Fan et al., 2016)) and LR model, we calculate the metric and confusion matrix on the testing sets of 4FGL-DR2 / 4LAC-DR2 and 3FHL. For 4FGL-DR2 / 4LAC-DR2, the $AUC$: 73% and $TPR$: 62% are both lower than those predicted by LR, up to 27 TeV blazars are misclassified as non-TeVs. For 3FHL, take off four sources without $\log\nu_{\rm p}^{\rm s}$ data, the $AUC$: 68% and $TPR$: 63% are also lower than those predicted by LR, up to 25 TeV blazars are misclassified as non-TeVs. The relevant confusion matrix is shown in Fig. 9 and Fig. 10. The results illustrated that our selection criterion is more effective.

So far, TeV candidate selection techniques privileged targets with high predicted TeV fluxes. In our analysis, instead, we identified targets that, according to their broadband features, have high chances to become promising TeV targets in specific activity states and that may have been missed so far due to their current status or to the lack of the proper instruments and observing strategies. As a data-driven algorithm, ML is inevitably sensitive to the sample composition, the ratio of the training set and test set, the value of k in cross-validation, etc. Our result is to ensure its generalization based on existing data. There is a potential bias in the fact that a significant fraction of TeVCat sources are Fermi-LAT detected blazars, as Fermi-LAT catalogs are browsed to list potential candidates for IACTs, and because Fermi-LAT monitoring is used to trigger Target of Opportunity programs. Also, our observation proposal does not consider the weather or available time of detectors. A feasible observation proposal should be more reliable and agreed upon by the observatory that considers more influencing factors.

Detecting extragalactic photons in the TeV band needs to face two major limitations, which are EBL and the performance of detectors (sensitivity, effective area, FoV, etc). The relationship between the improvement of detector sensitivity and the detection of more TeV sources is far from simple, because if it is true that an improvement in sensitivity indeed helps detect more sources but, since the attenuation rapidly increases with energy and redshift, the overall effect is rather mild (see for example Gilmore et al. 2012). The redshifts of the 81 blazars listed in TeVCat are taken from the four catalogs and literature, noting that for 4FGL J2243.9+2021 (or RGB J2243+203) we took the upper limit of the redshift range of 0.75$\leq z\leq$1.1 (Sahu et al., 2019). Then we calculate the $\tau\ (E,z)$ fixing $E=$ 1 TeV and we evaluate $\tau\ (E,z)\leq$ 13.1 for all 81 blazars. The results are also listed in Tab. 10, where Col. (1) and Col. (2) give the TeVCat name and 4FGL-DR2 name; Col. (3) and Col. (4) are the celestial equator coordinates (in degrees); Col. (5) the redshift; Col. (6) the $\tau\ (E,z)$ provided by Saldana-Lopez et al. (2021) at $E=$ 1 TeV.

All in all, the validation of TeV candidates is counting on IACTs and EAS arrays, while each of the IACTs and EAS arrays has advantages and disadvantages. IACTs have better energy and angular resolution. Besides, the energy threshold of IACTs could reach the GeV band while EAS arrays are generally sensitive above $\sim$ 10 TeV, but at $\sim$ 30 TeV, LHAASO sensitivity is comparable to that of CTAO, and LHAASO at $\sim$ 100 TeV is the most sensitive instrument (Cao et al., 2019). On the other hand, EAS array observation time is not limited to moonless nights and has a much larger FoV, and the exposure time of IACTs is typically quoted for 50 hr whereas for 1 or 5 years of EAS arrays. Although IACTs are better at capturing transient events, larger FoVs and longer exposure time make up the gap for the EAS array. Overall, for IACTs and EAS arrays, the differences in their operation modes and performance parameters allow them to complement each other. Furthermore, The sensitivity of LHAASO and the next-generation IACT array: CTAO, has been improved by an order of magnitude compared with the current IACTs and EAS arrays. They have pushed the detection limit of the TeV band to higher than 100 TeV (The CTA Consortium et al., 2011; Cao et al., 2019), enabling us to detect the whole energy range of the TeV sky in the future.

In this paper, we implement a machine learning algorithms called Logistic Regression to search the TeV blazar candidates from 4FGL-DR2 / 4LAC-DR2, 3FHL, 3HSP, and 2BIGB. Furthermore, we filter out potential observation targets for IACTs and EAS arrays. The main conclusions are enumerated below:

1. 1.

   For the 4 catalogs, Logistic Regression provides an empirical formula to find blazars that could be detected at TeV energies with $logit^{\prime}\geqslant\rm 0$.

   For 4FGL-DR2 / 4LAC-DR2:

   $$logit^{\prime}=4.808\Gamma+2.809V_{F}+3.889\log f^{ph}_{\rm 7}-3.34z+0.857\log\nu_{\rm p}^{\rm s}+20.244,$$

   with $AUC=$ 98%, and the $FPR=$ 11%, $TPR=$ 99%, and the $p_{\rm thre}=$ 40%. 40 out of 150 non-TeV blazars have $logistic\geqslant$ 80% and are thus expected to be high-confidence candidates.

   For 3FHL:

   $$logit^{\prime}=0.116\log f^{ph}_{\rm 2}-0.628z+1.028,$$

   with $AUC=$ 88%, $FPR=$ 27%, $TPR=$ 88%, and the $p_{\rm thre=}$ 52%. 24 out of 126 TeV candidates are common with the high-confidence ones of 4FGL-DR2 / 4LAC-DR2;

   For 3HSP:

   $$logit^{\prime}=0.306FOM-3.861z-0.173,$$

   with $AUC=$ 98%, $FPR=$ 2%, $TPR=$ 91%, and $p_{\rm thre=}$ 61%. 11 out of 40 TeV candidates are common with the high-confidence ones of 4FGL-DR2 / 4LAC-DR2;

   For 2BIGB,

   $$logit^{\prime}=-0.016E_{\rm P}+0.248FOM-4.395z-0.122,$$

   with $AUC=$ 97%, $FPR=$ 8%, $TPR=$ 96%, and $p_{\rm thre=}$ 48%. 14 out of 83 TeV candidates are common with the high-confidence ones of 4FGL-DR2 / 4LAC-DR2;

2. 2.

   For the 40 high-confidence candidates, we independently fit the two bumps of the SEDs, infer the visibility in the sky. For EAS arrays and IACTs in 2023, 7 candidates are potential targets: 6 CTAO targets, 2 H.E.S.S targets, 5 MAGIC targets, 1 VERITAS, 1 LHAASO targets, and 0 HAWC targets, respectively.

3. 3.

   We get 2 common sources with Costamante & Ghisellini (2002), 9 ones with Massaro et al. (2013), 4 ones with Foffano et al. (2019), and none with Chiaro et al. (2019), respectively.