Hunting Very High-Energy (\(>\)​100 GeV) Emitting High-Synchrotron Peaked Blazars


Abstract

Very-high energy (VHE; \(>\)​100 GeV) \(\gamma\)-ray emission originates via some of the most extreme particle acceleration processes in the universe. Considering beamed active galactic nuclei, i.e., blazars, only a small fraction, mainly high synchrotron peak BL Lacs, have been detected in the VHE band with the ground-based Cherenkov telescopes. We utilized \(\sim\)​16 years of Fermi-Large Area Telescope (LAT) observations in the 0.1\(-\)​2 TeV energy range to systematically search for potential VHE emitters in a sample of high synchrotron peaked (\(\nu^{\rm peak}_{\rm syn}>10^{15}\) Hz) BL Lac sources. We identified, for the first time, 92 VHE emitting blazars at \(\geq 5\sigma\) confidence level. A significant VHE emission was also detected from 52 objects, which have been previously reported to be a VHE blazar. Comparing with the general blazar population, these VHE emitting blazars are found to be located at low redshifts (mean \(z=0.2 \pm 0.1\)) and exhibit bright synchrotron emission (\(\log F^{\rm peak}_{\rm syn}=-11.2 \pm 0.4\), in erg cm\(^{-2}\) s\(^{-1}\)). We also investigated the coincidence of VHE photon arrivals with the source activity states and found that Fermi-LAT has detected VHE photons during both quiescent and elevated activity epochs. These VHE emitting blazars represent promising targets for current and next-generation ground-based Cherenkov telescopes, and provide valuable laboratories for probing particle acceleration in relativistic jets, testing multi-messenger connections, and constraining extragalactic background light models.

1 Introduction↩︎

Blazars are a subclass of radio-loud active galactic nuclei (AGN) distinguished by relativistic jets oriented closely along the line of sight to the observer [1], [2]. This alignment leads to strong Doppler boosting of the jet emission, resulting in brightness enhancement and extreme variability across the entire electromagnetic spectrum, with timescales ranging from years down to minutes. Their broadband spectral energy distributions (SEDs) in \(\nu\,-\nu\,f_{\nu}\) representation exhibit a characteristic double-humped structure, indicative of two distinct non-thermal emission processes. The low-energy hump, extending from the far-infrared to the soft X-ray regime, is attributed to synchrotron radiation from relativistic electrons in the jet. The high-energy component, extending from hard X-rays to \(\gamma\)-rays and occasionally into the TeV range, is generally ascribed to inverse Compton scattering of relativistic electrons by photons internal or external to the jet or to hadronic emission processes [3][9].

The synchrotron peak in the SED varies across blazars and serves as a key indicator of the jet’s particle acceleration efficiency [10]. Based on the synchrotron peak frequency (\(\nu^{\rm peak}_{\rm syn}\)), blazars are typically categorized into three subclasses: low-synchrotron-peaked (LSP; \(\nu^{\rm peak}_{\rm syn} < 10^{14}\) Hz), intermediate-synchrotron-peaked (ISP; \(10^{14} \lesssim \nu^{\rm peak}_{\rm syn} < 10^{15}\) Hz), and high-synchrotron-peaked (HSP; \(\nu^{\rm peak}_{\rm syn} > 10^{15}\) Hz) blazars [3]. A subset of HSP sources visibly showcase synchrotron peaks extended into the keV regime – commonly referred to as extreme blazars (EHSP) [11][13]. HSP and EHSP BL Lac objects stand out for their characteristically hard \(\gamma\)-ray spectral profiles, effective particle acceleration mechanisms, and frequent detection in the very high-energy (VHE; E \(\geq\)​100 GeV) \(\gamma\)-ray band [14], [15].

Observationally, HSPs and EHSPs are prime candidates for VHE observations. This arises because relativistic electrons in HSPs are accelerated to higher energies than in LSPs, enhancing the likelihood of inverse Compton upscattering to produce VHE \(\gamma\)-rays [16][21]. An updated list of known VHE detected sources, many of which are HSPs, is maintained in the web-based TeVCat1 catalog [22]. BL Lac objects with radiative output upto TeV energies serve as powerful probes of extreme astrophysical environments, offering insights into particle acceleration mechanisms [23], [24], potential associations with high-energy neutrinos and cosmic rays [25], [26], and enabling measurements of the extragalactic background light [27][29] and constraints on key aspects of \(\gamma\)-ray cosmology, such as the intergalactic magnetic field [30][32].

Identifying potential VHE \(\gamma\)-ray candidates has become increasingly important with the rise of multimessenger astrophysics and the improved sensitivity expected from the next generation of VHE observatories for a more targeted survey of sky. However, current VHE astronomy is constrained by a limited number of known TeV-emitting sources – as reflected in the relatively sparse entries of the TeVCat catalog and by observational challenges inherent to imaging atmospheric Cherenkov telescopes (IACTs) such as Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescope, Very Energetic Radiation Imaging Telescope Array System (VERITAS), High Energy Stereoscopic System (H.E.S.S.), and Whipple observatory. These instruments usually have low-duty cycle, as observations can only be conducted during clear, moonless nights, limiting data collection to roughly 10–15% of the year. Additionally, weather conditions like humidity and atmospheric instability can degrade data quality. IACTs also have a narrow field of view, restricting their ability to conduct wide-area surveys and necessitating targeted monitoring strategies [33].

In contrast, the Fermi Large Area Telescope (Fermi-LAT) has been performing continuous all-sky monitoring since 2008 in the energy range of 100 MeV to 2 TeV. Although its sensitivity decreases at higher energies, it benefits from a large field of view and continuous exposure, which helps mitigate background contamination and makes it well-suited for identifying promising candidates in the VHE domain (\(E>100\,\rm{GeV}\)). The overlap between Fermi-LAT’s upper energy range and the VHE regime offers a complementary approach for uncovering new VHE-emitting blazars, thus providing potential targets for follow-up with current and future ground-based VHE instruments [34], [35].

In this work, we take advantage of Fermi-LAT’s high-energy sensitivity and its extensive all-sky monitoring since 2008 to search for potential VHE-emitting blazars. Given that HSP blazars are the best candidates for VHE detection, we use the population of bonafide HSPs as a seed sample for further investigation. The strategy for candidate identification and sample selection is outlined in Section 2. Details of the data reduction and analysis procedures are provided in Section 3. The list of newly identified candidates is presented in Section 4, and the results of the population-level study are discussed in Section 5.

2 Sample selection↩︎

The fourth data release of the fourth Fermi-LAT \(\gamma\)-ray source catalog includes 7194 objects detected over 14 years of observations in the 50 MeV to 1 TeV energy range [36], [37]. A subset of these sources are recognized as VHE emitters, with counterparts listed in the TeVCat online catalog [22]. Currently, the identification of VHE sources largely depends on follow-up observations triggered by detections from ground-based and space-based multi-wavelength facilities. Given the continuously growing Fermi-LAT dataset and its all-sky monitoring capabilities, it presents a valuable instrument to identify new VHE candidates and support future efforts in VHE \(\gamma\)-ray astronomy.

Among the most promising VHE source classes are HSP blazars [38], which are considered the best targets for follow-up VHE observations [21], [39]. The most comprehensive and up-to-date sample of these sources is provided by the 3HSP catalog, which includes 2013 HSP blazars [40].

The primary objective of this work is to identify potential VHE-emitting sources among \(\gamma\)-ray detected HSP sources, using the 4FGL-DR4 catalog. To achieve this, we cross-matched the 4FGL-DR4 and 3HSP catalogs using a search radius of 5 arcseconds, leading to the selection of 1004 HSP sources. Additionally, recognizing that about 30% 4FGL-DR4 sources lack multiwavelength counterparts, we also considered 41 3HSP blazars that lie within the 95% positional uncertainty regions of unassociated Fermi sources, and hence likely to be the low-frequency counterpart of the unassociated \(\gamma\)-ray sources. This results in a parent sample of 1045 \(\gamma\)-ray emitting 3HSP sources, which we further studied for signatures of VHE emission.

3 Data Reduction↩︎

We utilized \(\sim\)​16 years of Fermi-LAT observations conducted from MJD 54683 to 60587 (2008 August 5 to 2024 October 4). We performed the standard binned likelihood method using fermiPy [41] with fermitools (v 2.2.0). The data was analyzed using the latest instrument response function P8R3_SOURCE_V3 in the energy range 100 GeV\(-\)​2 TeV and utilized photons within \(2^{\circ}\) around the source location, considering the PSF of the telescope above 100 GeV [42]. The SOURCE class events (evclass\(=\)​128) were spatially binned with \(0.05^{\circ}\) per pixel and divided into ten logarithmically spaced bins per energy decade. Additionally, to ensure high-quality data corresponding to the good time intervals, we employed the filter expression recommended by the LAT team DATA_QUAL \(>\) 0 && LAT_CONFIG==1. We also adopted a zenith angle cut of \(105^{\circ}\) to eliminate contamination from secondary \(\gamma\)-rays from Earth’s albedo.

The cleaned events data was then modeled by considering all 4FGL-DR4 sources lying within \(3^{\circ}\) of the source location. Bright sources with Test statistics (TS) \(>25\) were allowed to vary, whereas sources with TS\(<1\) were removed from the model. We additionally accounted for contributions from galactic and isotropic background by using gll_iem_v07 and iso_P8R3_SOURCE_V3_v1 template [43]. Considering the low photon statistics, the \(\gamma\)-ray emission from the source of interest was modeled with a power law spectral shape. We also computed the time of arrival of photons having more than \(95\%\) probability of association with the source using the tool gtsrcprob. For sources with TS\(>25\), we repeated the likelihood fitting using the EBL model provided by [44] to estimate the intrinsic, i.e., EBL-attenuation corrected, spectral parameters.

4 Results↩︎

The Fermi-LAT data analysis of 1045 sources led to the identification of 144 \(\gamma\)-ray sources detected at \(\gtrsim 5\sigma\) confidence level, i.e., TS\(\geq\)​25. Additionally, 187 sources were detected with \(9\leq{\rm TS}<25\). Cross-referencing the parent sample of 1045 objects with the TeVCat, we found that 60 of them have already been reported as VHE emitters in previous works. Of these, 52 are significantly detected in our analysis with TS \(\geq\)​25, while 7 sources have \(9\leq{\rm TS}<25\) (4FGL J1103.6\(-\)​2329 (TS=24), 4FGL J1518.0\(-\)​2731 (TS=22), 4FGL J1230.2+2517 (TS=20), 4FGL J1442.7\(+\)​1200 (TS=20), 4FGL J1958.3\(-\)​3010 (TS=17)(TS=18), 4FGL J0847.2\(+\)​1134 (TS=16), 4FGL J0013.9\(-\)​1854 (TS=16)). The only known VHE emitting blazar missing in our list is the HSP source 4FGL J0152.6+0147 (associated with RGB J0152+017, \(z=0.08\)) with a TS of 2.5. The source was detected in VHE energy range only in 2007 October and November by H.E.S.S., i.e., before the launch of the Fermi satellite [45]. Therefore, it is likely that during the period of the Fermi-LAT operation, the source remained in quiescence and hence was not detected with the satellite. All in all, we have identified 92 new VHE emitting HSP blazars using the \(Fermi\)-LAT data. Their spectral parameters estimated from the Fermi-LAT data analysis are provided in Table [tab:new95candidates951], whereas those for already known VHE emitters are listed in Table [tab:known95VHE95candidates]. Sources detected at low significance (\(9\leq{\rm TS}<25\)) are tabulated in the appendix (Table [tab:ts95betwen95995and9525]).

lccccccccc J0014.7\(+\)​5801 & 27 & \(0.77 \pm 0.50\) & \(3.33 \pm 1.27\) & \(-\) & \(-\) & \(-\) & 3 & 160.8 & 54957.0
J0015.6\(+\)​5551 & 39 & \(1.57 \pm 0.83\) & \(2.80 \pm 0.77\) & 40 & \(1.53 \pm 0.74\) & \(2.01 \pm 0.88\) & 4 & 230.0 & 58384.7
J0022.0\(+\)​0006 & 41 & \(1.61 \pm 0.86\) & \(3.34 \pm 1.06\) & 41 & \(1.56 \pm 0.77\) & \(2.36 \pm 1.24\) & 3 & 353.8 & 57494.7
J0047.9\(+\)​5448 & 26 & \(1.28 \pm 0.91\) & \(2.16 \pm 0.79\) & \(-\) & \(-\) & \(-\) & 4 & 605.9 & 56692.3
J0051.2\(-\)​6242 & 30 & \(0.45 \pm 0.33\) & \(8.35 \pm 4.07\) & 30 & \(0.46 \pm 0.34\) & \(8.08 \pm 4.28\) & 3 & 151.8 & 58968.8

llllllllll J0033.5\(-\)​1921 & 36 & \(1.49 \pm 0.93\) & \(3.19 \pm 1.14\) & \(-\) & \(-\) & \(-\) & 3 & 170.0 & 54956.5
J0035.9\(+\)​5950 & 396 & \(8.54 \pm 1.75\) & \(2.97 \pm 0.34\) & 398 & \(8.04 \pm 1.43\) & \(1.11 \pm 0.44\) & 29 & 655.3 & 57381.2
J0136.5\(+\)​3906 & 82 & \(1.50 \pm 0.62\) & \(4.29 \pm 1.22\) & \(-\) & \(-\) & \(-\) & 7 & 250.6 & 56683.5
J0214.3\(+\)​5145 & 61 & \(2.62 \pm 1.08\) & \(3.01 \pm 0.67\) & 61 & \(2.61 \pm 1.05\) & \(2.86 \pm 0.68\) & 5 & 179.2 & 58240.6
J0232.8\(+\)​2018 & 30 & \(1.81 \pm 1.14\) & \(2.53 \pm 0.82\) & 31 & \(1.75 \pm 1.03\) & \(2.00 \pm 0.89\) & 3 & 308.1 & 59508.4

There are 456 HSP sources from which at least one VHE photon was detected by the Fermi-LAT. All 144 significantly detected objects (with TS \(\geq\)​25) have been found to emit two or more VHE photons with \(\geq 95\%\) association probability, thus robustly confirming them to be VHE emitting blazars [46]. We also found 312 sources that emitted at least one VHE photon; however, their detection significance is low (TS \(<\) 25). This includes 187 sources with TS values between 9 and 25, and 125 sources with TS \(<\)​9. Though their VHE detection cannot be claimed at high significance, these objects can be considered as candidate VHE-emitting blazars. In the appendix, we provide the relevant VHE photon information (Table [tab:ts95betwen95995and9525] and 1).

Out of 41 unassociated \(\gamma\)-ray sources, four were detected with TS\(\geq\)​9, and one VHE photon with \(\geq\)​95% association probability was detected from ten sources (see Table [tab:unassociated95sources]). This finding suggests these HSP blazars to be the likely counterpart of the unassociated VHE emitting \(\gamma\)-ray sources.

VHE \(\gamma\)-rays traveling cosmological distances undergo pair production when interacting with EBL photons, leading to energy-dependent flux suppression that becomes increasingly significant at redshifts \(z \gtrsim 0.1\), especially for \(\gamma\)-rays with energies above 100 GeV [47]. We applied EBL absorption correction to recover the intrinsic spectra utilizing the spectroscopic redshift for the identified candidates [48], [49]. EBL-corrected spectral parameters are provided in Tables [tab:new95candidates951] and  [tab:known95VHE95candidates]. A significant fraction of these sources exhibit EBL-corrected photon indices below 2.0, indicating them to have intrinsically hard \(\gamma\)-ray spectrum [14].

5 Discussion↩︎

The origin of VHE emission in blazars remains a relatively unexplored problem in high-energy astrophysics. Although the number of confirmed VHE blazars is limited, they offer valuable insight into extreme particle acceleration and the physical conditions that enable such high-energy processes. Expanding the sample of VHE candidates is particularly important in view of the upcoming Cherenkov Telescope Array Observatory (CTAO), which will greatly enhance detection sensitivity in this regime. A well-defined candidate list is therefore essential for guiding targeted observations and advancing our understanding of extreme particle accelerators. In this work, we report, for the first time, the identification of 92 VHE \(\gamma\)-ray emitting blazars and confirm VHE detections of 52 sources previously identified in this energy range.

Figure 1: Left: Distribution of synchrotron peak frequency for the 4HSP-Fermi non-VHE population (blue) and VHE sources identified in this study (yellow). Middle: Redshift distribution for the same populations. Right: Distribution of synchrotron peak flux for the same populations.)

5.1 Comparison with non-VHE Emitting HSP Blazars↩︎

To understand the origin of the VHE emission, we compared the brightness, redshift, and synchrotron peak properties of the identified VHE sources and non-VHE-detected HSP sources using data from the 4FGL-DR4 catalog (see Figure 1). The synchrotron peak frequency distribution of the VHE emitters (\(\log \nu^{\rm peak}_{\text{syn}} = 16.5 \pm 0.9\), in Hz) peaks at slightly higher frequency compared to the non-VHE 3HSP–Fermi sources (\(\log \nu^{\rm peak}_{\text{syn}} = 16.1 \pm 0.8\), in Hz) though the dispersion is large. Furthermore, most of the VHE blazars are located in the nearby universe, with a mean redshift of \(0.2 \pm 0.1\), which is lower than the non-VHE population (\(0.3 \pm 0.2\)), though with a large spread. Additionally, the synchrotron peak flux (\(\log \nu F_{\nu,\text{syn}}\)) for the VHE emitters is systematically higher, peaking at \(-11.2 \pm 0.5\) (in erg cm\(^{-2}\) s\(^{-1}\)), compared to \(-11.9 \pm 0.4\) (in erg cm\(^{-2}\) s\(^{-1}\)) for non-VHE sources, hinting at a potential brightness-driven bias in the VHE detection. This is because, considering one-zone leptonic models, a brighter synchrotron emission peaking at X-rays corresponds to the inverse Compton spectrum peaking at hundreds of GeV energies, making the source brighter at VHE energies. Indeed, the detection of hard \(\gamma\)-ray spectra at MeV-GeV energies from the HSP blazars supports this hypothesis [50].

Figure 2: Top: Variation of the photon index with energy flux. Middle: The variability index versus the power law index. Bottom: This plot shows the variability index versus energy flux for 4FGL-DR4 sources (orange circles), 3HSP–Fermi sources (red circles), new VHE emitters (empty blue diamonds), and known VHE sources (filled blue squares).

5.2 Comparison with General Gamma-ray Emitting Blazar Population↩︎

In Figure 2, we show the variations of the power law photon index, variability index, and energy flux with respect to each other, by adopting these values from the 4FGL-DR4 catalog. Such comparisons enable us to explore the distribution of VHE-emitting blazars within the spectral and variability parameter space and reveal clear distinctions from the broader \(Fermi\)-LAT detected blazar population. As illustrated in Figure 2, newly identified VHE blazars preferentially exhibit hard photon spectra (median photon index \(1.8 \pm 0.1\)) and moderate-to-high energy fluxes in the 0.1\(–\)​100 GeV range (\(\gtrsim 1.1 \times 10^{-12}\,\rm erg\,cm^{-2}\,s^{-1}\)), whereas the general 4FGL-DR4 sample includes a substantial fraction of softer, fainter sources. Both newly identified VHE sources and already known VHE emitters share a similar parameter space in all the plots.

Interestingly, a majority of the VHE emitting sources tend to display relatively low-to-moderate variability indices (median value of \(35\)) for sources of comparable flux. This observation can be explained as follows: the detection of a hard \(\gamma\)-ray spectrum (0.1\(-\)​100 GeV photon index \(<\)​2) implies that the inverse Compton peak lies at energies \(>\)​100 GeV. Therefore, the 0.1\(-\)​100 GeV emission is mainly produced by the low-energy electron population whose cooling time is longer, hence lesser variability, compared to high-energy electrons producing \(>\)​100 GeV to TeV radiation [51], [52].

Figure 3: Cosmic \gamma-ray horizon: Highest-energy photons from identified VHE-emitting HSPs plotted against their spectroscopic redshift, based on Fermi-LAT analysis in the 0.1–2 TeV energy range. Sources are color-coded according to the EBL optical depth values derived from [53].

5.3 Cosmic Gamma-ray horizon↩︎

In addition to the jet energetics and spectral properties of the source, the detectability of VHE \(\gamma\)-ray sources is significantly influenced by attenuation due to the EBL, which limits \(\gamma\)-ray propagation over cosmological distances. To evaluate this effect for the newly identified VHE emitting HSP blazars, we plotted the energy of the maximum energy photon as a function of the blazar redshift in Figure 3. We also show the cosmic \(\gamma\)-ray horizon, i.e., EBL optical depth (\(\tau_{\rm EBL}=1\)) with a dashed line following the EBL model proposed by [53].

The majority of VHE emitting blazars lie below or near this horizon, indicating that they are situated in regions where the Universe remains transparent to VHE \(\gamma\)-rays. This reinforces their viability as detectable targets for current and future ground-based \(\gamma\)-ray observatories. The highest redshift source in our analysis - 4FGL J0506.0\(+\)​6113 [40] has the highest energy photon \(\sim148\) GeV and lies below the horizon. Fermi-LAT, with its limited sensitivity, detected \(>\)​1 TeV photons from 5 sources (highlighted in bold in Table [tab:known95VHE95candidates]). In Figure 3, a few sources appear to be located with EBL opacity \(\tau >2\). These sources, although fewer, may serve as valuable probes of the \(\gamma\)-ray opacity of the Universe. We briefly comment on these objects below.

4FGL J0035.9\(+\)​5950 (1ES 0033\(+\)​595): This object is a known VHE emitter with a synchrotron peak frequency of \(\log \nu^{\rm peak}_{\rm syn}=18.2\) (in Hz). Two spectroscopic redshift measurements are reported: \(z=0.467\) [54] and \(z=0.086\) [55]. In our work, we adopted \(z=0.467\) since this redshift measurement was based on a higher signal-to-noise ratio optical spectrum [54]. Fermi-LAT registered 29 photons above 100 GeV, with the highest-energy photon of 655.3 GeV.

4FGL J0349.4\(-\)​1159 (1ES 0347\(-\)​121): This extreme blazar is a known VHE emitter at \(z=0.19\) [48] and \(\log \nu^{\rm peak}_{\rm syn}=18\) (in Hz). The energy of the highest energy photon detected by the Fermi\(-\)LAT is 1250.2 GeV, corresponding to an EBL optical depth of \(\tau_{\rm EBL} = 2.24\).

4FGL J0507.9\(+\)​6737 (1ES 0502\(+\)​675): This object is also a known VHE emitting BL Lac object at the redshift of \(z=0.34\) [56] and \(\log \nu^{\rm peak}_{\rm syn}=17.9\) (in Hz). The energy of the highest-energy photon is 609.5 GeV, corresponding to an EBL optical depth \(\tau_{\rm EBL} = 2.75\). Our analysis yielded a test statistic (TS) of 433, with a total of 37 photons detected above 100 GeV.

4FGL J2055.4\(-\)​0020 (1RXS J205528.2\(-\)​002123): Identified as a new VHE emitting blazar for the first time in this work, this object is located at \(z=0.44\) [56]. The source has \(\log \nu^{\rm peak}_{\rm syn}=18\) (in Hz) and three photons were detected above 100 GeV at \(>\)​95% association probability, with the highest energy photon of 357 GeV.

5.4 Comparing with other catalogs↩︎

Recently, [57] reported a catalog of VHE emitting AGN candidates at high Galactic latitudes using a photon-counting clustering technique, identifying a total of 175 high-confidence candidates, which include 71 previously known sources, and 127 HSP blazars. In contrast, we employed the full likelihood optimization technique to determine the detection significance and spectral parameters and found 144 VHE emitting HSP blazars. Among these, 106 objects are common in both works, including 42 established IACT-detected VHE-emitting blazars. The remaining 21 objects reported as VHE emitters by [57] have \(9<{\rm TS}<25\) as per our Fermi-LAT data analysis.

Figure 4: (Top left) Dependence of number of VHE photons detected and the variability index. (Top right) Dependence of number of VHE photons and detection significance. (bottom) Variability Index and its dependence on detection significance of the sources in 4FGL Catalog overlaid with FSRQ and Bl lac values and VHE sources identified in this work.

5.5 Source Variability and VHE Time of arrival↩︎

Figure 5: Weekly (black) and Monthly (red) lightcurve of sources highlighted in Figure 4 with a variability index >​70 and atleast 10 VHE photon amongst the known and newly identified VHE sources presented in this work. The time of arrival of VHE photons for each source is highlighted with vertical green dashed lines. The median flux level is marked as a baseline.

a

b

Figure 6: Continued..

Variability provides a crucial diagnostic of the emission processes in blazars, particularly in the VHE domain where photons are expected to originate from compact regions during episodes of enhanced particle acceleration [58]. Since the EGRET era, spectral hardening during flares has highlighted the enhanced detectability of VHE photons in such states. However, the observational strategies of current ground-based VHE facilities — typically triggered by flares — leave it uncertain whether VHE emission is confined to active states or also present at lower, quiescent flux levels. The latter possibility was proposed by [59] to account for the diffuse extragalactic \(\gamma\)-ray background ([60]; [61]), and has been supported by dedicated observations such as the extended low-state campaign of PKS 2155\(–\)​304 from 2005–2007 [62]. Yet, such detections remain limited by sensitivity and the long exposures required.

With Fermi-LAT, it is possible to probe the VHE domain continuously, independent of flaring triggers. In our analysis, we find that the number of detected VHE photons shows a positive correlation with both the source variability index and its overall detection significance, with a Spearman correlation coefficient of 0.48 (Fig. 4). This suggests that sources exhibiting stronger variability tend to yield more VHE photons, consistent with the expectation that highly variable systems undergo more frequent episodes of efficient particle acceleration. At the same time, the variability index and detection significance are themselves tightly correlated, reflecting the fact that brighter sources are more likely to exhibit measurable variability. Disentangling the roles of intrinsic variability and observational bias is therefore critical. We also note that an apparent upper envelope exists in the number of VHE photons at a given variability index.

To investigate whether VHE photon arrival is preferentially associated with active source states, we examined the weekly and monthly binned light curves of sources with the 4FGL-DR4 variability index \(>70\) and at least 10 detected VHE photons (Fig. 5). In these plots, the arrival times of individual VHE photons are marked as vertical green lines.

We find that while several VHE photons coincide with pronounced high-energy flaring episodes, e.g., in 4FGL J1104.4+3812 or Mrk 421, a considerable fraction are instead detected during epochs consistent with the blazar being in a low activity. This indicates that VHE emission is not exclusively tied to strong flares, but can also persist at lower flux states. The presence of such quiescent VHE photons may reflect a steady acceleration component in the jet consistent with earlier claims of baseline VHE activity in blazars [62], [63]. Conversely, the clustering of photon arrivals during flares reinforces the link between efficient acceleration of energetic particles and enhanced VHE production.

These results suggest that both transient flares and underlying quiescent processes contribute to the observed VHE emission, highlighting the need for continuous monitoring to disentangle their relative roles.

5.6 Prospects for followup with Cherenkov Telescopes↩︎

Figure 7: Very-high-energy (VHE) \gamma-ray spectra of significant sources identified in this work, color-coded by flux levels as defined in the inset. Overlaid are the differential sensitivity curves of CTAO, MAGIC, VERITAS, and H.E.S.S. for an observation time of 50 hr.

The newly identified VHE emitting blazars span a flux range covering approximately two orders of magnitude in the 100 GeV\(–\)​2 TeV energy band. To characterize their distribution, we divided the total flux range into three bins:

  1. Low flux: 108 sources with fluxes \(<3 \times 10^{-6}\,\mathrm{MeV\,cm^{-2}\,s^{-1}}\),

  2. Intermediate flux: 30 sources with fluxes in the range \((3\)\(20) \times 10^{-6}\,\mathrm{MeV\,cm^{-2}\,s^{-1}}\),

  3. High flux: 6 sources with fluxes in the range \((2\)\(15) \times 10^{-5}\,\mathrm{MeV\,cm^{-2}\,s^{-1}}\).

In Figure 7, we show the 0.1\(-\)​2 TeV power law spectrum of all sources and overplot the sensitivity limits of several Cherenkov telescopes. As can be seen, blazars present in the intermediate- and high-flux categories should be easily detectable with currently operating Cherenkov facilities. Furthermore, a major fraction of all sources should be detectable with the upcoming CTAO.

6 Summary↩︎

In this work, we analyzed the 0.1\(-\)​2 TeV Fermi-LAT data covering \(\sim\)​16 years of the satellite operation, to identify VHE emitting HSP blazars. We summarize our findings below:

  • We identified 92 new VHE-emitting HSP blazars for the first time, detected at a high confidence level (TS \(\geq\) 25, equivalent to \(\geq\) 5\(\sigma\)). The study also confirmed VHE emission from 52 previously known VHE blazars listed in the TeVCat catalog. Additionally, 187 sources were detected at a lower significance (9 \(\leq\) TS \(<\) 25) and can be considered as promising VHE candidates.

  • A comparison with non-VHE emitting HSP blazars revealed that the identified VHE sources are generally located at lower redshifts (mean \(z = 0.2 \pm 0.1\)), and exhibit brighter synchrotron emission (\(\log F^{\rm peak}_{\rm syn} = -11.22 \pm 0.45\), in erg cm\(^{-2}\) s\(^{-1}\)). These characteristics suggest a potential Malmquist bias, where brighter and closer objects are more easily detected.

  • VHE blazars stand out from the general \(\gamma\)-ray population by their hard spectra and moderate brightness. Interestingly, VHE blazars tend to exhibit relatively low-to-moderate variability indices for their flux level, consistent with expectations from spectral considerations: a hard LAT spectrum (\(\Gamma<2\)) implies an inverse Compton peak at hundreds of GeV, where the LAT band is dominated by lower-energy electrons with longer cooling times, naturally leading to reduced variability compared to the TeV regime.

  • We found that the number of detected VHE photons correlates positively with both the variability index and detection significance, consistent with more variable (and brighter) systems undergoing frequent episodes of efficient particle acceleration. However, this correlation is partly influenced by observational bias, as brighter sources are easier to detect and characterize. By examining the arrival times of VHE photons against source light curves of strongly variable sources, the study shows that while many VHE photons arrive during major flares, a significant fraction are also detected during quiescent states. This dual behavior indicates that VHE emission is not confined to flaring episodes, but can also arise from a steady, baseline acceleration component in jets.

7 Software and third party data repository citations↩︎

In Table [tab:ts95betwen95995and9525], we provide the spectral parameters of HSP blazars detected with \(9<{\rm TS}<25\). The remaining HSP blazars and unassociated \(\gamma\)-ray sources from which at least one VHE photon was detected, are tabulated in Table 1 and [tab:unassociated95sources], respectively.

lcccccc J0013.9\(-\)​1854 & 16 & \(0.71 \pm 0.68\) & \(2.66 \pm 1.37\) & 2 & 314.6 & 58670.1
J0039.1\(-\)​2219 & 21 & \(1.49 \pm 0.96\) & \(2.75 \pm 0.95\) & 2 & 309.1 & 57207.6
J0040.3\(+\)​4050 & 18 & \(1.05 \pm 0.95\) & \(2.42 \pm 1.07\) & 2 & 159.1 & 60008.0
J0045.3\(+\)​2128 & 17 & \(1.10 \pm 1.07\) & \(2.18 \pm 1.03\) & 2 & 246.8 & 54998.7
J0047.9\(+\)​3947 & 13 & \(0.81 \pm 0.70\) & \(2.52 \pm 1.14\) & 2 & 292.2 & 58842.6

Table 1: The remaining HSP blazars from which at least one VHE photon was detected.
Source Name No. E\(_{\rm max}\) Time\(_{\rm arr}\)
\([1]\) \([2]\) \([3]\) \([4]\)
J0018.4+2946 1 127.3 56743.5
J0021.9-5140 1 194.8 60513.5
J0033.3-2040 1 219.7 57500.7
J0035.2+1514 1 264.0 57134.3
J0057.9+6326 1 249.4 54922.0

lcccc J0057.9\(+\)​6326 & 1 & 249.4 & \(54922.0\)
J0357.7\(-\)​6808 & 1 & 110.6 & \(54687.0\)
J0438.0\(-\)​7329 & 1 & 155.1 & \(58435.8\)
J1146.0\(-\)​0638 & 1 & 113.4 & \(56170.9\)
J1452.0\(-\)​4148 & 1 & 346.3 & \(57104.2\)

References↩︎

[1]
Urry, C. M., &Padovani, P. 1995, , 107, 803,.
[2]
Blandford, R., Meier, D., &Readhead, A. 2019, , 57, 467,.
[3]
Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, , 716, 30,.
[4]
Maraschi, L., Ghisellini, G., &Celotti, A. 1992, , 397, L5,.
[5]
Dermer, C. D., Finke, J. D., Krug, H., & Böttcher, M. 2009, The Astrophysical Journal, 692, 32,.
[6]
Sikora, M., Begelman, M. C., &Rees, M. J. 1994, , 421, 153,.
[7]
Mannheim, K. 1993, , 269, 67,.
[8]
Rachen, J. P. 2000, AIP Conf. Proc., 515, 41,.
[9]
Mücke, A., &Protheroe, R. J. 2001, Astroparticle Physics, 15, 121,.
[10]
Giommi, P., Polenta, G., Lähteenmäki, A., et al. 2012, , 541, A160,.
[11]
Ghisellini, G. 1999, Astroparticle Physics, 11, 11,.
[12]
Costamante, L., Ghisellini, G., Giommi, P., et al. 2001, , 371, 512,.
[13]
Foffano, L., Prandini, E., Franceschini, A., &Paiano, S. 2019, , 486, 1741,.
[14]
Paliya, V. S., Domı́nguez, A., Ajello, M., Franckowiak, A., &Hartmann, D. 2019, , 882, L3,.
[15]
Biteau, J., Prandini, E., Costamante, L., et al. 2020, Nature Astronomy, 4, 124,.
[16]
Stecker, F. W., de Jager, O. C., &Salamon, M. H. 1996, , 473, L75,.
[17]
Sambruna, R. M., Maraschi, L., &Urry, C. M. 1996, , 463, 444,.
[18]
Costamante, L., &Ghisellini, G. 2002, , 384, 56,.
[19]
Colafrancesco, S., & Giommi, P. 2006, Chinese Journal of Astronomy and Astrophysics, 6, 47,.
[20]
Weekes, T. C. 2008, in American Institute of Physics Conference Series, Vol. 1085, American Institute of Physics Conference Series, ed. F. A. Aharonian, W. Hofmann, & F. Rieger(AIP), 3–17,.
[21]
Arsioli, B., &Chang, Y. L. 2017, , 598, A134,.
[22]
Wakely, S. P., &Horan, D. 2008, in International Cosmic Ray Conference, Vol. 3, International Cosmic Ray Conference, 1341–1344.
[23]
Böttcher, M. 2007, , 309, 95,.
[24]
Saugé, L., & Henri, G. 2004, The Astrophysical Journal, 616, 136,.
[25]
Murase, K., Inoue, Y., &Dermer, C. D. 2014, , 90, 023007,.
[26]
Suray, A., & Troitsky, S. 2023, Monthly Notices of the Royal Astronomical Society: Letters, 527, L26,.
[27]
Biteau, J., & Williams, D. A. 2015, The Astrophysical Journal, 812, 60,.
[28]
Giommi, P., & Padovani, P. 2015, Monthly Notices of the Royal Astronomical Society, 450, 2404,.
[29]
Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, , 440, 1018,.
[30]
Neronov, A., & Vovk, I. 2010, Science, 328, 73,.
[31]
Finke, J. D., Reyes, L. C., Georganopoulos, M., et al. 2015, The Astrophysical Journal, 814, 20,.
[32]
Aharonian, F., Aschersleben, J., Backes, M., et al. 2023, The Astrophysical Journal Letters, 950, L16,.
[33]
H. E. S. S. Collaboration, Aharonian, F., Ait Benkhali, F., et al. 2025, , 695, A261,.
[34]
Tomar, G., Paliya, V. S., Saikia, D. J., &Stalin, C. S. 2026, Journal of High Energy Astrophysics, 49, 100454,.
[35]
Paliya, V. S., Böttcher, M., Wani, K., et al. 2025, , 991, L8,.
[36]
Abdollahi, S., Acero, F., Baldini, L., et al. 2022, , 260, 53,.
[37]
Ballet, J., Bruel, P., Burnett, T. H., Lott, B., &The Fermi-LAT collaboration. 2023, arXiv e-prints, arXiv:2307.12546,.
[38]
Arsioli, B., Fraga, B., Giommi, P., Padovani, P., &Marrese, P. M. 2015, , 579, A34,.
[39]
Chang, Y. L., Arsioli, B., Giommi, P., &Padovani, P. 2017, , 598, A17,.
[40]
Chang, Y. L., Arsioli, B., Giommi, P., Padovani, P., &Brandt, C. H. 2019, , 632, A77,.
[41]
Wood, M., Caputo, R., Charles, E., et al. 2017, PoS, ICRC2017, 824,.
[42]
Ajello, M., Atwood, W. B., Axelsson, M., et al. 2021, , 256, 12,.
[43]
Acero, F., Ackermann, M., Ajello, M., et al. 2016, , 223, 26,.
[44]
Domı́nguez, A., Primack, J. R., Rosario, D. J., et al. 2011, , 410, 2556,.
[45]
Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2008, , 481, L103,.
[46]
Tanaka, Y. T., Cheung, C. C., Inoue, Y., et al. 2013, , 777, L18,.
[47]
Hauser, M. G., &Dwek, E. 2001, , 39, 249,.
[48]
Arsioli, B., Chang, Y. L., &Musiimenta, B. 2020, , 493, 2438,.
[49]
Paliya, V. S., Domı́nguez, A., Ajello, M., Olmo-Garcı́a, A., &Hartmann, D. 2021, , 253, 46,.
[50]
Ajello, M., Baldini, L., Ballet, J., et al. 2022, , 263, 24,.
[51]
Chiaberge, M., &Ghisellini, G. 1999, , 306, 551,.
[52]
Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2020, , 248, 29,.
[53]
Saldana-Lopez, A., Domı́nguez, A., Pérez-González, P. G., et al. 2021, Mon. Not. Roy. Astron. Soc., 507, 5144,.
[54]
Paiano, S., Landoni, M., Falomo, R., et al. 2017, , 837, 144,.
[55]
Koss, M. J., Ricci, C., Trakhtenbrot, B., et al. 2022, , 261, 2,.
[56]
Shaw, M. S., Romani, R. W., Cotter, G., et al. 2013, , 764, 135,.
[57]
Neronov, A., &Semikoz, D. 2025, arXiv e-prints, arXiv:2506.08497,.
[58]
—. 2007, , 664, L71,.
[59]
Stecker, F. W., &Salamon, M. H. 1996, , 464, 600,.
[60]
Fichtel, C. 1996, , 120, 23.
[61]
Sreekumar, P., Bertsch, D. L., Dingus, B. L., et al. 1998, The Astrophysical Journal, 494, 523,.
[62]
H. E. S. S. Collaboration, Abramowski, A., Acero, F., et al. 2010, , 520, A83,.
[63]
Dmytriiev, A., Sol, H., &Zech, A. 2021, , 505, 2712,.

  1. https://www.tevcat.org/↩︎