GeV Variability Properties of TeV Blazars Detected by Fermi-LAT

Variability is a prominent observational feature of blazars. The high-energy radiation mechanism of jets has always been important but is still unclear. In this work, we performed a detailed analysis using Fermi-LAT data across 15 yr and obtained GeV light-curve information for 78 TeV blazars detected by Fermi. We provided annual GeV fluxes and corresponding spectral indices for the 78 TeV blazars and thorough monthly GeV fluxes for a subsample of 41 bright blazars. Our results suggest a strong correlation between the γ-ray photon index and logLγ for the flat spectrum radio quasars (FSRQs) and high-energy peaked BL Lacs. Fourteen sources in our sample show significant GeV outbursts/flares above the relatively stable, low-flux light curve, with six of them showing a clear sharp peak profile in their 5 day binned light curves. We quantified the variability utilizing the fractional variability parameter F var, and found that the flux of the FSRQs showed significantly stronger variability than that of the BL Lacs. The 41 bright blazars in this work are best fit by a log-normal flux distribution. We checked the spectral behavior and found 11 out of the 14 sources show a bluer-when-brighter trend, suggesting this spectral behavior for these TeV blazars at the GeV band arises from the mechanism in which the synchrotron-self Compton process dominates the GeV emission. Our research offers a systematic analysis of the GeV variability properties of TeV blazars and serves as a helpful resource for further associated blazar studies.


INTRODUCTION
Blazars are the most powerful active galactic nuclei (AGNs) sources, which show very extreme observational properties, including variability over almost the whole electromagnetic waveband, high and variable polarization, strong γ-ray emissions, and apparent superluminal motion, which are believed to be associated with a relativistic beaming effect of the jet (Urry & Padovani 1995;Villata et al. 2006;Fan et al. 2014;Gupta et al. 2016;Xiao et al. 2019Xiao et al. , 2022b;;Fan et al. 2021).Blazars are usually divided into two subclasses: flat spectrum radio quasars (FSRQs) with strong emission lines, and BL Lac objects (BL Lacs) that have weak or even no emission lines (Scarpa & Falomo 1997).The broadband emission of blazars ranges from radio band to very-high-energy (VHE), which is generally dominated by non-thermal radiation.The broadband spectral energy distribution (SED) of blazar shows two humps, which is generally accepted that the lower energy hump peak is dominated by the synchrotron mechanism.The higher energy hump peak could be produced by inverse Compton (IC) scattering of synchrotron photons (synchrotron-self Compton, SSC, Bloom & Marscher 1996;Finke et al. 2008) and external photons (external Compton, EC, Sikora et al. 1994;Kang et al. 2014) in the framework of leptonic models.Meanwhile, the hadronic model suggests that the higher energy hump is attributed to the proton synchrotron radiation and secondary particle cascade (Mücke & Protheroe 2001;Dimitrakoudis et al. 2012;Diltz et al. 2015;Cerruti et al. 2015;Xue et al. 2021;Wang et al. 2022c).The hadronic model seems to be promising following the detection of extragalactic neutrino events from the blazar TXS 0506+056 (IceCube Collaboration et al. 2018a).
The discovery of the first TeV blazar, Mrk 421, was a surprise when it was detected by the Whipple telescope in 1992 because it was barely seen in the γ-ray band (Punch et al. 1992).The following detection of more TeV blazars, e.g., Mrk 501, 1ES 2344+514, PKS 2155-304 (Quinn et al. 1996;Catanese et al. 1998;Chadwick et al. 1999), started the era of the TeV blazar study.There are 252 sources associated with TeV emission, and 81 of them are confirmed as blazars according to TeVCat1 , the detection of TeV emissions mainly relies on ground-based Cherenkov telescopes, e.g., Major Atmospheric Gamma Imaging Cherenkov Telescopes (MAGIC), High Energy Stereoscopic System telescopes (H.E.S.S), Very Energetic Radiation Imaging Telescope Array System (VERITAS).Our understanding of blazar TeV emission is limited by several issues, e.g., the sample size of TeV blazars, the lack of TeV light curves due to the observation mode of Cherenkov telescopes, the effective absorption of extragalactic background light (EBL).A multi-wavelength study is usually employed to investigate the emission properties of TeV blazars, however, this method can only be applied to several individual sources.Otherwise, we can also study this subject at other bands, for instance, the GeV γ-ray band.
Fermi -LAT, which has been launched since 2008, scans the entire sky every three hours ranging from 20 MeV to above 300 GeV (Atwood et al. 2009).During the last 15 years, there are 5 generations of released Fermi catalogs with the latest one being the fourth Fermi -LAT source catalog (4FGL, Abdollahi et al. 2020).More than 5000 sources have been observed, about 60% of which are confirmed as blazars and blazars have been established to be the dominant γ-ray sources in the extragalactic sky (Ackermann et al. 2015a;Ajello et al. 2020).Based on these observations, people made significant progress in blazar studies, e.g., the classification that depends on the synchrotron peak frequency (Abdo et al. 2010a;Fan et al. 2016;Yang et al. 2022), the 'blazar sequence' (Fan et al. 2017;Ghisellini et al. 2017;Ouyang et al. 2023), the blazar central engine (Paliya et al. 2021;Xiao et al. 2022a).More studies focus on individual sources, study the properties of flares or outbursts, and put constraints on the blazar emission mechanism, such as the flare of 3C 279 (Shukla & Mannheim 2020;Tolamatti et al. 2022;Wang et al. 2022a), the neutrino TXS 0506+056 (IceCube Collaboration et al. 2018b), variability and spectral properties for 3C 279, Ton 599, and PKS 1222+216 (Adams et al. 2022a), light curve study of PKS 1510+089 (Prince et al. 2017).To obtain information on blazar emission variability, periodicity, and spectrum.Long-coverage observations on different timescales and spectral analysis can be carried out by taking advantage of the all-sky monitoring capabilities of Fermi -LAT.Recently, the Fermi -LAT light curve repository (LCR), which provides a publicly available, continually updated library of gamma-ray light curves of Fermi sources, is released (Abdollahi et al. 2023).However, this library provides light curves binned only on timescales of 3 days, 7 days, and 30 days based on the 10-year Fermi -LAT point source (4FGL-DR2) catalog (Ballet et al. 2020).
In this work, we aim to provide detailed GeV γ-ray variability information for the TeV blazars based on 15-year 4FGL-DR3 data.We described the sample selection and Fermi data analysis in Section 2. The results are reported in Section 3. The discussion and conclusions are presented in Section 4 and Section 5, respectively.

Sample selection
We collected 78 blazars, including 66 BL Lacs, 8 FSRQs, and 4 blazar candidates of uncertain type (BCUs) by cross-matching TeVCat and the latest LAT 12-year source (4FGL-DR3) catalog (Abdollahi et al. 2022).These sources are listed in Table 1, in which columns (1) and ( 2) give source 4FGL name and associated name; column (3) gives the redshift obtained from Chen (2018); column (4) gives the classification that is determined based on the synchrotron peak frequency information and criterion in Fan et al. (2016); We also show the redshift distribution of each type of these blazars in Figure 1.

Fermi-LAT observations and data reduction
LAT is one of the main instruments on board Fermi.LAT scans the whole sky every three hours in the energy range from 20 MeV to >300 GeV (Atwood et al. 2009).We selected LAT data from the Fermi Pass 8 database in the time period from 2008 August 4 15:43:36 (UTC) to 2023 Mar 9 03:03:00 (UTC), with an energy range of 1-300 GeV.Following the recommendations of the LAT team2 , we selected events with zenith angles less than 90 deg to prevent possible contamination from the Earth's limb.The LAT science tool Fermitools 2.0.8 and the instrument response function (IRF) P8R3 SOURCE V2 were used.For the selected samples, a 20 • × 20 • square region of interest (ROI) centered at their positions given in 4FGL-DR3 was selected.The normalization parameters and spectral indices of the sources within 5 deg from the target, as well as sources within the ROI with variable index ≥ 72.44 (Acero et al. 2015), were set as free parameters.All other parameters were fixed at their catalog values in 4FGL-DR3.We used the original spectral models in 4FGL-DR3 for the sources in the source model when performing binned maximum likelihood analysis with gtlike.A simple power-law (dN/dE ∝ E −Γ , where Γ is the photon index) spectral type was used for each blazar when deriving its light curve.We checked through the likelihood analysis results assuming a power-law model comparing with a log-parabola (dN/dE ∝ (E/E 0 ) −α−β log(E/E0) , where α and β are spectral parameters) for the samples in the source models, and found that a log-parabola is not significantly preferred over a power-law for the samples, except for J0035.9+5950 and J0221.1+3556.Therefore we changed the two sources' spectral parameters accordingly in Table 1.The comparison was conducted by calculating −2 log(L pl /L logP ), where L pl and L logP are the maximum likelihood values obtained from a power-law and a log-parabola, respectively (Abdo et al. 2013).In addition, the background galactic and extragalactic diffuse emission models were added to the source model using the spectral model file gll iem v07.fits and iso P8R3 SOURCE V2 v1.txt, respectively.The normalizations of the two diffuse emission components were set as free parameters in the analysis.We constructed light curves binned in 90-day time intervals by performing standard binned maximum likelihood analysis, calculated flux (F γ ) and photon spectral index (Γ) for the energy range of 1-300 GeV spectrum and listed them in columns ( 5) and ( 6) of Table 1.

Annual and monthly intensity at the GeV band
Fermi -LAT has conducted observations at γ-ray energy bands over 15 years by scanning the whole sky every three hours.As we aim to provide detailed GeV spectral behaviors, study the GeV variability of the TeV blazars, and put constraints on the blazar emission model.We calculated the annual GeV fluxes and corresponding photon spectral indices for the 78 TeV blazars in our sample and listed them in Table 2, in which the annual (360-day time interval) maximum, the minimum, and the mean fluxes and the corresponding photon spectral indices are given for the past 15 years since the launch of Fermi (MJD 54683).

The GeV luminosity and spectral photon index
where is a luminosity distance (Komatsu et al. 2011) and (1 + z) (Γ−2) stands for a K-correction.We calculated the γ-ray luminosity of the 74 blazars that have redshift information, using the energy flux derived from the binned likelihood analysis, and studied the correlations between the GeV γ-ray luminosity and photon index in Figure 3.It is found that FSRQs occupy the upper-right region, then the IBLs occupy the middle region, and the HBLs occupy the lower-left region of Figure 3.This result suggests that the TeV blazars show a decrease in the GeV γ-ray luminosity and photon spectral index with the increase of synchrotron peak frequency, and indicates a 'blazar sequence' that was initially proposed by Fossati et al. (1998).In addition, we calculated the linear regressions between Γ and log L  3 and suggest a strong correlation between Γ and log L γ for the FSRQs and HBLs, while no correlation for IBLs.We also conducted statistics on weighted Kendall's tau (Shieh 1998) and Spearman's coefficients.The r value obtained through weighted Kendall's tau analysis was −0.74, −0.60, −0.28 for FSRQs, HBLs, and IBLs respectively.The coefficient obtained using Spearman's statistics is r = −0.86,p = 6.53 × 10 −3 , r = −0.59,p = 4.96 × 10 −5 and r = −0.26,p = 0.26 for FSRQs, HBLs and IBLs respectively.These results have supported that of Pearson's analysis mentioned above.Ackermann et al. (2015b); Ajello et al. (2020) showed the LAT photon index versus the γ-ray luminosity for the different blazar subclasses of the whole sample in the Third LAT AGN Catalog (3LAC) and the Fourth LAT AGN Catalog (4LAC) blazars.The trend of softer spectra with higher luminosity reported in previous catalogs is also confirmed.However, 4LAC noted that the correlation between photon index and γ-ray luminosity is significant overall for blazars, but much weaker when the different classes are taken independently.While 3LAC γ-ray luminosity results were computed from the 4-year Fermi -LAT point source (3FGL) catalog energy flux between 100 MeV and 100 GeV.They also mentioned that due to the bias in the selection criteria for the 57 BL Lacs with both lower and upper limits on their redshifts or only upper limits, the HBLs with both limits were found to be more luminous on average than those with measured redshifts.

GeV flux and spectral photon index in flares
We note that 28 of our TeV blazar samples were reported in the second Fermi all-sky variability analysis catalog (2FAV, Abdollahi et al. 2017).The analysis of 2FAV was ran in weekly time bins using the first 7.4 years of Fermi data in two independent energy bands, 100-800 MeV and 0.8-300 GeV.We have checked these light curves to find the GeV outbursts/flares that meet the criterion, which is that a source shows flare flux more than 10 times larger than its flux in quiescent states and the significance compared to the quiescent light curves is more than 5σ simultaneously.There are 14 sources that show significant outbursts/flares at the GeV band during the Fermi campaign, these sources are listed in Table 5.All of these significant flares have been reported in 2FAV, except for J1422.3+3223.The photon indices and fluxes of these 14 blazars with bright flares in the 1-300 GeV band are shown in Figure 6, for which only the flux data points with TS > 9 were selected for the plot.The insets in Figure 6 display the photon index resulting from an analysis where photons were sorted in five bins in 5-day flux, plotted versus the 5-day flux.Fluxes and photon indices during their flaring states are listed in columns ( 4) and (5) of Table 5, respectively.
For these 14 blazars with bright flares, we constructed light curves in the 5-day bin in their flaring states.There are 6 blazars that showed a clear single sharp peak profile contained in the flare that meet the criterion, which is that its flare flux shows more than 12 times larger than its flux in quiescence and the significance compared to the quiescent states is more than 4σ simultaneously.We also searched intra-day flares and only found 4FGL J1256.1−0547(3C 279) had minute-scale variability in 2018, and this result has been reported in our previous work (Wang et al. 2022a).We determined the properties of the 6 single sharp peak cases by fitting their profiles with a formula given by where F c and F 0 are the constant flux and height of a peak, respectively, t 0 is the flux peak time, T r and T d are used to measure the rise and decay time in units of day.We show the flare profiles in Figure 4 for flares with a single sharp peak, and the distribution of rise and decay time in Figure 5.We can calculate the flare asymmetry parameter following Chatterjee et al. (2012) as the results are listed in column (9) of Table 5,   We have checked the coincidence between Fermi -LAT GeV detections and TeV detections of the sample.There are 22 sources in Table 1 that have been detected TeV emission during the flaring states observed by Fermi -LAT , and 56 sources were in low-state.J0509.4+0542(TXS 0506+056) was detected at VHE by MAGIC and VERITAS (Ansoldi et al. 2018;Abeysekara et al. 2018;Acciari et al. 2022).It was in an active flaring state around the arrival of the highenergy neutrino IceCube-170922A (IceCube Collaboration et al. 2018b).While Garrappa et al. (2019) found another blazar GB6 J1040+0617, in spatial coincidence with a neutrino in this sample and the chance probability of 30% after trial correction, indicating the source of this neutrino remains unknown.J1015.0+4926 was detected in a flaring state at VHE by MAGIC during February−March 2014 (Ahnen et al. 2016), Fermi -LAT observation was coincident with the TeV detection, and the GeV flux reached a level of 6.5 times higher than its low-state.J1058.6+2817,Fermi -LAT and MAGIC successively reported its flaring activity during March−April 2021 (Angioni 2021;Blanch 2021).J1217.9+3007, its multiwavelength observations with VERITAS and Fermi -LAT showed a well-connected high flux state in February 2014 (Abeysekara et al. 2017a).J1728.3+5013,Archambault et al. (2015) reported the first detection of γ-ray flaring activity at VHE from this blazar, the flaring flux is about five times higher than its low-state.Fermi -LAT detected this source with mild flare and it was observed a photon of energy more than 300 GeV as reported in MAGIC Collaboration et al. (2020a).VERITAS detected VHE emission from J1813.5+3144 with the similar flux      where S 2 is the variance of the flux, σ 2 err is the mean square value of uncertainties, ⟨F γ ⟩ is the mean photon flux.Negative values of F var indicate very small or absent variability and/or slightly overestimated errors.We derived the mean values of F var are 1.54 ± 0.02, 0.12 ± 0.15, 0.65 ± 0.06, 1.07 ± 0.04 for the FSRQs, HBLs, IBLs, and LBLs, respectively.The resulting values indicate that the flux of the FSRQs showed significantly stronger variability than that of the BL Lacs.As the synchrotron peak frequency decreases, the F var value generally becomes larger.Here we presented a histogram of F var values for the FSRQs, HBLs, IBLs, and LBLs in Figure 7. Bhatta & Dhital (2020) presented an analysis of a sample of 20 powerful blazars (12 BL Lacs and 8 FSRQs) with 10 yr Fermi -LAT data, they obtained that the mean F var value of BL Lacs is 0.58 and that of the FSRQs is 0.96.The results show that in general FSRQs are more variable than BL Lac sources in their sample, which is compatible with ours.Similar future studies involving larger samples should be carried out for a stronger conclusion.For the individual source, our result of S5 0716+714 is consistent with that reported in Bhatta et al. (2016), the F var values are 0.65, 0.57, 0.58, 0.53 for BVRI filters versus our 0.59.
Besides, we found 14 TeV blazars (8 FSRQs, 1 LBL, 4 IBLs, and 1 BCU) with outbursts/flares, 6 out of the 14 flares show sharp peak profiles in flares.Based on the sharp peak profiles, we notice 4FGL J0303.4−2407 and 4FGL J0739.2+0137show a fast-rise-slow-decay subflare.This asymmetry can be related to the particle acceleration mechanism in the jet, a fast rise could result from an effective particle acceleration at the shock front and slow decay may be interpreted as the weakening of the shock (Sokolov et al. 2004;Tolamatti et al. 2022) or from the injection of energetic particles on a shorter timescale than the cooling process timescales (Acharyya et al. 2021).While 4FGL J1751.5+0938shows a slow-rise-fast-decay subflare, which may be associated with an efficient cooling process.

Flux distributions
The analysis of flux distribution helps us to determine whether the variability is caused by multiplicative or additive mechanisms.Evidence for log-normality in blazars in γ-ray on different timescales has been reported for different sources (e.g., Kushwaha et al. 2017;Sinha et al. 2017;Bhatta & Dhital 2020).Similarly, the log-normal flux distribution of blazars was seen in 3LAC (Ackermann et al. 2015b)).Shah et al. (2018) studied the flux distribution features of the selected 38 brightest Fermi blazars using the data collected during more than 8 years and found that the flux distribution for 35 blazars supports a log-normal distribution, implying a multiplicative perturbation linked with the emission process.Using a large sample of 1414 variable blazars from the Fermi -LAT LCR catalog, Wang et al. (2023) thoroughly investigated the γ-ray flux distribution and statistical properties, and compared the flux distributions with normal and log-normal distributions.Their results showed that the probability of not rejecting log-normal is 42.05%.We constructed histograms of the observed LAT GeV flux and fitted them to two different probability density functions (PDFs), a normal distribution and a log-normal distribution, and compared the results of χ 2 .To ensure sufficient data points for fitting the flux distribution, we selected the 41 bright blazars mentioned in the subsection 3.1.According to the chi-squared values from the fit, our results show that all of the bright blazars support a log-normal distribution rather than a normal distribution, which is also consistent with the results of previous studies.As the consistency between TeV detection and LAT observation that discussed in the subsection 4.1, the TeV detections correspond to the outlier periods of the flux distribution.The parameters of the considered two distributions fitting results and the source flux histograms are shown in Table 6 and Figure 8, only several items are presented here.

Flare duty cycle
The flaring state lasts only a fraction of the observation.Here we define the flaring state when any of the light curve's flux points exceeds a certain threshold following the method in Yoshida et al. (2023), f th γ , which is given by where f q γ is the quiescent level of γ-ray fluxes, f err γ is the average uncertainty of the γ-ray fluxes, and s denotes the significance above the quiescent level in standard deviation units of σ.Here we use s = 6 in this work, and the flaring threshold levels are plotted with dashed grey lines in Figure 2. From the light curves, we calculated the flare duty cycle (i.e., fraction of time spent in flaring states) for each flare.The flare duty cycle, is defined as where T tot is the total observation time, f γ is the γ-ray photon flux, and T is the time spent at the respective flux level.We find that our duty cycle results of the monthly-binned light curves show values ranging from 0.0 to 0.36, and there Abdo et al. (2010b) conducted an analysis of the first 11 months of the LAT Bright AGN sample (LBAS), and revealed that the average β values of the brightest 22 FSRQs and of the 6 brightest BL Lacs is 1.5 and 1.7, respectively.While Ackermann et al. (2011b) using 24 months of data and found the β value is ∼ 1.15 ± 0.10, which is somewhat flatter than the results deduced from the LBAS sample.Tarnopolski et al. (2020) presented a comprehensive analysis of the Fermi -LAT 10 yr long light curve modeling of 11 selected blazars by employing various methods.They found that the power-law slope index β calculated from the Fourier and LSP modeling falls in the range 1 ≲ β ≲ 2 mostly.Our results of PKS 1510−089, PKS 2155−304, and Mrk 421 are consistent with Sobolewska et al. (2014).They analyzed the γ-ray variability of 13 blazars with a linear superposition of OU processes, for which they found slopes mostly to be β ≲ 1. Prokhorov & Moraghan (2017) obtained β = 0.67 for PKS 2155−304, while we obtain β = 0.65 ± 0.03.Also, our result of 3C 279 is similar to the PSD slopes found by Meyer et al. (2019).Chatterjee et al. (2012) mentioned the average slope of the PSD in R-band of 6 blazars is similar to that found by the Fermi team, our result was in agreement for PKS 1510−089, but they obtained clearly steeper power-law fits than we did (2.3, 2.2 for 3C 279 and PKS 2155−304 versus our 0.75, 0.65).Compared to these recent results of selecting the several brightest sources, our PSD result at the GeV band is slightly flatter and has a larger range.The discrepancies can be caused by the difference in the analysis methods, different binning schemes, sampling interval, and total observation duration of the analyzed light curves or methods of their generation between the works.

The periodic behaviors
The periodogram of the light curves can be characterized by a single power-law PSD.However, if we closely observe the structures of the periodogram, we may occasionally find peaks at certain frequencies indicating the possible presence of (quasi-) periodic signals in the observations.The periodic oscillation in γ-ray band of blazar PG 1553+113 was reported by Ackermann et al. (2015c), this source is also contained in our sample and its light curve at the GeV band shows clear periodicity, and has been explained in mechanisms invokes a supermassive binary black hole system (Cavaliere et al. 2017;Sobacchi et al. 2017).Several studies have systematically searched γ-ray QPOs based on 3FGL (e.g., Prokhorov & Moraghan 2017;Peñil et al. 2020).Peñil et al. (2022) made a search for periodicity in a sample of 24 blazars by using 12 well-established methods applied to Fermi 12-year data, and found six out of the 24 sources show light curve periodicity with global significance greater than 3σ.Among our samples, some showed quasi-periodic oscillation (QPO) characteristics in their γ-ray light curves.There are 12 blazars have been reported to have γ-ray QPOs according to Table 2 in Wang et al. (2022b), while nearly 30 blazars have been reported to show possible QPOs with high-significance based on Fermi -LAT data so far.We note that various analysis methods can be affected by several caveats or effects that may have an impact when analysing time series, and lead to the overestimation of signal significance.The caveats remind us of the complexity of the QPO analysis in AGNs, and the importance of correction for trials when computing probabilities.Otero-Santos et al. ( 2023); Ren et al. (2023) provided a detailed discussion of some of the caveats.

The spectral behavior
Variability is one of the main characteristics of blazars, the variability time scale spans from years to hours and even to minutes.The variability of flux is always accompanied by the variation of spectra.The correlation between the spectral index and flux has been investigated for individual sources and also for large samples (Fiorucci et al. 2004;Gu et al. 2006;Dai et al. 2009;Bonning et al. 2012;Yuan et al. 2017;Raiteri et al. 2017;Meng et al. 2018;Xiong et al. 2020;Safna et al. 2020).In general, this correlation was mainly discussed at the optical band and demonstrates 'bluer-when-brighter (BWB)' behavior for BL Lacs, and shows 'redder-when-brighter (RWB)' behavior for FSRQs, except in some special cases e.g., 14 out of 29 Sloan Digital Sky Survey (SDSS) FSRQs show BWB trend (Gu & Ai 2011), 2 out of 40 Fermi FSRQs exhibit BWB trend and 7 out of 13 BL Lacs exhibit RWB trend (Zhang et al. 2022).Various models have been proposed to explain blazar optical spectral behavior, shock-in-jet model (Rani et al. 2010), two-components (one variable + one stable) or one synchrotron model (Fiorucci et al. 2004), the energy injection model (Spada et al. 2001;Zhang et al. 2002), and the also the vary of beaming effect (Larionov et al. 2010).Recently, Zhang et al. (2022) suggests a universal two component-model to interpret these two spectral behaviors, in which the observed optical emission of blazars consists of a stable or less-variable thermal emission component ( F ther ) primarily coming from the accretion disk, and a highly variable non-thermal emission component (F syn ) coming from the jet.The stronger the thermal emission component the bluer the color is, the weaker the thermal emission the redder the color is.
However, the spectral behavior at higher energy bands seems monochrome.We found a universal BWB trend at γ-ray band for the TeV blazars in our sample, especially the LBLs and FSRQs showing strong anti-correlation between the photon index and the GeV γ-ray luminosity.For the individual sources, Hayashida et al. (2012) performed a broadband study of 3C 279 flare and found BWB trend at the X-ray band and γ-ray bands.And this BWB trend was found again at the X-ray band for the same source during a phase of increased activity from 2013 December to 2014 April (Hayashida et al. 2015).Moreover, Aleksić et al. (2014b) made multi-frequency observations of PKS 1510-089 in early 2012 and reported a BWB trend at the X-ray band.Prince et al. (2017) studied the long-term light curve of PKS 1510-089 at GeV bands and reported the BWB trend during flares at different campaigns.There are 14 outbursts/flares of individual TeV blazars that have been analyzed and their spectral behavior has been illustrated in Figure 6.11 out of the 14 sources show the BWB trend, according to the insets in Figure 6, except 4FGL J0221.1+3556,4FGL J0303.4−2407, and 4FGL J1512.8−0906.We suggest this spectral behavior for blazars at the GeV band arises from the same mechanism, which is that the synchrotron-self Compton (SSC) process dominates the GeV emission for these TeV blazars.Considering the non-thermal electrons that produce the observed inverse Compton emission with an energy distribution where γ is the Lorentz factor of electrons, γ min and γ max are the minimum and the maximum values of Lorentz factor at the time of particle injection, N 0 is related to the total particle density N tot by N 0 = N tot (1 − α)/(γ 1−α max − γ 1−α min ), α is the electron spectral index.Then, the SSC emissivity (j ssc ) is related to the electron spectral index by j ssc (ϵ) ∼ ϵ −(2+α)/4 , ( where ϵ = hν/m e c 2 (Chiang & Böttcher 2002).From equation 9, we can see that the spectral behavior at the GeV band for blazars is mainly determined by the shape of the electron spectrum, which means a harder electron spectrum results in a corresponding harder emission spectrum.In this case, we can obtain the electron spectral index through the GeV γ-ray photon index via −(Γ − 1) = −(2 + α)/4 for the 11 outbursts/flares and list the results in column (6) of Table 5.

SUMMARY
This paper aims to provide detailed GeV variability of the TeV blazars and study the GeV spectral behaviors.We performed an analysis using the LAT data across 15 years and offered annual GeV fluxes and corresponding photon spectral indices for the 78 TeV blazars of our sample.We calculated the detailed monthly flux and corresponding photon spectral index of the 41 bright TeV blazars to further investigate the spectral behavior.A series of variability property analyses were conducted on the fractional variability, flux distribution, flare duty cycle, PSDs, and periodic properties.
Our main conclusions are as follows: (1) We investigated the possible correlation between GeV luminosity and spectral photon index.The results suggest a strong correlation between the log L γ and γ-ray photon index for the FSRQs and HBLs, while no correlation for the IBLs.
(2) There are 14 sources out of our sample that show significant flares, of which 6 exhibit a clear sharp peak profile in their 5-day binned light curves.4FGL J0303.4−2407 and 4FGL J0739.2+0137show a fast-rise-slow-decay subflare.This asymmetry can be related to the particle acceleration mechanism in the jet.While 4FGL J1751.5+0938shows a slow-rise-fast-decay subflare, which may be associated with an efficient cooling process.
(3) We quantified the variability utilizing the fractional variability parameter F var and the results indicate that the flux of the FSRQs showed significantly stronger variability than that of the BL Lacs.As the synchrotron peak frequency decreases, the F var value generally becomes larger.
(4) We constructed histograms of the observed GeV light curves and fitted them to two different PDFs, a normal distribution and a log-normal distribution.The results show that all of the bright sources in this work support a log-normal distribution.
(5) Our duty cycle results of the monthly-binned light curves show values ranging from 0.0 to 0.36, while the vast majority of the values are in the range of 0.0 to 0.2 except for three blazars.
(6) We found that the periodograms are consistent with a power-law form with the slope index β ranging between 0.22−1.98.Our PSD result at the GeV band is slightly flatter and has a larger range compared with the previous studies.In addition, 12 blazars in our sample have been reported to have high-significance γ-ray QPOs.
(7) Through checking the spectral behavior, we found 11 out of the 14 sources show a 'bluer-when-brighter' trend, which suggests this spectral behavior at the GeV band arises from the mechanism that the synchrotron-self Compton process dominates the GeV emission for these TeV blazars.
Thanks are given to the reviewer for the constructive comments and suggestions that helped us to make the paper more thorough.This work is supported by the National Natural Science Foundation of China (Grants Nos. 11975072, 11835009, and 11805031), the National SKA Program of China (Grants Nos.2022SKA0110200 and 2022SKA0110203), the 111 Project (Grant No. B16009), and the China Postdoctoral Science Foundation No. 2023M730523.

Figure 1 .
Figure1.The redshift distribution of each type of blazar in the sample.The histogram is illustrated in 5 bins, which are 0 ∼ 0.2, 0.2 ∼ 0.4, 0.4 ∼ 0.6, 0.6 ∼ 0.8, 0.8 ∼ 1.0.The blue bar stands for BCU, the orange bar stands for FSRQ, the green bar stands for HBL, the red bar stands for IBL, and the violet bar stands for LBL.

Table 2 .
Annual GeV fluxes and photon indices for the 78 TeV blazars of our sample.

Figure 2 .
Figure 2. The monthly binned light curves for the 41 bright blazars in our sample.Only six items are presented here.The complete figure set (41 images) is available in the online journal.

Figure 3 .
Figure 3.The correlation between 1-300 GeV photon index and luminosity of 74 blazars.The solid lines are the linear regressions fitting results.

Figure 4 .
Figure 4. Flare profiles for 6 TeV blazars with a single sharp peak in their 5-day binned light curves, an analytic function (dashed red curve) was used to fit the profile.
while k < 0 indicates a fast-rise exponential-decay (FRED) type flare.Approximately, k < −0.3 indicates faster rise than decay, k > 0.3 indicates faster decay than rise, while −0.3 < k < 0.3 indicates a symmetric profile, k = 0 for exactly symmetric flares.Among the 6 sharp peak flares, 4FGL J0303.4−2407 and 4FGL J0739.2+0137show FRED behavior, 4FGL J1751.5+0938shows the opposite, and the other three show symmetric profiles.Chatterjee et al. (2012) showed the distribution of the flare asymmetry parameter (k) for the optical and γ-ray flares with a sample of six blazars, which indicated that most of the flare profiles are symmetric at both wave bands.Abdo et al. (2010b) gave a systematic analysis of a larger sample of 106 objects by using the first 11 months of data of the Fermi survey and found only 4 sources with markedly asymmetric flares.
The connection with the TeV band

Figure 5 .
Figure 5. Distribution of rise and decay time for the 6 sharp peak flare profiles fitted to the data.

FluxFigure 6
Figure 6.1-300 GeV photon indices and fluxes of blazars with bright flares.Only data points with TS > 9 are plotted.The blue dashed horizontal lines indicate the average photon indices of those data.The insets show the photon index resulting from an analysis where photons were sorted in five bins using 5-day fluxes plotted vs. the 5-day flux (red points).

Figure 7 .
Figure 7. Distribution of the fractional variability Fvar for the light curves of FSRQs, HBLs, IBLs, and LBLs.

Figure 8 .
Figure 8. Flux distribution of bright blazars in our sample in the GeV band.The black and red curves correspond to normal and log-normal fits respectively.Only three items are presented here.The complete figure set (41 images) is available in the online journal.

Figure 9 .
Figure 9. PSD fits with power-law for LSPs of the bright blazars.The black curve is the raw LSP, and the red dashed line is the best fit.Only three items are presented here.The complete figure set (41 images) is available in the online journal.

Table 1 .
Sample of TeV blazars

Table 4 .
Monthly GeV photon indices of the 41 bright TeV blazars

Table 5 .
Fluxes, photon indices, fitting results for the flaring sharp peaks, and asymmetry parameters of the bright GeV flares of 14 TeV blazars The rise time fitting results for the flaring sharp peaks in units of day; (8): The decay time fitting results for the flaring sharp peaks in units of day; (9): The flare asymmetry parameter.

Table 6 .
Parameters of normal and log-normal distribution fitting for the γ-ray flux distribution of the Fermi-LAT sources.Here β slope gives the slope index result of the periodograms.Table6is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.