Flare Duty Cycle of Gamma-Ray Blazars and Implications for High-energy Neutrino Emission

We select 145 gamma-ray bright blazars (Table A1), consisting of 106 flat spectrum radio quasars (FSRQs), 31 BL Lac objects, and 8 blazar candidates of uncertain type (BCUs). Among them, the 144 blazars are selected from the Fermi LAT Monitored Source List, ⁸ which we identified as blazars with the Fourth Catalog of AGN detected by the LAT (4LAC-DR2, hereafter 4LAC; Ajello et al. 2020). The Fermi-LAT team monitors the light curves of a number of bright and transient sources that have flares during the mission to facilitate follow-up multiwavelength observations of flaring sources. When sources exceed the monitoring flux threshold of 1 × 10⁻⁶ cm⁻² s⁻¹, they are added to the list of monitored sources. ⁹ Since blazars are known for their intense flares, the list gives the appropriate blazar samples that have flaring fluxes with high statistics. In addition, we added one blazar of TXS 0506+056, whose flaring state is indicated to be in coincidence with IceCube-170922A (IceCube Collaboration et al. 2018a). Our sample comprises the brightest blazars in the 4LAC catalog as shown in Figure 1. In this figure, we show the time-average energy flux in the energy range of 0.1–100 GeV based on 10 yr of data as a function of redshift for our sample and all blazars in the 4LAC. In this work, we analyze data collected by the LAT for the period 2009–2018 as described in the next section.

Zoom In Zoom Out Reset image size

The Fermi LAT is a pair conversion telescope sensitive to gamma rays in the 20 MeV to >300 GeV energy range (Atwood et al. 2009). In this work, we analyze the LAT data from MJD 54682 to MJD 58739 using the Fermi Science Tools version 11-05-03 ¹⁰ and Fermipy version 0.17.3 (Wood et al. 2017). Following the standard selection recommendation, ¹¹ we select photons with energies between 100 MeV and 316 GeV that were detected within a radius of 20° centered on the position of the source of interest, while photons detected with a zenith angle larger than 90° with respect to the spacecraft were discarded to reduce the Earth-limb gamma-ray contamination.

The contribution from isotropic and Galactic diffuse backgrounds is modeled using the parameterization provided in the iso_P8R3_SOURCE_V2_v01.txt and gll_iem_v07.fits files, respectively, leaving their all parameters free to vary in the spectral fit. Sources in the 4FGL catalog within a radius of 15° from the source position are included in the model with their spectral parameters fixed to their catalog values (Abdollahi et al. 2020), while the normalization of sources within 3° is allowed to vary freely during the spectral fit. The spectral fit is performed with a binned likelihood method using the P8R3_SOURCE_V2 instrument response function with parameterization from 4FGL to characterize the spectral parameters of the source of interest in the 100 MeV–316 GeV energy range.

The gamma-ray light curve for each source is built with an end-to-end analysis in each time bin using the same processing steps as the spectral analysis. The period between MJD 54682 and MJD 58739 is divided into ∼570 7 days-long bins. A model spectral fit ¹² is performed in each time bin with the free model parameters. The best-fit parameters are then used to calculate the gamma-ray flux in the 0.1–316 GeV energy range for each time bin of the light curve.

Figure 2 shows examples of gamma-ray light curves for four blazars, including the first plausible cosmic neutrino source, TXS 0506+056 . In all light curves, we show flux points with their 1σ uncertainties (black symbols), if the source is detected with a test statistic larger than 9 (corresponding to a 3σ excess) and the flux in one bin is larger than its uncertainty. Otherwise, we show the flux upper limits at the 95% confidence level.

Zoom In Zoom Out Reset image size

3.1.1. Quiescent Flux Level

3.1.2. Flaring Threshold Level

Various definitions of flaring states in blazar emission exist in the literature (e.g., Resconi et al. 2009; Nalewajko 2013; Meyer et al. 2019; Stathopoulos et al. 2022). Here, we define, as a gamma-ray flare, any flux point of the weekly binned light curve that exceeds a certain threshold, ${F}_{\gamma }^{\mathrm{th}}$, given by

$$\begin{eqnarray}&&{F}_{\gamma }^{\mathrm{th}}={F}_{\gamma }^{q}+s\langle {F}_{\gamma }^{\mathrm{err}}\rangle ,\end{eqnarray} \tag{ 1 }$$

where ${F}_{\gamma }^{q}$ is the gamma-ray quiescent flux level, $\langle {F}_{\gamma }^{\mathrm{err}}\rangle$ is the average uncertainty of the gamma-ray fluxes, ¹³ and s corresponds to the significance above the quiescent flux level in units of the standard deviation σ. Unless otherwise noted, we use s = 6 in this work. The flaring threshold levels are indicated with dotted red lines in Figure 2.

The fraction of time spent in the flaring state, i.e., the flare duty cycle, is defined as

$$\begin{eqnarray}&&{f}_{\mathrm{fl}}=\displaystyle \frac{1}{{T}_{\mathrm{tot}}}{\int }_{{F}_{\gamma }^{\mathrm{th}}}{{dF}}_{\gamma }\displaystyle \frac{{dT}}{{{dF}}_{\gamma }},\end{eqnarray} \tag{ 2 }$$

where T_tot is the total observation time, F_γ is the gamma-ray energy flux, and T is the time spent at the respective flux level. The fraction of energy emitted in the flaring state, i.e., the gamma-ray flare energy fraction, is given by Murase et al. (2018):

$$\begin{eqnarray}&&{b}_{\mathrm{fl}}^{\gamma }=\displaystyle \frac{1}{{F}_{\gamma }^{\mathrm{ave}}{T}_{\mathrm{tot}}}{\int }_{{F}_{\gamma }^{\mathrm{th}}}{{dF}}_{\gamma }{F}_{\gamma }\displaystyle \frac{{dT}}{{{dF}}_{\gamma }}\end{eqnarray} \tag{ 3 }$$

where ${F}_{\gamma }^{\mathrm{ave}}$ is the average gamma-ray energy flux over the whole observation period. Note that Murase et al. (2018) introduced the flare duty cycle and energy fraction with gamma-ray luminosity L_γ .

Figure 3 shows the cumulative distributions of the 0.1–316 GeV gamma-ray flare duty cycles (f_fl) and flare energy fractions (${b}_{\mathrm{fl}}^{\gamma }$) of 145 blazars (106 FSRQs, 31 BL Lac objects, and 8 BCUs). The f_fl and ${b}_{\mathrm{fl}}^{\gamma }$ values for TXS 0506+056 are indicated with vertical dashed lines. TXS 0506+056 belongs to the 30% (40%) of the cumulative distribution of the flare duty cycle (flare energy fraction). Thus, TXS 0506+056 is not an extraordinary source in terms of its gamma-ray flaring activity. We also perform a two-sample Kolmogorov–Smirnov (K-S) test to check if the cumulative distributions of either f_fl or ${b}_{\mathrm{fl}}^{\gamma }$ are the same for FSRQs and BL Lacs. Under the null hypothesis that the FSRQ and BL Lac distributions are derived from the same parent distribution, the probabilities are p = 0.039 for the flare duty cycles distributions with the K-S statistic D = 0.279, and p = 0.521 for the gamma-ray flare energy fraction distributions, with D = 0.159, respectively. Hence, we can reject the null hypothesis that the distributions of flare duty cycles for FSRQs and BL Lac objects are drawn from the same parent population at a significance level of 5%. On the other hand, we cannot exclude the null hypothesis for the flare energy fractions of FSRQs and BL Lac objects at a significance level of 10%.

Zoom In Zoom Out Reset image size

Figure 4 shows a strong correlation between the flare duty cycles (f_fl) and the flare energy fractions (${b}_{\mathrm{fl}}^{\gamma }$) of 141 blazars (103 FSRQs, 30 BL Lac objects, and 8 BCUs); here, the nonflaring sources are excluded at the 6σ threshold. We find an almost linear relation between the two quantities for f_fl ≲ 0.1, which however becomes sublinear for higher flaring duty cycles. The side panels of Figure 4 show histograms of the logarithm of the flare duty cycle and energy fraction. It can be shown that both distributions are well represented by power-law functions, namely ${dN}/{{df}}_{\mathrm{fl}}\propto {f}_{\mathrm{fl}}^{-{\beta }_{f}}$ for f_fl ≳ 0.01 and ${dN}/{{db}}_{\mathrm{fl}}\propto {({b}_{\mathrm{fl}}^{\gamma })}^{-{\beta }_{b}}$ for ${b}_{\mathrm{fl}}^{\gamma }\gtrsim 0.03$, with β_f ∼ 1.7 and β_b ∼ 1.4. In other words, the blazars with lower flare duty cycles and flare energy fractions are more common in the sample compared to sources that are commonly flaring in gamma rays. Figure 3 also depicts the cumulative power-law functions with β_f = 1.7 for 0.006–1.0 flare duty cycles, and β_b = 1.4 for 0.02–1.0 flare energy fractions.

Zoom In Zoom Out Reset image size

So far, we have presented results for f_fl and ${b}_{\mathrm{fl}}^{\gamma }$ obtained using the 6σ threshold for the definition of flares. To check how our results depend on the choice for s, we construct the cumulative histograms of both quantities for integer values of s ranging from 2 to 9. Our results are shown in Figure 5. By lowering the threshold for the flaring flux, both f_fl and ${b}_{\mathrm{fl}}^{\gamma }$ distributions shift to higher values as expected. The duty cycle distribution tends to saturate only for s ≥ 8. Still, for s ≥ 6, we find small differences in the cumulative distributions of the duty cycle for f_fl ≳ 0.05. This suggests that our default choice of s = 6 is sufficient for describing the high-end of the duty cycle distribution. The impact on the flare energy fraction distribution is smaller, as indicated by the smaller spread of the curves. For completeness, we also show the cumulative histograms for s = 6 and two choices of the false alarm probability (p₀ = 0.01 and 0.10) of the Bayesian block algorithm that we applied to the Fermi-LAT light curves (see Section 3.1.1). The distributions are not sensitive to the choice of p₀ at least for the values we consider.

Zoom In Zoom Out Reset image size

Having determined the flare duty cycle for different blazar subclasses and having discussed the effects of the flare threshold value, we can now check if there is a correlation between the flare duty cycle and gamma-ray luminosity. In other words, we would like to check if the (quiescent) gamma-ray luminosity is a good indicator of the flaring activity of a blazar. We compute the quiescent gamma-ray luminosity as ${L}_{\gamma }^{q}=4\pi {d}_{L}^{2}{F}_{\gamma }^{q}$ without the correction for extragalactic background light absorption, where d_L is the luminosity distance calculated with the redshift z, and ${F}_{\gamma }^{q}$ is the 0.1–316 GeV gamma-ray quiescent flux (see Section 3.1.1). According to the 4FGL catalog, the blazars for z > 1.0 show spectral turnovers less than 100 GeV, indicating that the extragalactic background light absorption is irrelevant. As an uncertainty in the flare duty cycle, we use the square root of the sum of squares of the flaring time-bin errors, which are half of the time-bin width, over the total observation time. As shown in Figure 6, there is no correlation between the gamma-ray flare duty cycle and the quiescent gamma-ray luminosity for neither blazar subclass. Similar results are found for different flare threshold values s = 2–9.

Zoom In Zoom Out Reset image size

One of the ways to scale the neutrino luminosity is through the gamma-ray luminosity. Such a scaling is reasonable in hadronic models for gamma-ray emission, but it is nontrivial in leptonic models. Nevertheless, it may approximately hold given that the Compton dominance parameter is similar. The quiescent muon neutrino flux is benchmarked with the quiescent gamma-ray flux, and the flaring muon neutrino flux is written as

$$\begin{eqnarray}&&{E}_{{\nu }_{\mu }}{F}_{{E}_{{\nu }_{\mu }}}^{\mathrm{fl}}={E}_{{\nu }_{\mu }}{F}_{{E}_{{\nu }_{\mu }}}^{q}{\left(\displaystyle \frac{{F}_{\gamma }^{\mathrm{fl}}}{{F}_{\gamma }^{q}}\right)}^{\gamma }={A}_{\gamma }\,\displaystyle \frac{{E}_{\gamma }{F}_{{E}_{\gamma }}^{q}}{3}{\left(\displaystyle \frac{{F}_{\gamma }^{\mathrm{fl}}}{{F}_{\gamma }^{q}}\right)}^{\gamma },\end{eqnarray} \tag{ 6 }$$

where A_γ is a normalization parameter, and the factor of 3 comes from neutrino flavor mixing in vacuum ν_e : ν_μ : ν_τ ≃ 1: 1: 1. Here, we relate the differential quiescent neutrino flux at 300 TeV to the quiescent gamma-ray flux at the pivot energies for each weekly time bin. The pivot energies, which are chosen close to the decorrelation energy to minimize the correlation between the fitted gamma-ray spectral parameters, range from 0.1 to 316 GeV even in one source. The typical values are, however, 0.2–0.3 GeV. Then, the quiescent gamma-ray flux ${E}_{\gamma }{F}_{{E}_{\gamma }}^{q}$ is estimated from the ${E}_{\gamma }{F}_{{E}_{\gamma }}$ light curve at the pivot energies in the same way as the quiescent gamma-ray flux—see Section 3.1.1.

In photomeson production interactions in the blazar jet, the neutrinos and other secondaries produced via the neutral and charged pion decay chains carry comparable energy fluxes (see, e.g., Equations (13) and (14) of Murase et al. 2018). As a result, the neutrino energy flux is accompanied by a comparable energy flux of electromagnetic cascade emission at lower energies. However, not all X-rays and gamma rays are produced in photomeson production interactions. They may as well originate from Bethe–Heitler pair production and leptonic processes that are not accompanied by neutrinos. Therefore, we adopt A_γ = 1 as an upper limit for the normalization parameter. A stronger upper limit on A_γ can be obtained from the nondetection of neutrinos from some of the blazars in our sample as illustrated in the following figures.

Figure 8 shows the estimated muon neutrino flare fluxes from scenario (1) with A_γ = 1.0, and γ = 1.5 as a function of $\sin (\delta )$, where δ is the source decl., neglecting the muon neutrino flux equal to or less than its uncertainty. The results for 1 week and 10 yr bins are shown in the left and right panels, respectively. ¹⁴ The neutrino flare fluxes for a 10 yr bin are the averages of the weekly neutrino flare fluxes, with zero fluxes assumed for the nonflaring states. For both binning choices, there are several blazar flares with predicted neutrino fluxes that exceed the IceCube sensitivity. The lack of neutrino detections from these sources can therefore be used to derive limits on the normalization A_γ for different values of the index γ (gray lines), as demonstrated in Figure 9. As expected, the weaker the relation between neutrino and gamma-ray energy fluxes is (i.e., lower γ values), the higher A_γ can be. The red, orange, brown, and purple lines indicate stringent limits obtained for a couple of sources in our sample, such as 3C 454.3, 3C 279, PKS 1502+106, and PKS 1510-089. The estimated neutrino flare fluxes with a 10 yr bin are lower than those with a 1 week bin because of the averaging. Meanwhile, the IceCube 90% sensitivity for a 10 yr bin is much lower than that for 1 week bin. Hence, the limits placed on A_γ from the 10 yr binning analysis are more stringent than those set by the weekly binned analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This case is motivated by modeling studies of the TXS 0506+056 multimessenger observations of TXS 0506+056 (e.g., Keivani et al. 2018; Petropoulou et al. 2020), which demonstrated that the neutrino fluxes are bound by the X-ray data. The quiescent muon neutrino flux is benchmarked with the quiescent X-ray flux, ${E}_{X}{F}_{{E}_{X}}^{q}$, and the muon neutrino flare flux is written as

$$\begin{eqnarray}&&{E}_{{\nu }_{\mu }}{F}_{{\nu }_{\mu }}^{\mathrm{fl}}={E}_{{\nu }_{\mu }}{F}_{{E}_{{\nu }_{\mu }}}^{q}{\left(\displaystyle \frac{{F}_{\gamma }^{\mathrm{fl}}}{{F}_{\gamma }^{q}}\right)}^{\gamma }={A}_{X}\,\displaystyle \frac{{E}_{X}{F}_{{E}_{X}}^{q}}{3}{\left(\displaystyle \frac{{F}_{\gamma }^{\mathrm{fl}}}{{F}_{\gamma }^{q}}\right)}^{\gamma }\end{eqnarray} \tag{ 7 }$$

where A_X is a normalization parameter. For γ = 1.0–2.0, we can set an upper limit on A_X so that our predictions are consistent with the nondetection of neutrinos from flares of individual blazars.

For the X-ray data, we use the Open Universe for Blazars, which provides blazar X-ray light curves based on 14 yr of Swift-XRT data (Giommi et al. 2019). The 0.5, 1.0, 1.5, 3.0, and 4.5 keV quiescent fluxes are estimated in the same way as the gamma-ray quiescent level (see Section 3.1.1). The X-ray quiescent flux, ${E}_{X}{F}_{{E}_{X}}^{q}$, is then defined as the average of the quiescent fluxes in these five energies. Figure D1 in Appendix D shows examples of the 1.0 keV X-ray light curves with the respective Bayesian blocks' representations and quiescent fluxes. By using Equation (7), we derive the neutrino flare fluxes of 129 blazars (97 FSRQs, 26 BL Lac objects, 6 BCUs), ignoring 16 sources without flaring states found at the 6σ threshold or without well-defined X-ray quiescent state.

In Figure 10, we plot the estimated muon neutrino flare fluxes from scenario (2) with A_X = 1.0, and γ = 1.5 as a function of $\sin (\delta )$, for two binning choices (panels (a) and (b)). The estimated neutrino flare fluxes in scenario (2) are generally lower by 1 or 2 orders of magnitude than those in scenario (1) (compare to Figure 8). This result is related to the different spectral energy distributions (SEDs) of blazars in our sample. The majority of them are FSRQs, which are typically low-synchrotron peaked sources, and their SED has a trough in the X-ray band covered by Swift-XRT. As a result, the X-ray fluxes used as a proxy to scale the neutrino luminosity are lower than the gamma-ray fluxes used in scenario (1) (see also Figure 7). This does not apply, however, to MeV blazars and extreme synchrotron peaked sources (Costamante 2020). Figure 11 shows the normalization parameter A_X in gray as a function of the indices γ based on the nondetection of neutrino events, compared to the IceCube 90% sensitivity for 1 week binning (left) and 10 yr binning (right). The red, blue, green, and orange lines show sources with the lowest normalization parameters (see inset legend).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also investigate the sensitivity of our results to a non-power-law shape of the neutrino spectrum. To this end, we assume a more realistic, albeit model-dependent neutrino spectrum motivated by leptohadronic models of TXS 0506+056. We compare the expected number of neutrinos to the one expected under the assumption of an ${E}_{\nu }^{-2}$ neutrino spectrum, which we have used in the previous sections. We choose the LMBB2b model of Keivani et al. (2018) as an illustrative example of a characteristic neutrino spectrum that might result from photomeson production processes in the jet of a blazar during a gamma-ray flare.

We use Equations (6) and (7) with a 10 yr binning to test this variant of scenarios (1) and (2) respectively, and to place upper limits on the normalization parameters, A_γ and A_X . To calculate the expected number of neutrinos, we assume the IC86 point-source effective area given in Aartsen et al. (2017b).

Columns (4)–(5) and (6)–(7) of Table 1 can also be used to see the effect of the more realistic spectral shape. For TXS 0506+056 and 3C 454.3, the expected number of neutrinos with the more realistic model is a factor of 2 lower than that of an ${E}_{\nu }^{-2}$ neutrino spectrum. Therefore, the respective upper limits on A_γ and A_X  of the more realistic model are a factor of 2 less stringent constraints. We thus conclude that the results obtained with our standard assumption of an ${E}_{\nu }^{-2}$ spectral shape may be up to a factor of a few more strict for A_γ and A_X than with a different spectral shape for a given total neutrino energy flux. Regardless of the neutrino spectral template, the upper limits on ${A}_{\gamma /X}^{\mathrm{up}}$ are several orders of magnitude more constraining for 3C 454.3 than for TXS 0506+056.

Table 1. Neutrino Counts Using the LMBB2b Model of Keivani et al. (2018; Columns 4, 5) and an ${E}_{\nu }^{-2}$ Neutrino Spectrum (Columns 6, 7)

			LMBB2b		${E}_{\nu }^{-2}$
Source	Scenario	γ	${N}_{{\nu }_{\mu }+{\bar{\nu }}_{\mu }}({A}_{\gamma /X}=1$)	${A}_{\gamma /X}^{\mathrm{up}}$ (90% CL)	${N}_{{\nu }_{\mu }+{\bar{\nu }}_{\mu }}({A}_{\gamma /X}=1$)	${A}_{\gamma /X}^{\mathrm{up}}$ (90% CL)
TXS 0506+056	1	1.0	2	1	4	1
TXS 0506+056	1	1.5	6	1	10	1
TXS 0506+056	1	2.0	10	1	20	1
TXS 0506+056	2	1.0	0.3	1	0.5	1
TXS 0506+056	2	1.5	0.9	1	1.0	1
TXS 0506+056	2	2.0	2	1	3	1
3C 454.3	1	1.0	20	0.1	70	0.04
3C 454.3	1	1.5	200	10−2	600	4 × 10−3
3C 454.3	1	2.0	2000	10−3	6000	4 × 10−4
3C 454.3	2	1.0	10	0.2	30	0.08
3C 454.3	2	1.5	100	2 × 10−2	300	8 × 10−3
3C 454.3	2	2.0	1000	2 × 10−3	3000	8 × 10−4

Note. The assumed upper limit is 2.44 neutrinos in the livetime of IceCube. The IC86 point-source effective area has been assumed.

Download table as:  ASCIITypeset image

There are two ways to examine the contribution of blazars to the all-sky neutrino flux. One is the stacking method considering the directions of IceCube neutrino events (Aartsen et al. 2017a; Hooper et al. 2019; Abbasi et al. 2022), and the other is the clustering analysis (Murase & Waxman 2016). Instead, our flux stacking method uses the estimated muon neutrino fluxes, not directly using the spatial correlation between the directions of IceCube neutrino events and blazars, but using the IceCube sensitivity at the decl. of each blazar in our sample.

As shown in Figures 8 and 10, some neutrino flare fluxes with A_γ = 1.0, and A_X = 1.0 exceed the IceCube 90% sensitivity. The normalization parameters of A_γ and A_X are set to 1.0 for the sources whose flares are all less than the sensitivity flux. However, the normalization parameter of the source is reduced, if the flux of at least one flare exceeds the sensitivity, in a way that the maximum neutrino flux of the source is equal to the sensitivity flux.

In this work, we estimate the upper limits on the diffuse isotropic neutrino flux through flux stacking analyses of blazar flares. Our sample consists of bright gamma-ray blazars from the 4LAC catalog as shown in Figure 1. There are fainter blazars in the 4LAC catalog, and there ought to be a lot of not-yet-detected blazars in the Universe.

To obtain the upper limits on the all-sky neutrino flux, the sum of the estimated neutrino fluxes from the flaring blazars in our sample is divided by 4π and then multiplied by two factors for the completeness correction. First, the correction factor of all 4LAC blazars and the blazars of our subsample is obtained by the fraction of the total gamma-ray energy fluxes in 0.1–100 GeV (raised to the power of γ) of 4LAC blazars to the sample blazars:

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{{N}_{c}}{(\mathrm{Energy}\ {\mathrm{Flux}}_{i})}^{\gamma }/\displaystyle \sum _{j=1}^{{N}_{s}}{(\mathrm{Energy}\ {\mathrm{Flux}}_{j})}^{\gamma }\end{eqnarray} \tag{ 8 }$$

where N_c is the number of blazars in the 4LAC catalog, N_s is the number of blazars in our sample, and Energy Flux_i,j are time-average gamma-ray energy fluxes in 0.1–100 GeV. The second correction factor for the blazar population, which includes too distant or too low in gamma-ray luminosity sources, and the 4LAC blazars, is obtained by completeness correction factors with γ, which are taken from Figure 1 of Yuan et al. (2020). Since the analysis of Yuan et al. (2020) originates from the 3LAC catalog, our results are slightly conservative when applying it to the 4LAC catalog. The correction factors are presented in Table 2.

Figure 12 shows the total fluxes of muon neutrinos estimated from our bright gamma-ray FSRQs and BL Lac objects, multiplied by correction factors of the contribution of the 4LAC blazars and still unresolved blazars that are too gamma-ray faint to be included in the 4LAC catalog. The total muon neutrino fluxes of FSRQs and BL Lac objects of scenarios (1) and (2) are plotted as a function of γ, compared to a reference flux of the isotropic diffuse muon neutrino flux (Aartsen et al. 2016) and the three upper limits derived by Hooper et al. (2019), which are reduced to 1/3 for muon neutrinos. As the neutrino energy fluxes are enhanced by the weighting factor γ, the sum of the neutrino energy fluxes of the sample blazars gradually increases with γ. On the other hand, the correction factors sharply decrease with γ up to γ ≃ 1.5 and stay nearly constant for γ ≳ 1.5 as shown in Table 2. Therefore, the total neutrino energy fluxes, namely the product of the sum of the neutrino energy fluxes of the sample blazars and the correction factors, decrease with γ up to γ ≃ 1.5 and turn to increase gradually with γ for γ ≳ 1.5 as shown in scenario (1) of Figure 12. On the other hand, the total neutrino energy fluxes are maintained nearly constant for γ ≳ 1.5 as shown in scenario (2) of Figure 12, since the gamma-ray normalization factors A_γ are strongly suppressed in tens of blazars shown in Figure 9 (right). The upper limits of scenario (1) are in agreement with those of Hooper et al. (2019) for γ ≳ 1.5 even though the analyses and assumptions are different. We predict weaker upper limits, though, for smaller values of γ. This is because the second correction factor should be larger for small values of γ (Yuan et al. 2020). We also note that the upper limits of scenario (2) are tighter than those of scenario (1) for all values of γ.

Zoom In Zoom Out Reset image size

So far, we have defined gamma-ray flares with the 6σ level (s = 6) above the quiescent gamma-ray flux. Figure 13 presents the upper-limit values of ${E}_{{\nu }_{\mu }}^{2}{{\rm{\Phi }}}_{{\nu }_{\mu }}$ as a function of the flare significance s above the quiescent level for FSRQs and BL Lac objects in both scenarios with γ = 1.5. The upper limits of FSRQs and BL Lac objects are similar both for scenarios (1) and (2) with each other. The upper limits decrease with the significance levels s, becoming ∼70% from s = 2–9 for FSRQs and BL Lac objects in scenario (1).

Zoom In Zoom Out Reset image size