1. INTRODUCTION

Blazars are the most variable and luminous type of active galactic nuclei (AGN). They are among the most powerful emitters in the Universe, and they are thought to be powered by material falling onto a supermassive 106-1010M⊙ black hole (BH) at the center of the host galaxy through an accretion disc. The strong, non-thermal electromagnetic emission is generally detected in all observable bands, from radio to gamma-rays, and is dominated by a relativistic jet pointing in the direction of the observer (when the jet is not pointing in this direction, they are catalogued as radio galaxies). Blazars thus reveal the energetic processes occurring in the center of active galaxies. The most compelling mechanism 251

for the emission in the range from radio to optical is synchrotron emission, while gamma-rays are produced via inverse Compton (IC) by the same relativistic electrons producing the synchrotron emission (^{Sikora et al. 1994}; ^{Sikora et al. 2009}).

Many blazars were first identified as irregular variable stars in our own Galaxy. Their luminosity changes over periods that range from days to years, but without a pattern. With the development of radio astronomy many bright radio sources were discovered, most of them having optical counterparts. This led to the discovery of quasars and blazars in the late 50’s. The first redshift measurement was done for Blazar 3C 273 (^{Schmidt 1963}). The Fermi satellite, launched in 2008, paved the way for the beginning of a new blazar era thanks to its Large Area Telescope (LAT), which detects sources in the 0.1 - 300 GeV energy range (^{Atwood et al. 2009}) with much higher sensitivity than its predecessors. Blazars present high variability in this energy range, with periods in which the flux increases considerably (^{Abdo et al. 2009}). Detailed studies of blazar spectra, and in particular their spectral variability (^{Urry 1996a},^b), are fundamental tools to determine the physical processes responsible for particle acceleration and emission in the jet.

Blazars are often divided into two subclasses: BL Lacs (^{Stickel et al. 1991}) and flat spectrum radio quasars (FSRQs) (^{Angel & Stockman 1980}), although up to date there is no complete agreement in the selection criteria. BL Lacs are named after their prototype, BL Lacertae. They are characterized by rapid, large-amplitude flux variability and significant optical polarization. One of their defining features is the weakness or absence of emission lines; therefore, their redshifts can only be determined from features in their host galaxies’ spectra. Besides, their spectra are dominated by a relatively featureless nonthermal emission continuum over the entire electromagnetic range. The nearby BL Lacs, which are associated to elliptical galaxies with typical absorption spectra, sometimes have narrow emission lines, while the more distant BL Lacs, whose host galaxies have never been detected, show broad emission lines. All known BL Lacs are associated with core-dominated radio sources.

Another type of source are flat spectrum radio quasars (FSRQ). They are usually more distant, more luminous, and have stronger emission lines than BL Lacs but, due to their similar continuum properties, they are collectively called blazars, together with BL Lacs. There are three suggestions regarding their classification: (i) FSRQ evolve into BL Lacs; (ii) they are different manifestations of the same physical object; and (iii) BL Lacs are gravitationally micro-lensed FSRQ (^{Urry & Padovani 1995}). An example of FSRQ is the source 3C 279 (also known as PKS 1253-05). It presents very high variability in the visible, radio and X-ray bands. It is also one of the brightest sources in the gamma-ray sky monitored by the Fermi Space Telescope.

Blazars have been strong candidate sources of astrophysical neutrinos since 2017, when IceCube reported the first neutrino EHE event (IC-170922A) (^{Kopper & Blaufuss 2017}) with energy ≈ 290 TeV associated with a flare from Blazar TXS 0506+056 (^{Aartsen et al. 2018}). This opened the possibility to establish a link between the neutrino signal and the electromagnetic emission at different energy ranges, which is of key importance in the search for the privileged target of hadronic interaction processes, as well as for a pion decay component at VHE.

However, the blazar contribution to the diffuse astrophysical flux measured by IceCube is still a matter of debate (^{Aartsen et al. 2017}). The diffuse Galactic contribution (^{Gaggero et al. 2015}) is constrained by recent analyses from IceCube and ANTARES (^{Aartsen et al. 2017a}; ^{Albert et al. 2017}) to just 8.5% of the full sky measured neutrino flux, and large room remains for other AGNs (^{Atoyan & Dermer 2001}; ^{Padovani et al. 2016}; ^{Murase & Waxman 2016}) and starburst galaxies (^{Tamborra et al. 2014}; ^{Peretti et al. 2020}). In this work we investigate whether the long term observation, which is mainly characterized by the quiescent states, is well described by a leptonic scenario. On the other hand, by concentrating on a possible hadronic component during the flaring states, we remark the importance of the estimation of a DC factor to obtain long term neutrino expectations.

In this paper we analyze a sample of 38 bright blazars: 19 FSRQ and 19 BL Lacs from the Fermi 3FGL catalog (^{Acero et al. 2015}), through 10 years of data (from 2008 to 2018). We made a first selection of the blazars with the highest flux values and, from them, we selected the ones with ‘good’ LC (not dominated by upper limits). The paper is organized as follows: in § 2 we analyse the LC for each source and identify the active states. Following two different procedures, we distinguish the steady-state phase from the active phase for each of the sources of the sample. In § 3 we calculate the duty cycle (DC) according to the baseline flux obtained with each procedure. In § 4 we discuss a possible negative correlation between the blazar’s synchrotron peak frequency and the gamma-ray luminosity. In § 5 we consider the DC value and the blazar sequence trend from § 4 as useful tools to better understand the physical processes at play in the potential neutrino sources. In § 6 we show our results. Finally, in § 7 we present our discussion and conclusions.

2. DATA ANALYSIS

We selected 38 of the blazars with highest flux values from the Fermi 3FGL catalog (^{Acero et al. 2015}) as listed in Table 1. Among these, there are 2 with no reported redshift. Figure 1 shows the redshift distribution of the sample. The redshift and spectral type of each source were taken from ^{Ackermann et al. (2015)}. In order to perform the data analysis we made use of the tool Enrico (^{Sanchez & Deil 2013}) from the Fermi Science Tools (FSSC) software, version v11r5p3¹⁰ (^{Asercion 2018}). The analysis was performed in the energy range of 0.1 - 300GeV, in a RoI (region of interest) with radius of 20^◦, considering the isotropic and Galactic diffuse emission components (iso P8R2 SOURCE v6 v06.txt and gll iem v06.fits¹¹) falling within the RoI. For the cases in which we considered the extragalactic background light (EBL) absorption we used the model by ^{Franceschini, Rodighiero, & Vaccari (2008)} as reference.

Fig. 1 Redshift distribution of the sample of blazars considered in this work. The color figure can be viewed online. 

We generated the LC for the 10-year period with two-month, one-month, two-week and one-week time bins. From these, we selected the two-week binning, since it allowed us to see periods of activity with good resolution while minimizing the number of upper limits (UL). Figure 2 shows the LC for the sources 3C 66A, a BL Lac in the Andromeda constellation, and Mkn 421, a BL Lac located in Ursa Major. We consider them as benchmark sources since they have been intensively studied in the last decade, they are relatively nearby sources [(z = 0.444 (^{Miller et al. 1978}) and z = 0.03 (^{de Vaucouleurs et al. 1991}), respectively)] and flux and spectral information is available over the entire electromagnetic spectrum.

Fig. 2 Flux light curves (considering EBL absorption) for the sources 3C 66A and Mkn 421 using a two-week binning. The red dashed-dotted line and the black solid line are the threshold flux F _thr calculated according to Vercellone’s criterion and the µ _G + 3σ criterion (Tluczykont). The color figure can be viewed online. 

3. DUTY CYCLE

The duty cycle is the level of activity of a source. It is obtained by calculating the fraction of time that the source spends in the active phase (or flaring state, with the flux above a given threshold) with respect to the total observation period T _obs as a function of the flux:

Here, t _i is each of the time intervals the source spends in the flaring state (^{Vercellone et al. 2004}). The baseline phase (or steady state) is defined as the fraction of time the flux remains below a given threshold.

In order to separate the active phase from the steady-state phase, we need to define a threshold value for the flux. There are different criteria to do this.

3.2. 3σ or Tluczykont Criterion

According to ^{Tluczykont et al. (2010)}), we organize the data into a flux histogram, as shown in Figure 3. A flux state is defined as the integrated flux value in each time bin, obtained from the source LC. We then organize the data into a flux histogram (flux vs. number of events with that value of the flux) and fit them with the sum of a Gaussian and a Log-Normal distribution¹², f _G + f _LN . The Gaussian,

fG(x)=AGσG2πexp(-(x-μG)22σG2), (4)

represents the baseline state, and the Log-Normal,

fLN(x)=ALNxσLN2πexp(-(log(x)-μLN)22σLN2), (5)

describes the active-state phases. The parameters µ _G , σ _G and A _G are the mean, standard deviation and normalization for the Gaussian distribution, while µ _LN , σ _LN and A _LN are the mean, standard deviation and normalization for the Log-Normal distribution. The parameters of the fits of the entire sample can be found in Appendix A. According to Tluczykont’s criterion, the source threshold flux F _thr is a variable in the range from 3σ _G to µ _G + 3σ _G . We define the flux F _baseline as the mean value of all the non-flaring states:

Fbaseline=∫0FthrxfG(x)dx∫0FthrfG(x)dx. (6)

Fig. 3 Flux distribution for the sources 3C 66A and Mkn 421 considering EBL absorption. The black dotted line is the fit with the Log-Normal distribution, the black solid line is the fit with the Gaussian distribution, the red solid line is the fit with the sum of both distributions and the blue dotted line represents the threshold value. The color figure can be viewed online. 

If the flux is > F _thr the source is considered to be in the active phase. We recall that

Tbaseline=Tobs-Tflare, (7)

with Tflare = Pi ti.

Following the work of ^{Abdo et al. (2014)}, we calculated the average flare flux

⟨Fflare⟩=∫FthrFmaxxfLN(x)dx∫FthrFmaxfLN(x)dx (8)

as a function of F _thr and F _max, the maximum flux detected.

In order to obtain the average flux of the flaring states we need to replace F _flare by hF _flare i. Substituting these quantities in equation (3) we obtain

DC=TflareT[obs]=F¯-FbaselineFflare-Fbaseline. (9)

Figure 4 shows the plots of the DC obtained for 3C 66A and Mkn 421. In particular, the value of the DC obtained for Mrk 421 (av DC= 24.7%) is consistent with the results reported by ^{Patricelli et al. (2014)} and ^{Abdo et al. (2014)} using 3 year Milagro data in the TeV energy range. The DC plots for the whole sample are shown in Appendix B.

Fig. 4 DC vs. flux with (left panels) and without (right panels) EBL absorption, for 3C 66A (upper panels) and Mkn 421 (lower panels). The blue (black) lines represent the DC range inferred from Tluczykont (Vercellone’s) criterion. The color figure can be viewed online. 

The threshold flux calculated for 3C 66A is F _thr = 1.4 × 10⁻⁷ phcm⁻² s⁻¹ (Tluczykont criterion, black solid line in Figure 2) and F _thr = 1.06 × 10⁻⁷ phcm⁻² s⁻¹ (Vercellone’s criterion, red dashed-dotted line in Figure 2). In the case of Mkn 421, F _thr = 3.87 × 10⁻⁷ phcm⁻² s⁻¹ (Tluczykont criterion) and F _thr = 3.52 × 10⁻⁷ phcm⁻² s⁻¹ (Vercellone’s criterion). We recall that these values of F _thrhave been calculated considering EBL absorption.

4. BLAZAR SEQUENCE

^{Fossati et al. (1998)}, ^{Kubo et al. (1998)} and ^{Donato et al. (2001)} observed a negative correlation between the peak frequency of the synchrotron emission from blazars and their bolometric luminosities (traced through their radio luminosities in their observations). It was named the observational blazar sequence, and was thought to be an evidence in favor of leptonic emission from the source. This was later revealed to be a result of selection bias, although the existence of a physical blazar sequence (^{Ghisellini et al. 1998}) was explained due to the different cooling experienced by the electrons under different environments in BL Lacs and FSRQs.

In the case of blazars, the gamma-ray luminosity can be a better indicator of the integrated bolometric luminosity than the radio luminosity (^{Fossati et al. 1998}; ^{Wang, Z. et al. 2017}). It is thus possible to explore this relationship between the two variables for our sample, using Fermi-LAT data accumulated over 10 years. Figure 5 shows the integrated gamma-ray luminosity (L _γ ) in the 0.1 - 300 GeV range, averaged over 10 years of observation, vs. the synchrotron peak frequency (ν ^s _peak) for each source. The value of ν ^s _peak has been obtained from the 3FHL catalog (^{Ajello et al. 2017}) when available. Luminosities have been calculated assuming a ΛCDM Universe with standard cosmological parameters (H ₀ = 67.8, Ω _m = 0.308, Ω _λ = 0.692).

Fig. 5 Correlation between the synchrotron peak frequency and gamma-ray luminosity (0.1 - 300 GeV) for our sample. BL Lacs are shown in blue and FSRQs in red. The black line represents the best fit for the combined sub-sets. The color figure can be viewed online. 

The two sub-samples of BL Lacs and FSRQs show a negative correlation trend, the FSRQs having higher luminosity and lower peak frequency on average as compared to the set of BL Lacs. The combined sample of BL Lacs and FSRQs has a

Pearson coefficient of −0.7, with a best-fit slope of −0.72 ± 0.10. This moderate correlation trend can be interpreted in the context of leptonic emission favored in the inner part of the jet, providing information on the sub-PeV neutrino emission from this region (^{Murase et al. 2014}). Other authors (^{Finke 2013}) argue that there is a correlation between Compton dominance and the peak frequency of the synchrotron component, which is also consistent with the blazar sequence.

5. THE MULTIMESSENGER CASE AND VERY HIGH ENERGY NEUTRINO EMISSION

In this work, two important multimessenger aspects of gamma-ray emission from blazars were analyzed: the DC value and the blazar sequence trend. In particular, they can be used to better understand the physical processes responsible for the gammaray emission and to constrain a possible VHE neutrino counterpart. The first parameter represents a duty factor of the blazar activity, and it is considered a crucial factor whenever a long term unblinded analysis is done with a neutrino telescope, taking into account one or more VHE emitting blazars. In particular, under the assumption that the hadronic emission can be amplified during the VHE flaring states, it is worth using this value to set the expected neutrino observations. Figure 6 shows the value of average DCs of the samples obtained with 10 years of Fermi-LAT data. Taking them into account, the known ULs of full-sky neutrino emission expected from blazars can change by a factor > 2. In fact, the full sky upper limit reported by IceCube collaboration for blazars considering the unblinded point-like searches (^{Aartsen et al. 2017}) does not consider the DC for the weighting factor (w _γ/ν ) when deducing possible neutrino flux from the gamma-ray emission.

Fig. 6 DC distribution of our samples according to Tluczykont’s (left panels) and Vercellone’s (right panels) criteria, considering EBL absortion (upper panels) and not considering it (lower panels). The color figure can be viewed online. 

In more detail, the contribution to the diffuse astrophysical neutrino flux observed by IceCube (^{Aartsen et al. 2020}, ²⁰¹⁶, ^2017b) attributed to blazars will decrease from 27% to 7% in ^{Aartsen et al. (2017)}, for energies > 10 TeV, and from 40-80% to 10-20% in ^{IceCube Collaboration et al. (2018)} for energies > 200 TeV, when an average DC of 26% is considered. This decreasing factor on the expected neutrino flux will be quite independent of the blazar gamma-ray luminosity reached during the flaring periods even though the possibility to observe it during such activity period is strictly dependent on that, under the assumption of an amplified hadronic emission of active states.

While the DC can represent a temporal constraint on the hadronic activity from BL Lacs and FSRQs, the blazar sequence trend can constrain the emitting process as well as the energy range of VHE neutrinos produced, limiting the production of subPeV neutrinos (^{Murase et al. 2014}) by the interaction of jet accelerated protons with the ambient photon fields. The trend shown in Figure 5 can be explained assuming that a self-synchrotron Compton (SSC) emission scenario takes place in the inner part of the jet. Assuming the results of ^{Murase et al. (2014)}, Figure 5 favors FSRQs as opposed to BL Lacs as the main candidates for neutrino production in the range 1-100 PeV. Moreover, it is worth highlighting the fact that the blazar sequence reported here takes into account the full Fermi-LAT data set, which is dominated by the quiescent periods and does not exclude a different emission scenario during the flaring periods. If that were the case, Figure 7 shows that we cannot discriminate which class of blazars, BL Lacs or FSRQs, can produce a higher neutrino fluence over a longer observational period.

Fig. 7 Relation between flare duration (T _flare ), duty cycle (DC) and flare luminosity (L _ff ) for BL Lac type (left panel) and FSRQ type blazars (right panel). The color figure can be viewed online. 

We searched for possible correlations between these variables. However, no correlation was found. Furthermore, we found no correlation between DC and the source redshift, independently of the criterion used.

6. RESULTS

Tables 2 and 3 show the DC values for all the sources with and without EBL absorption, according to Tluczykonts and Vercellone’s criteria. We calculated the DC in the 0.1-300 GeV energy range for ≈ 10 years of Fermi-LAT data (2008 - 2018). The table includes: for the Tluczykont criterion, the DC obtained considering the maximum baseline flux (DC _bl %), the DC obtained assuming 0 as the baseline flux (DC₀%), and the average DC for each source (DC _av %). For Vercellone’s criterion, the DC value and its error. The distribution of the values presented in the tables is shown in Figure 6.

Following the Tluczykont criterion, the average DC for the sample when considering EBL absorption (and not) is 22.74% (23.08%). According to Vercellone’s criterion we obtained values of 20.09% (22.27%). The blazars PKS 0426-380 and BL Lacertae have the highest DC values (considering EBL absorption), with an average DC of ≈ 42% and 43%, respectively (Tluczykont’s criterion). On the other

hand, the lowest value of DC obtained corresponds to Blazar S5 0716+71 with a DC _av value of ≈ 8% considering EBL absorption (Tluczykont’s criterion).

At this point, we must consider the following: choosing a flare flux threshold in terms of the time average flux (Vercellone’s criterion) does not allow a direct comparison of DC between sources to be made, because the time average flux is influenced by the level of activity. In the case of a highly active source, the time average flux, and consequently its flare flux threshold, would be much higher than the baseline flux, so the DC only refers to the highest flux states. On the contrary, for a less active source the time average flux is close to the baseline flux, so the DC refers to almost all the flaring states. Since the DC of an active source includes only the highest flux-flaring states, it is possible to obtain a DC value smaller than the one for a less active source and erroneously conclude that the latter is more active than the former. In general, Tluczykont criterion leads to a higher DC average value than Vercellone’s criterion.

Figure 7 shows the average DC and flare luminosity (Tluczykont criterion) as a function of the major flare duration, together with flare luminosity vs DC. As seen from the figure, there is no significant difference between the values of DC for BL Lacs and FSRQ. We selected the major flares as the ones presenting the greatest number of consecutive temporal bins with the highest fluence in the 10-year dataset.

7. DISCUSSION AND CONCLUSION

We have presented the results of 10 years of FermiLAT data in the energy range of 0.1-300 GeV for a sample of 38 bright blazars following the two criteria presented in § 3. In particular, we have estimated the gamma-ray variability of the BL Lacs and FSRQs present in the sample, as well as the gamma-ray fluence during the activity periods. We have searched for the possible origin of this emission when looking at the entire observation time, dominated by quiescent periods, advocating a more enhanced hadronic component during the flaring periods.

The average DC for the whole sample following Tluczykont’s criterion varies in the range [≈8-43]% considering EBL absorption. This implies that blazars present a wide range of activity levels. Ongoing work aims to make a complete spectral analysis of these sources. When correlating the gamma-ray activity to the hadronic emission from the sources, the obtained DC becomes relevant also for the calculation of the expected neutrino signal over a long observational period, as discussed in ^{Marinelli et al. (2019)}. In particular, the average DC values reported in Figure 7 can be used as a scaling factor for the expected neutrino fluence. It is noteworthy, however, that the variance of the DC inferred considering EBL and not considering it is negligible, which would imply that the DC is not affected by the EBL absorption at these energies and distances. Even though the sample is small, a quick analysis shows that there is no correlation between DC and redshift (nor between DC and luminosity) for this sample of blazars. However, further analysis is required in order to make a stronger statement. Finally, we find that, for this sample, there is no notorious difference between the DC values inferred for FSRQ and BL Lacs (see Figure 7).