Very High-energy (>50 GeV) Gamma-Ray Flux Variability of Bright Fermi Blazars

Understanding the high-energy emission processes and variability patterns are two of the most challenging research problems associated with relativistic jets. In particular, the long-term (months to years) flux variability at very high energies (VHE >50 GeV) has remained an unexplored domain so far. This is possibly due to the decreased sensitivity of the Fermi Large Area Telescope (LAT) above a few GeV, hence low photon statistics, and observing constraints associated with the ground-based Cherenkov telescopes. This paper reports the results obtained from the 0.05−2 TeV Fermi-LAT data analysis of a sample of 29 blazars with the primary objective to explore their months-to-year-long very high-energy (VHE) flux variability behavior. This systematic search has led to, for the first time, the detection of significant flux variations in five blazars at the >99% confidence level, whereas eight of them exhibit variability, albeit at a lower confidence level (∼95%–99%). A comparison of the 0.05–2 TeV flux variations with that observed at 0.1–50 GeV band has revealed similar variability behavior for most of the sources. However, complex variability patterns that are not reflected contemporaneously in both energy bands were also detected, thereby providing tantalizing clues about the underlying radiative mechanisms. These results open up a new dimension to unravel the VHE emission processes operating in relativistic jets, hence sowing the seeds for their future observations with the upcoming Cherenkov Telescope Array.


INTRODUCTION
The high-energy emission and the erratic, rapid, and large amplitude variability observed in all accessible spectral regimes (radio-to-γ-ray) are two of the main defining properties of blazars (e.g., Marscher 2016), a class of radio-loud active galactic nuclei (AGN) that have powerful relativistic jets oriented close to the line of sight to the observer.They are the most variable (on years-to-minute timescales), the most luminous (L bol.> 10 46 erg s −1 ), the most polarized (even > 30% in the optical and radio bands), and thus probably the most exciting type of AGN.The emitted radiation from the jet is strongly amplified due to relativistic effects, so-called Doppler boosting.This flux enhancement makes blazars one of the very few types of astrophysical objects detectable throughout the accessible electromagnetic spectrum and the most abundant sources in the extragalactic γ-ray sky Corresponding author: Vaidehi S. Paliya vaidehi.s.paliya@gmail.com(e.g., Ajello et al. 2022).Blazars are categorized into flat spectrum radio quasars (FSRQs) and BL Lac objects based on the rest-frame equivalent width (EW) of their broad optical emission lines, with BL Lac sources having EW<5 Å (cf.Stickel et al. 1991).
The spectral energy distribution (SED) of a blazar is dominated by non-thermal jet emission and characterized by two broad humps.The low energy peak, produced by the synchrotron process, lies in the sub-millimeter (mm) to the UV/X-ray energy regime, whereas, the high energy hump from inverse Compton mechanism peaks in the MeV/GeV range, assuming leptons to be responsible for the observed radiation.Alternatively, hadronic processes (π 0 decay and subsequent cascade radiation) have also been proposed to explain the high-energy emission from blazars (e.g., Böttcher et al. 2013).Interestingly, hadronic models predict blazar jets as promising candidates for neutrino emission (cf.Mannheim & Biermann 1989;Reimer et al. 2019).This indicates that relativistic jets can be among the most efficient particle accelerators.Furthermore, based on the location of the syn-chrotron peak frequency (ν peak syn ), blazars are also classified as low-(ν peak syn <10 14 Hz), intermediate-(10 14 < ν peak syn <10 15 Hz) and high-(ν peak syn >10 15 Hz) synchrotron peaked sources (Abdo et al. 2010b).
Flux variability is probably the most studied topic in jet physics, yet we do not have a clear consensus about the underlying physical processes, the behavior of the variability pattern across various blazar sub-classes and connection with the accretion and central black hole, the variability timescales, and multi-band flux correlations.There have been several studies done to understand the blazar flux variability across the electromagnetic spectrum, e.g., at radio (cf.Fuhrmann et al. 2016), infrared (e.g., Mao et al. 2018;Anjum et al. 2020), optical-ultraviolet (cf. Edelson 1992;Stalin et al. 2004;Bonning et al. 2012;Paliya et al. 2017), X-ray (Pian 2002;Rani et al. 2017), and γ-ray energies (cf.Abdo et al. 2010a).Additionally, multi-wavelength studies focused both on short-(≲days) and long-term (months-to-years) flux variations have been carried out to understand the radiative processes, jet composition, and kinematics (e.g., Böttcher & Dermer 1998;Sokolov & Marscher 2005;Nalewajko et al. 2011;Marscher 2014;Tavani et al. 2015;Böttcher & Baring 2019).
At very high energies (VHE, >50 GeV), blazar observations are often triggered by elevated activity states identified at lower frequencies (e.g., Adams et al. 2022).These followup observations with Cherenkov telescopes usually cover the source activity for days to weeks, thus providing an opportunity to study the event on short timescales.The long-term, i.e., months-to-years, VHE flux variations of blazars, on the other hand, remained almost an unexplored territory with studies have been limited so for to the brightest sources (e.g., Acciari et al. 2014;Abeysekara et al. 2017;Valverde et al. 2020).This is possibly due to telescope time availability, weather constraints, and/or source visibility from the groundbased facilities.This work attempts to address this outstanding research problem using VHE observations carried out with the Fermi-Large Area Telescope (Fermi-LAT, Atwood et al. 2009).Due to the surveying capabilities of the Fermi-LAT, almost uninterrupted monitoring of the whole γ-ray sky is possible at energies as high as ≳2 TeV.Though the sensitivity of the Fermi-LAT starts decreasing above a few GeV energies, so the background contamination.This enables an improved source localization and, thus, the more accurate association of the γ-ray photons with a point source.Admittedly, the temporal VHE flux variations cannot be studied using the Fermi-LAT on shorter timescales due to the faintness of blazars at these energies, hence low photon statistics.However, integrating the LAT data on longer timescales allows one to study the months-to-year flux variations.This is the primary objective of this work.In Section 2, and 3, the details of the sample selection and data reduction steps are described.The derived results are presented in Section 4, and noteworthy points about individual blazars are highlighted in Section 5.The overall findings are summarized in Section 6.

SAMPLE SELECTION
The third data release of the fourth catalog of Fermi-LAT detected γ-ray sources (4FGL-DR3, Abdollahi et al. 2022) includes 172 objects with detection significance ≳4.2σ (test statistic TS >25, Mattox et al. 1996) in the energy range of 100 GeV-to-1 TeV.Considering only blazars or blazar candidates of uncertain type (BCU), the sample size was reduced to 125 sources.The Fermi-LAT data reduction covering the first ∼15 years of the mission was carried out for all of them using the methodology described in Section 3.Only those 29 objects were retained whose overall detection significance was estimated to be ≳10σ (TS >100).This criterion ensures that a meaningful temporal study can be done on the selected γ-ray sources.This final sample contains 28 BL Lac objects and one BCU listed in Table 1.All sources have been flagged as variable in the 4FGL-DR3 catalog (see also Abdollahi et al. 2023).

FERMI-LAT DATA REDUCTION
The standard data reduction steps were followed to analyze the ∼15 years of the Fermi-LAT data (modified Julian date or MJD 54683−60220, i.e., 2008 August 5 to 2023 October 2).In the energy range of 0.1−2 TeV, only SOURCE class events were selected from a circular region of radius 5 • which passed the zenith angle-based cut (z max < 105 • ).The Pass 8 data includes division of photons by pointspread function (PSF) type with the worst (PSF0) and best (PSF3) best angular reconstructions, respectively.Therefore, a component-wise LAT data analysis for each PSF was performed.The final fitting was carried out using the SUMMED likelihood method of the FermiTools.For the likelihood fitting, all γ-ray sources present in the 4FGL-DR3 catalog and lying within 7 • from the position of interest were considered to represent the γ-ray sky.The model file also includes latest diffuse background models1 (e.g., Acero et al. 2016).Due to low photon statistics, the γ-ray spectrum of the source of interest was assumed to follow a power law spectral shape.The first round of the model optimization was performed using fermiPy (Wood et al. 2017) to estimate the detection significance, i.e., TS, of all sources and the spectral parameters of those with TS >25 were left free to vary during the likelihood fitting.TS maps were also generated to search for any unmodeled γ-ray objects present in the data but not in the model since the time period covered in the analysis is larger than that used to produce the 4FGL-DR3 catalog.Once identified, if any, they were included in the model assuming a Table 1.The spectral information of blazars in 0.1−2 TeV energy range derived from the analysis of ∼15 years of the Fermi-LAT observations.Column information are as follows: (1) 4FGL name; (2) association; (3) redshift; (4) (5), and (6) observed photon flux (in units of 10 −11 ph cm −2 s −1 ), photon index, and test statistic, in 0.1−2 TeV energy range; (7) (8), and (9) EBL-corrected photon flux (in units of 10 −11 ph cm −2 s −1 ), photon index, and test statistic, in 0.1−2 TeV energy range; (10) whether the source was detected with ground-based Cherenkov telescopes, Y: yes, N: no; (11) variability test statistic; and (12) status whether a source is variable (V), probably variable (PV), or non-variable (NV) in 0.1−2 TeV energy range.

4FGL name
Other Since the extragalactic background light (EBL) attenuation of γ-ray photons becomes important at the considered energy range, the spectral fitting was also performed including the EBL absorption to determine the EBL-corrected spectral parameters.The EBL model proposed by Domínguez et al. (2011) was used for this purpose.
After the ∼15 years averaged data analysis, 29 objects were selected whose overall TS was estimated to be >100.The 0.1−2 TeV spectral parameters of these sources are provided in Table 1.To generate the γ-ray light curves, 3 months time binning was adopted and the minimum energy considered in the analysis was changed to 50 GeV.This was done to improve the photon statistics while ensuring that the analysis probes the VHE domain which is the main objective of the work.The EBL-corrected spectral parameters, i.e., normalization and photon index, of the source of interest were allowed to vary during the likelihood fitting carried out in the individual time bin.Whenever the TS of the target object was found to be <4, the flux upper limit at 95% confidence level was derived.This choice is slightly different from that conventionally adopted while analyzing 0.1−300 GeV data, i.e., TS< 9.However, all blazars under consideration are welldetected with the Fermi-LAT and several of them have also been detected with the ground-based Cherenkov telescopes (Table 1).Therefore, a slightly reduced detection significance can be adopted to maximize the number of data points in the light curves.
To compare the flux variations observed in 0.05−2 TeV energy range (hereafter TeV band) with that at lower γ-ray energies, the light curves were also generated in 0.1−50 GeV energy band (hereafter GeV band) using the same time binning.An ROI of 15 • radius was chosen, and all 4FGL-DR3 sources lying within 20 • from the position of interest were considered for the likelihood fitting.The spectral shapes were directly taken from the 4FGL-DR3 catalog to model the 0.1−50 GeV energy spectrum.

RESULTS
The light curves for both TeV and GeV bands are shown in Figure 1.In the low-energy GeV band, most of the sources exhibited strong flux variability, however, only a small fraction appeared variable at very high frequencies.The extent of temporal flux variations was quantified by estimating the variability test statistic (TS var ) following a maximumlikelihood approach (cf.Nolan et al. 2012;Abeysekara et al. 2017): where L i (F i ) is the likelihood with F i being the best-fit flux value for time bin i.The parameter L i (F avg ), is the likelihood with F avg being is the best-fit flux value for the full time period derived from the likelihood fitting.The parameter TS var is distributed as χ2 with 61 degrees of freedom if the null hypothesis is correct, and a value of TS var > 80.2 or > 89.6 corresponds to a source being variable at >95% or >99% confidence level, respectively.The TS var values estimated for all sources are reported in Table 1.A blazar was considered variable (V) for >99% confidence level and probably variable (PV) if it had a value between 95% and 99%.The remaining sources were considered non-variable (NV).This exercise led to the identification of 5 objects showing significant flux variability in the TeV band, including Mkn 421, which is the most variable source in the sample.There are 8 blazars whose TeV light curves hint for the existence of flux variability, albeit at a lower confidence level.Remaining 16 objects did not show any significant variability.Even though several blazars do not show strong VHE flux variations, likely due to low photon statistics, interesting variability patterns emerge when comparing the light curves in the TeV and GeV energy bands.To quantify the observed patterns, a time-lag analysis was carried out using the z-transformed discrete correlation function (z-DCF) technique and the uncertainties were estimated using a likelihood method (Alexander 1997(Alexander , 2013)).The z-DCF method provides a conservative estimate for the cross-correlation as a function of time-lag between unevenly sampled light curves.Furthermore, blazars often exhibit erratic flux variations with the same object some time showing correlated flux changes and uncorrelated variability during other epochs (e.g., 3C 279, Hayashida et al. 2015;Paliya et al. 2015Paliya et al. , 2016)).Additionally, only flux upper limits could be derived during several epochs for various objects under consideration and that may somewhat bias the result as those upper limits would cover mainly the low states of the source.Therefore, rather than running on the full dataset, z-DCF analysis was carried out separately on specific epochs.The time periods were selected in such a manner to cover the interesting activity epochs while maximizing the number of data points.The results are discussed in the next section.
The Bayesian block representation of the time series data allows to identify the periods of different activity states.An advantage of this technique is that it is able to identify significant changes in light curve independently of variations in exposure or gaps (Scargle et al. 2013;Wagner et al. 2022).Therefore, the best-fit piece wise constant representations of the Fermi-LAT data were derived using this method.The false-alarm probability was chosen as 25% and 5% for the TeV band and all other (GeV band and hard X-ray) light curves, respectively.For 22 out of 29 objects, this exercise resulted in the identification of blocks of different TeV band flux states.The results are shown in Figure 1.
In the conventional leptonic radiative models, the VHE radiation observed from high-synchrotron peaked BL Lac objects is produced by the highest energy electrons via inverse Compton scattering of synchrotron photons (e.g., Tavecchio et al. 2010).At lower frequencies, the same electron population is expected to radiate synchrotron emission usually lying in the hard X-ray band (>10 keV, Costamante et al. 2018;Foffano et al. 2019).Therefore, it is imperative to compare the TeV band flux variations with that observed at hard X-rays.For this purpose, Crab weighted, monthly binned light curves of 11 Swift-Burst Alert Telescope (BAT) detected objects were retrieved from the 157-month Swift-BAT hard X-ray survey website 2 .The Swift-BAT data overlaps with the Fermi-LAT observations in the time period MJD 54683−58103, i.e., up to 2017 December 16.For the analysis, only those hard X-ray data points were considered that have the signal-to-noise ratio >1.The light curves for both TeV and hard X-ray bands are shown in Figure 2.    1ES 0033+595: The TeV band light curve of this Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescopes detected source (Aleksić et al. 2015a) revealed a flaring activity in the time period MJD 56392−57562.Interestingly, the pattern of the GeV flux variations suggests the flare to peak at a later time than the TeV flare.However, the results of the z-DCF analysis indicated the absence of a significant lead/lag due to large uncertainties (Table 2).The Bayesian block analysis suggested a high flux state overlapping in both TeV and GeV bands.The source is overall nonvariable in the TeV band as per the TS var calculation.Furthermore, a comparison of the TeV and 14−195 keV light curves revealed the source to be in an elevated hard X-ray activity state at the time of the TeV flare (Figure 2).

NOTES ON INDIVIDUAL OBJECTS
3C 66A: This well-known Very Energetic Radiation Imaging Telescope Array System (VERITAS) detected blazar was in a high GeV activity state during the first year of the Fermi-LAT operation (Abdo et al. 2011b).The high activity was also reflected in the TeV band light curve (Figure 1; Aleksić et al. 2011).The z-DCF analysis suggested the contemporaneous flux enhancements in both TeV and GeV bands (Table 2).Later, the object went into a low activity state but was frequently detected in the TeV band at a flux level comparable to its ∼15 years averaged value.The source can be considered moderately variable based on the estimated TS var value.The Bayesian block analysis carried out on the TeV band light curve hinted for the presence of variable flux activity states.
PKS 0447−439: The TeV emission from this highsynchrotron peaked blazar was first detected with the High Energy Stereoscopic System (H.E.S.S.) during 2009-2010 (H.E. S. S. Collaboration et al. 2013) consistent with the low-activity detection in the TeV band light curve (Figure 1).The source has shown moderate variability in the TeV band.It exhibited its brightest TeV flare in the time bin MJD 59002−59092.This activity was also reflected in the GeV band light curve.Interestingly, a GeV flare of similar amplitude was observed from this object earlier in the time bin MJD 57652−57742, however, the source was detected at its ∼15-years averaged flux level in the TeV band.To ascer- tain this, the z-DCF analysis was employed on two time periods: MJD 56752−57922 and MJD 58282−59272.The z-DCF coefficient for the first period was found to be 0.29 +0.42   −0.50   indicating the two events to be only marginally correlated.
For the second period, TeV and GeV flux variations were strongly correlated (z-DCF coefficient = 0.89 +0.09 −0.14 ) with no significant lead/lag.Moreover, the Bayesian block analysis revealed the existence of variable flux activity states and the source is found as mildly variable based on the estimated TS var value.
1ES 0502+675: A ∼2.5-year long TeV flaring event was observed by the Fermi-LAT during the initial years of its operation, and the source was also detected with VERITAS in this time-period (Figure 1; Ong 2009).Due to large uncertainties in both TeV and GeV bands data, however, the results of the z-DCF analysis were inconclusive.The source went into quiescence after this flare and was sporadically detected in the TeV band at a flux level comparable to its ∼15 years averaged value as also confirmed by the Bayesian block analysis.Overall, significant TeV flux variability at >99% confidence level was detected from this blazar.The source appeared to be bright in hard X-rays during the elevated TeV activity state which is confirmed by the z-DCF analysis.There are hints for the X-ray flare leading the TeV event, however, it cannot be claimed due sparseness of the data points and large uncertainties in the derived parameters (Table 2).
TXS 0518+211: This high-synchrotron peaked BL Lac object was detected with VERITAS during MJD 55126−55212 (Archambault et al. 2013).A flaring activity in the TeV band was identified in the time period MJD 56302−56752 which was reflected in the GeV band and also followed by the Cherenkov telescopes (Adams et al. 2022).Later, the source was detected in the TeV band on a few occasions and found to exhibit moderate flux variations as quantified by the TS var analysis.
RX J0543.9−5532: Brown et al. (2014) reported the first VHE detection of this source using the Fermi-LAT data above 100 GeV though it was not detected with the Cherenkov telescopes (H.E. S. S. Collaboration et al. 2014).The GeV and TeV band light curves indicated that the object was in quiescence most of the time since the Fermi-LAT operation began in 2008.According to the Bayesian block analysis, the enhanced GeV band activity was identified starting with the time bin MJD 58642−58732 and the source was also detected in the TeV band.However, no significant TeV flux variability was identified as revealed by the TS var analysis.
RX J0648.7+1516:This blazar was detected only on a few occasions in the TeV band during the ∼15 years of the Fermi-LAT operation and can be considered non-variable from the estimated TS var .The GeV light curve of this object also does not exhibit any significant flux variability.
1ES 0647+250: This object is one of the few blazars detected at TeV energies by MAGIC telescopes even during non-flaring activity states (Figure 1; MAGIC Collaboration et al. 2023).No significant flux variations were noticed in the TeV band light curve which was also confirmed by the derived TS var value.
S5 0716+71: This γ-ray bright blazar has exhibited exceptional GeV flaring activities and is also detected at TeV energies by the MAGIC telescopes (Anderhub et al. 2009;MAGIC Collaboration et al. 2018;Geng et al. 2020).The large amplitude flux variability can also be seen in the GeV band light curve, whereas, the flux variations observed in the TeV band data were minor.The largest TeV band flare was observed in the time period MJD 56032−56212 which was also reflected at lower energies.The elevated TeV flux activity was also detected by Cherenkov telescopes on several occasions which are consistent with the TeV band detections (Figure 1; Mirzoyan 2017;MAGIC Collaboration et al. 2018).Overall, the estimated TS var has indicated that the source is partially variable.This object is also detected with Swift-BAT and a comparison of the TeV and hard X-ray light curves revealed the lack of hard X-ray counterpart during the TeV flare (MJD 56032−56212).Though it was not possible to quantify the lack of correlation due to sparseness of the data points, this observation can be understood from the SED behavior of the source.It is an intermediate synchrotron peaked BL Lac object whose 'valley' between the two broad humps lies in the X-ray band (cf.MAGIC Collaboration et al. 2018).In such a scenario, the hard X-ray (>10 keV) radiation is primarily dominated by the inverse Compton emission produced by the low-energy electron population, thus, less variable than the VHE emission produced by the highest energy electrons.
1ES 0806+524: This object was detected in the TeV band on several occasions at a flux level consistent with its ∼15 years averaged value.The TS var analysis has revealed it to be non-variable and the Bayesian block analysis has also not found any significant variable flux states.The GeV band light curve highlighted a flaring activity in the time period MJD 55582−55852 in which the TeV band flux shows an enhancement as well.The source was also detected with MAGIC telescopes during this flaring episode (Figure 1; Aleksić et al. 2015b).
1H 1013+498: This blazar was detected with the MAGIC telescopes coincident with the observation of an optical outburst (Albert et al. 2007).From the Fermi-LAT data analysis, only a moderate level of flux variability was noticed in both the GeV and TeV band light curves, also confirmed by the Bayesian block routine.The source is non-variable according to the estimated TS var value.A comparison of the TeV and GeV bands light curves revealed that the elevated flux activities noticed in the former were also reflected in the latter.Furthermore, the source was detected with Cherenkov telescopes during these epochs (Ahnen et al. 2016;Aleksić et al. 2016).
GB6 J1037+5711: The TeV light curve of this object does not show any significant flux variability as confirmed by the TS var value.It was detected only on a few occasions in the TeV band and the Bayesian block analysis did not reveal any significant variable flux states.The High Altitude Water Cherenkov Observatory reported the detection of a TeV transient event from this blazar during MJD 58298−58300 (Weisgarber & HAWC Collaboration 2018).This epoch lies in the time bin MJD 58282−58372, in which GB6 J1037+5711 was significantly detected in the TeV band (Figure 1).
Mkn 421: This object is probably the brightest and the most studied TeV blazar to date (cf.Abdo et al. 2011a).It was detected in the TeV band throughout the ∼15 years of the Fermi-LAT operation and also by several Cherenkov telescopes (Figure 1; Bartoli et al. 2016;Garcia-Gonzalez & Martinez 2019;Abeysekara et al. 2020;Arbet-Engels et al. 2021).The pattern of the flux variability in both the GeV and TeV bands appears similar, with the high/low TeV activity states coinciding with that observed in the GeV band.However, the object also exhibited complex variability behavior on a few occasions.For example, the high activity noticed in the TeV band in the time bin MJD 56122−56212 was also seen in the GeV band, and the estimated GeV flux was found to be the highest observed from the source during the first ∼15 years of the Fermi-LAT operation.On the other hand, the GeV flux corresponding to a TeV flare of similar amplitude, observed in the time bin MJD 56482−56572, was considerably lower than that observed during the previous TeV flare described above.The z-DCF analysis also indicates the existence of a mild correlation (Table 2).The observation of such uncorrelated or weakly correlated GeV−TeV flux variations from Mkn 421 are not uncommon (e.g., Bartoli et al. 2016;Acciari et al. 2020).Interestingly, during the both TeV flaring events, Mkn 421 was too close to the Sun to be followed up by the multi-wavelength observatories implying that two of the brightest TeV outbursts detected by the Fermi-LAT could not be followed up.This finding also highlights the importance of an all-sky surveying γ-ray mission to capture the extraordinary flaring episodes, some of which, otherwise, may remain unexplored due to practical constraints for other facilities.Interestingly, a comparison of the TeV and hard X-ray light curves reveal the source to be in quiescence in the 14−195 keV band at the time of the TeV flare observed in the time bin MJD 56122−56212 (Figure 2).This uncorrelated variability pattern was also confirmed with the z-DCF analysis (Table 2).This is probably a surprising result since the VHE and X-ray observations are often strongly correlated (e.g., González et al. 2019).
RBS 0970 and Mkn 180: The TeV band light curves of these two blazars do not exhibit any significant flux variability and detected only on a few occasions.Furthermore, the GeV flux variations were modest, as seen in their light curves.The TS var analysis revealed both of them to be nonvariable.
1ES 1215+303: The long-term VHE flux variability behavior of this source was recently studied by Valverde et al. (2020) using VERITAS observations and significant flux variations were detected (see also Aleksić et al. 2012a).The TeV band light curve has reveled several γ-ray flaring episodes which were also found by the Bayesian block algorithm.The high activity was reflected in the GeV band light curve.Based on the estimated TS var parameter, this object can be considered partially variable.
PG 1218+304: The TeV band light curve of this highsynchrotron peaked blazar showed a low-amplitude γ-ray flare in the first year of the Fermi-LAT operation which coincided with the detection by VERITAS and observation of day-scale flux variability (Figure 1; Acciari et al. 2010a).The source was detected several times displaying a moderate level of TeV flux variability which is confirmed by the derived TS var value (Table 1).The Bayesian block analysis also revealed the presence of variable flux states.A compari-son of the TeV and hard X-ray light curves indicated that the low-level flux variability is present in both bands.
PKS 1424+240: At the redshift of z = 0.60 (Paiano et al. 2017), this is the most distant blazar among the sources considered in this work.This blazar was significantly detected in several time bins of the TeV band light curve and also by the VERITAS and MAGIC telescopes (Acciari et al. 2010b;Aleksić et al. 2014b;Archambault et al. 2014).The TeV flux variations were accompanied by that observed in the GeV band light curve which was confirmed by the overlapping Bayesian blocks.The results of the z-DCF analysis carried out during the high-activity period (MJD 55357−56707) also supported this finding though a strong claim cannot be made due to the large uncertainties in the z-DCF coefficient (Table 2).The estimated TS var parameter remained below the variability threshold thereby rendering this object as nonvariable.
PKS 1440−389: This blazar has not exhibited any significant flux variability in the TeV band.It was sporadically detected throughout the ∼15 years of the Fermi-LAT operation and and flux enhancement was seen mainly in the later half of the covered time period.The GeV band flux showed a similar behavior and the brightest GeV flux identified in the time bin MJD 59002−59092 had a TeV band counterpart.The Bayesian block analysis has found variable flux states.This source was detected with the H.E.S.S. during a low activity state (Abdalla et al. 2020).
PG 1553+113: This bright VHE emitting blazar has exhibited quasi-periodic variability in the γ-ray band (Ackermann et al. 2015).It was detected almost throughout all the analyzed periods and has been detected with several Cherenkov telescopes (Figure 1; Aleksić et al. 2012b;Aliu et al. 2015;Abramowski et al. 2015;H. E. S. S. Collaboration et al. 2017b;Gueta 2019).Based on the TS var value, this object is partially variable in the TeV band.Several episodes of moderate amplitude, i.e., a factor of ∼2 with respect to the ∼15 years averaged flux, flaring activity have been noticed which were also reflected in the GeV band light curve.In the time period MJD 58462−59542, two similar amplitude TeV flaring events peaking in the time bins MJD 58642−58732 and MJD 59272−59362 were identified.Though the GeV band light curve showed flux enhancement at both epochs, the amplitude of flux variations were different with the second flare having higher GeV flux which is also confirmed by the Bayesian block analysis.The source was detected with MAGIC telescopes during this flare (Blanch 2021).The z-DCF analysis hinted an overall positive correlation though the uncertainties in the estimated z-DCF coefficient are large (Table 2).
Mkn 501: This well-studied, bright TeV blazar has exhibited violent flux variability across the electromagnetic spectrum (e.g., Abdalla et al. 2019).The TeV band light curve of this source has revealed several episodes of VHE outbursts, with some of them also being detected with various Cherenkov telescopes (cf.Ahnen et al. 2018;Abdalla et al. 2019;Abe et al. 2023).The estimated TS var parameter suggests the source to be variable in the TeV band.Furthermore, it remained mostly in quiescence during the later half of the Fermi-LAT observation period.During the periods of high TeV activity, e.g., in the time interval MJD 55762−57562, the GeV band flux showed similar variations within uncertainties, which was also confirmed by the z-DCF analysis (Table 2).This source is detected with the Swift-BAT and a comparison of the TeV band and hard X-ray light curves indicated contemporaneous flux variations in both energy bands which was supported by the z-DCF analysis.
I Zw 187: This blazar was significantly detected only on a few occasions in the TeV band light curve and has also been detected with Cherenkov telescopes (Figure 1; Aleksić et al. 2014a;Archambault et al. 2015).A high TeV band activity was noticed in the time period MJD 56572−57112 and was found to have coincided with an elevated GeV flaring episode as also confirmed by the Bayesian block analysis.The source is non-variable in the TeV band as per the computed TS var value (Table 1).
1H 1914−194: This object was detected in the TeV band only on a few occasions, and the TS var parameter revealed it to be non-variable.The GeV band light curve shows only a moderate level of flux variability.The Bayesian block analysis did not find any significant variable flux states.
1ES 1959+650: This source is a bright VHE emitting blazer (e.g., Aliu et al. 2013), and the detection of an orphan TeV flare was also reported (Krawczynski et al. 2004).This is one of the few sources studied in this work that were consistently detected in the TeV band throughout the ∼15 years of the Fermi-LAT operation.The TeV band light curve of this blazar has indicated significant flux variability, which was also accompanied by the flux variations seen in the GeV band light curve.In particular, a ∼4 years long (∼MJD 57022−58462) γ-ray flare was identified in both the TeV and GeV band light curves which was confirmed by the Bayesian block analysis.The multiwavelength observations close in time to the peak of this γ-ray flare have been reported by MAGIC Collaboration et al. (2020b).In the time period covered by the Fermi-LAT and Swift-BAT, similar flux variations were observed in TeV and hard X-ray bands light curves.
PKS 2005−489: This blazar was in an elevated activity state during the first year of the Fermi-LAT operation as revealed by its TeV, GeV, and hard X-ray light curves.The source was also detected with the H.E.S.S. observatory in this time period (H.E. S. S. Collaboration et al. 2011, Figure 1;[).After that, the object went into quiescence for the rest of the observing period, exhibiting no significant variability and was sporadically detected at TeV energies.Over-all, it is non-variable in the TeV band based on the estimated TS var value.
RGB J2056+496: This object was occasionally detected in the TeV band at a flux level consistent with the ∼15 years averaged value.The calculated TS var parameter suggested this blazar to be non-variable.This object was also detected with the VERITAS (Mukherjee & VERITAS Collaboration 2016).The hard X-ray band light curve of this object showed an enhanced flux state around MJD 57524 (Figure 2), however, the presence/absence of the TeV counterpart could not be quantified due to sparseness of the data points.
PKS 2155−304: This bright VHE-detected blazar has exhibited significant flux variability at GeV energies and has been studied from several Cherenkov telescopes (Figure 1 H. E. S. S. Collaboration et al. 2017a,b;Wierzcholska et al. 2019).The detection of the rapid, ∼minute-scale, TeV flux variability was also reported (Aharonian et al. 2007).In the TeV band, it was found to be partially variable based on the derived TS var value.As seen in its TeV band light curve, it was almost regularly detected throughout the time period considered in this work.The elevated TeV activity states were also reflected in the GeV band.
BL Lacertae: This eponymous blazar has exhibited several γ-ray flares since the Fermi-LAT observations began in 2008, with the most prominent one observed during 2020-2021 (cf.Pandey & Stalin 2022).In the TeV band, the source remained in quiescence for most of the observing run prior to the 2020-2021 flaring episode though a rapid TeV flare was detected in 2015 June (MAGIC Collaboration et al. 2019).It was regularly detected in the TeV band during MJD 58822−60082, i.e., the period of the high activity, and exhibited huge TeV flares, consistent with its detection with MAGIC telescopes in this period(e.g., Blanch 2020;Blanch et al. 2021).The GeV band light curve also revealed large amplitude flares contemporaneous with the elevated TeV activity which was confirmed with the z-DCF analysis (Table 2).During the time period common to Fermi-LAT and Swift-BAT, the source remained in quiescence in the hard Xray and TeV bands.
1ES 2344+514: This source was sporadically detected in the TeV band light curve and no significant flux variability was noticed (Table 1).In the GeV band also, it has exhibited a moderate level of flux variations.This source has been detected several times with various Cherenkov telescopes (Aleksić et al. 2013;Allen et al. 2017;MAGIC Collaboration et al. 2020a).A comparison of the TeV and hard X-ray bands light curves indicated no significant flux variability.

SUMMARY
This work presents the first results of the 0.05−2 TeV LAT data analysis to explore the long-term VHE flux variations of bright Fermi-LAT detected blazars.Among the sample of 29 objects, 5 sources were found to exhibit statistically significant TeV band flux variability, whereas, 8 of them have shown moderate levels of VHE flux variations.The remaining 16 blazars were found to be non-variable, however, some of them have exhibited minor TeV flaring activities in the time period covered in this work.Overall, these results lead to the conclusion that blazars exhibit long-term flux variations at TeV energies.Another important finding of this work is the observation of complex variability patterns when comparing the TeV and GeV band light curves.In some cases, e.g., Mkn 421, the TeV flaring events had weak/no GeV counterpart and the GeV emission peaked before/after the peak of the TeV flaring activity in some other sources.Admittedly, a strong claim could not be made due to sparse data points in the TeV band light curves for the majority of sources.Furthermore, in the standard one-zone leptonic emission model, the VHE and hard X-ray radiations produced by the jet are expected to be originated from the same electron population of the highest energies.Indeed, the TeV band variability patterns observed from the Swift-BAT detected blazars present in our sample were often found to be similar to that identified in the 14−195 keV, monthlybinned, light curves indicating co-spatiality of the produced emission.Interestingly, the TeV band light curve of the bright blazar Mkn 421 revealed a flare that had no hard X-ray coun-terpart.These findings provide clues about the complex radiative mechanisms operating in blazar jets and are aligned with the results obtained in several recent multiwavelength campaigns focused on TeV blazars where complex variability and correlation patterns have been found (cf.Acciari et al. 2020;MAGIC Collaboration et al. 2021).In addition, further temporal and spectral investigation of these flaring episodes may also enable one to probe the radiative processes and possible interaction of the VHE emission with the extragalactic background light (e.g., Domínguez & Ajello 2015;Biteau et al. 2020).All in all, the findings presented in this work should be considered as the first step to unravel the long-term VHE variability behavior of blazars, which will be explored in greater detail with the upcoming Cherenkov Telescope Array due to its unprecedented sensitivity.

Figure 1 .
Figure 1.Three-months binned γ-ray light curves of VHE detected blazars studied in this work.In each figure, the top panel shows the flux variations in the 0.05−2 TeV energy range, whereas, the 0.1−50 GeV light curve is plotted in the bottom panel.Upper limits are derived at 95% confidence level and plotted with downward arrows.The grey bars represent the TS value in each time bin.In both panels, the blue dashed and the green solid lines represent the ∼15 years averaged γ-ray flux and Bayesian blocks, respectively, for the source under consideration.The vertical cyan lines highlight the periods of detection with Cherenkov telescopes.The hatched areas refer to the periods for which the z-DCF analysis was carried out.

Table 2 .
The z-DCF analysis results.The column information are as follows: (1) 4FGL name of the source; (2) time period considered for the z-DCF analysis, in MJD; (3) the z-DCF coefficient value; (4) the time lag between the two light curves in days with a positive value indicates the GeV band flux to lag behind TeV band flux, in days.