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

Kenji Yoshida  Department of Electronic Information Systems, Shibaura
Institute of Technology, 307 Fukasaku, Minuma-ku, Saitama, Japan  Maria
Petropoulou  National and Kapodistrian University of Athens, 157 72, Athens,
Greece  Institute of Accelerating Systems & Applications, University Campus
Zografos, Athens, Greece  Kohta Murase  Department of Physics; Department of
Astronomy & Astrophysics; Center for Multimessenger Astrophysics, Institute
for Gravitation and the Cosmos, The Pennsylvania State University, University
Park, PA, USA  School of Natural Sciences, Institute for Advanced Study,
Princeton, NJ, USA  Center for Gravitational Physics and Quantum Information.
Yukawa Institute for Theoretical Physics, Kyoto, Kyoto, Japan  Foteini
Oikonomou  Physics Department, Norwegian University of Science and Technology,
7491, Trondheim, Norway

###### Abstract

Gamma-ray flares of blazars may be accompanied by high-energy neutrinos due to
interactions of high-energy cosmic rays in the jet with photons, as suggested
by the detection of the high-energy neutrino IceCube-170922A during a major
gamma-ray flare from blazar TXS 0506+056 at the $\sim 3\sigma$ significance
level. In this work, we present a statistical study of gamma-ray emission from
blazars to constrain the contribution of gamma-ray flares to their neutrino
output. We construct weekly binned light curves for 145 gamma-ray bright
blazars in the Fermi Large Area Telescope (LAT) Monitored Source List adding
TXS 0506+056. We derive the fraction of time spent in the flaring state (flare
duty cycle) and the fraction of energy released during each flare from the
light curves with a Bayesian blocks algorithm. We find that blazars with lower
flare duty cycles and energy fractions are more numerous among our sample. We
identify a significant difference in flare duty cycles between blazar sub-
classes at a significance level of 5 %. Then using a general scaling relation
for the neutrino and gamma-ray luminosities,
$L_{\nu}\propto(L_{\gamma})^{\gamma}$ with a weighting exponent of
${\gamma}=1.0-2.0$, normalized to the quiescent gamma-ray or X-ray flux of
each blazar, we evaluate the neutrino energy flux of each gamma-ray flare. The
gamma-ray flare distribution indicates that blazar neutrino emission may be
dominated by flares for $\gamma\gtrsim 1.5$. The neutrino energy fluxes for
one-week and 10-year bins are compared with the declination-dependent IceCube
sensitivity to constrain the standard neutrino emission models for gamma-ray
flares. Finally, we present the upper-limit contribution of blazar gamma-ray
flares to the isotropic diffuse neutrino flux.

High energy astrophysics — Gamma-ray astronomy — Neutrino astronomy — Blazars

##  1 Introduction

In 2013 the IceCube Neutrino Observatory discovered an all-sky flux of
astrophysical neutrinos with energies from $\sim 100$ TeV to several PeV
(Aartsen et al., 2013a, b). Since then, the spectrum of this high-energy
neutrino flux has been measured to a higher precision with different analyses
(e.g., Aartsen et al., 2015a, 2016a, 2020a). The nearly isotropic distribution
of the neutrino arrival directions constrains a Galactic contribution to the
measured IceCube flux to the level of $\sim 10$% (Albert et al., 2018)
supporting an extragalactic origin (e.g., Aartsen et al., 2014, 2015b).

Many proposals have been put forward to explain the diffuse neutrino flux as
the cumulative emission from a population of sources, with several of them
even predating its discovery (for recent reviews, see Ahlers & Halzen, 2017;
Mészáros, 2017). Blazars, namely active galactic nuclei (AGN) with
relativistic jets closely aligned to our line of sight (Urry & Padovani,
1995), belong to this long-list of candidate sources. Due to their powerful
magnetized jets and copious radiation fields, blazars are promising sites for
high-energy neutrino production (see reviews, e.g., Murase, 2017; Murase &
Stecker, 2022, and references therein).

In 2017 IceCube detected a muon neutrino event with most probable energy of
$\sim 290$ TeV, IceCube-170922A, in spatial and temporal coincidence with a
flaring gamma-ray blazar, TXS 0506+056, at the statistical significance level
of $3{\sigma}$ (IceCube Collaboration et al., 2018a). Prior to the 2017
flaring event, the IceCube Collaboration also reported $3.5{\sigma}$ evidence
for a neutrino excess on top of the atmospheric background expectation from
the direction of TXS 0506+056 with an archival search of past IceCube data
(IceCube Collaboration et al., 2018b). In addition, some blazars may be
powerful enough to produce a strong point-source neutrino signal during flares
and have been associated with high-energy neutrinos (Dermer et al., 2014;
Kadler et al., 2016; Petropoulou et al., 2016; Gao et al., 2017; Oikonomou et
al., 2019; Paliya et al., 2020; Petropoulou et al., 2020; Rodrigues et al.,
2021; Oikonomou et al., 2021; Sahakyan et al., 2022).

Clustering limits (Murase & Waxman, 2016; Capel et al., 2020) and stacking
limits (Aartsen et al., 2017; Neronov et al., 2017; Hooper et al., 2019; Smith
et al., 2021; Abbasi et al., 2022) suggest that blazars are unlikely to be the
dominant sources of the diffuse astrophysical neutrino flux detected by
IceCube (Murase et al., 2018; Palladino et al., 2019; Yuan et al., 2020).
However, they may still appear in stacking analyses of radio-bright blazar
arrival directions (Hovatta et al., 2021; Plavin et al., 2021; Zhou et al.,
2021; Buson et al., 2022), which may support that the dominant neutrino
sources are gamma-ray hidden or opaque (Murase et al., 2016) or the
correlation level may be explained by a handful of the brighest blazars
(Murase et al., 2014). The role of blazar flaring activity in neutrino-source
associations is not yet fully understood. Moreover, the results of most
stacking analysis are relevant to the time-average blazar emission. Murase et
al. (2018) estimated the contribution of blazar flares like the 2017 multi-
messenger flare of TXS 0506+056 to the diffuse neutrino flux to be up to $\sim
10\%$ based on the Fermi All-sky Variability Analysis (FAVA), but the
constraint on the contribution from all blazar flares is very uncertain
because it is affected by different assumptions (Yuan et al., 2020).

In this paper we expand upon the work of Murase et al. (2018) by computing the
duty cycle and the energy output fraction of gamma-ray ray flares using a much
larger sample of blazars observed by the Large Area Telescope (LAT) on board
the _Fermi_ satellite. To estimate the neutrino flux during gamma-ray flares
we use a generic relationship between the all-flavor neutrino and gamma-ray
fluxes, i.e. $F^{\rm fl}_{\nu}\propto\left(F_{\gamma}^{\rm
fl}\right)^{\gamma}$, where $1\leq\gamma\leq 2$ (see also Yuan et al., 2020).
We explore two cases depending on the definition of the neutrino flux in non-
flaring periods (i.e. quiescence). First, we study the case where the neutrino
flux outside flares is benchmarked with the quiescent gamma-ray flux, which
has been widely adopted prior to the TXS 0506+056 observations (e.g., Murase
et al., 2014; Petropoulou et al., 2015; Padovani et al., 2015). Second,
motivated by the modeling studies of TXS 0506+056 we consider a case where the
non-flaring neutrino flux is bound by the quiescent X-ray flux. Besides the
usual $E_{\nu}^{-2}$ power-law neutrino spectrum, we consider for both
scenarios the effects of realistic neutrino spectra, using a spectral template
derived for the multi-messenger flare of TXS 0506+056 . By comparing our
neutrino predictions for each blazar in the sample to the time-integrated
10-year IceCube sensitivity and the IceCube sensitivity for transient searches
(scaled to one-week) we constrain the aforementioned scenarios. We finally
compute the contribution of blazar flares to the all-sky diffuse neutrino
flux.

This paper is structured as follows. In Sec. 2 we describe the sample and the
_Fermi_ -LAT analysis for generating the light curves. We describe next how we
define gamma-ray flaring states using the binned LAT light curves and present
our results about the duty cycle of flares and their energy output fraction in
Sec. 3. We discuss the implications of our results on the high-energy neutrino
emission from flaring blazars in Sec. 4. We conclude in Sec. 5 with a summary
of our results. In this work we adopt a flat $\Lambda$CDM cosmology with
$H_{0}=70$ km s-1 Mpc-1 and $\Omega_{\rm m}=0.30$.

##  2 Gamma-ray data and analysis

###  2.1 Data selection

We select 145 gamma-ray bright blazars (Table LABEL:tab:source_list),
consisting of 106 flat spectrum radio quasars (FSRQs), 31 BL Lac objects (BL
Lacs), and 8 blazar candidates of uncertain type (BCUs). Among them, the 144
blazars are selected from the _Fermi_ Large Area Telescope (LAT) Monitored
Source List111https://fermi.gsfc.nasa.gov/ssc/data/access/lat/msl_lc/, 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 multi-wavelength observations of
flaring sources. When sources exceed the monitoring flux threshold of
$1{\times}10^{-6}$ cm-2s-1, they are added to the list of monitored sources
222https://fermi.gsfc.nasa.gov/ssc/data/access/lat/msl_lc/. 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 years 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.

Figure 1: Scatter diagram of the time-average gamma-ray energy flux and
redshift for the 145 blazars in our sample (106 FSRQs in blue, 31 BL Lacs in
orange, and 8 BCUs in green). All blazars in the 4LAC are overlaid (grey
markers) for comparison reasons. Side panels show the projected histograms on
each axis. The energy fluxes are taken from 4LAC and the redshifts are from
4LAC and SIMBAD (Oberto et al., 2019).

###  2.2 Fermi LAT analysis

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-03333http://fermi.gsfc.nasa.gov/ssc/data/analysis/software and
Fermipy version 0.17.3 (Wood et al., 2017). Following the standard selection
recommendation444https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data_Exploration/Data_preparation.html,
we select photons with energies between 100 MeV and 316 GeV that were detected
within a radius of $20^{\circ}$ 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^{\circ}$ 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^{\circ}$ 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
parametrization 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 $\sim$570 7 day-long
bins. A model spectral fit555The spectral models used are from the 4FGL, e.g.,
PowerLaw, LogParabola, and PLSuperExpCutoff. In LogParabola the free
parameters are norms, $\alpha$, and $\beta$. 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\sigma$ uncertainties (black
symbols), if the source is detected with a test statistic (TS) larger than 9
(corresponding to a $3\sigma$ excess) and the flux in one bin is larger than
its uncertainty. Otherwise, we show the flux upper limits at the 95 %
confidence level.

##  3 Gamma-Ray flare distributions

In this section, we start by defining gamma-ray flaring states using the
weekly-binned _Fermi_ LAT light curves (Section 3.1). We then derive the flare
duty cycle (i.e., fraction of time spent in flaring states) and the flare
energy fraction (i.e., fraction of energy released in flaring states) for all
sources in our sample, and discuss potential differences between BL Lac
objects and FSRQs (Section 3.2). A comparison of our results with those
reported by Murase et al. (2018) can be found in Appendix E.

Figure 2:  Indicative examples of weekly-binned blazar gamma-ray light curves
(black points) in the 0.1-316 GeV energy range. Grey triangles indicate 95%
upper limits. The quiescent flux and flare threshold levels are plotted with
dashed blue and red lines, respectively (for a definition of both flux
thresholds, see Sec. 3.1). The Bayesian blocks representation of each light
curve is overplotted with orange solid lines. The 1-316 GeV light curves of
the same sources can be found in Appendix B.

###  3.1 Definition of flares

####  3.1.1 Quiescent flux level

Although there is no rigorous way to define the quiescent flux level of blazar
emission, it is reasonable to assume that the quiescent state has a relatively
low and stable flux. In this study, we adopt the following procedure to derive
the quiescent gamma-ray flux level of a blazar:

  1. 1.

We apply a Bayesian blocks algorithm to each gamma-ray light curve with a
false alarm probability $p_{0}=0.05$, adopting the so-called point
measurements fitness function for the Bayesian blocks algorithm (Scargle et
al., 2013; Astropy Collaboration et al., 2013, 2018).

  2. 2.

The Bayesian blocks algorithm provides the partitions of the light curve into
blocks where the flux is regarded to be constant. This is an optimal step-
function representation of the light curve.

  3. 3.

We consider the minimum average flux as a representative of the source
quiescent flux, under the following condition: the number of data points
within the selected block is more than or equal to a half of the average
number of data points of all other blocks. This condition prevents us from
using an insufficient amount of data points to estimate the quiescent flux
level. We note that the quiescent flux levels hardly depend on the choice of
$p_{0}$ (in the range $0.01-0.1$) as described in Sec. 3.1.2.

Figure 2 shows indicative examples of weekly-binned _Fermi_ LAT light curves
from our sample (black points) and the Bayesian blocks representations (solid
orange lines). The quiescent flux levels derived using the method described
above are also indicated with dashed cyan lines. To verify the quiescent flux
level being identified properly, we derived the gamma-ray quiescent fluxes by
using coarser time binning for the light curves. See Appendix C for details.

Figure 3:  Cumulative histograms (normalized to one) of the $0.1-316$ GeV
gamma-ray flare duty cycles (a) and flare energy fractions (b) of 145 blazars
(black solid histogram), 106 FSRQs (blue solid histogram), and 31 BL Lacs
(orange solid histogram). In both panels, the vertical dashed line indicates
the values for TXS 0506+056 .The dashed gray lines depict cumulative power-law
functions (see text for more details).

####  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}^{\rm th}$, given by

$F_{\gamma}^{\rm th}=F_{\gamma}^{q}+s{\langle}F_{\gamma}^{\rm err}{\rangle},$
(1)

where $F_{\gamma}^{q}$ is the gamma-ray quiescent flux level,
${\langle}F_{\gamma}^{\rm err}{\rangle}$ is the average uncertainty of the
gamma-ray fluxes666${\langle}F_{\gamma}^{\rm
err}{\rangle}=(\sum_{i=1}^{N}\sigma_{i})/N$, where $N$ is the number of time
bins of a light curve and $\sigma_{i}$ is the flux uncertainty of the $i$-th
time bin., and $s$ corresponds to the significance above the quiescent flux
level in units of the standard deviation ${\sigma}$. Unless otherwise noted,
we use $s=6$ in this work. The flaring threshold levels are indicated with
dotted red lines in Figure 2.

###  3.2 Flare duty cycle and flare energy fraction

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

$f_{\rm fl}=\frac{1}{T_{\rm tot}}\int_{F_{\gamma}^{\rm
th}}dF_{\gamma}\frac{dT}{dF_{\gamma}}$ (2)

where $T_{\rm tot}$ is the total observation time, $F_{\gamma}$ 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)

$b^{\gamma}_{\rm fl}=\frac{1}{F_{\gamma}^{\rm ave}T_{\rm
tot}}\int_{F_{\gamma}^{\rm th}}dF_{\gamma}F_{\gamma}\frac{dT}{dF_{\gamma}}$
(3)

where $F_{\gamma}^{\rm 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_{\gamma}$.

Figure 4:  Scatter plot of the $0.1-316$ GeV gamma-ray flare duty cycles and
flare energy fractions of 141 blazars (103 FSRQs in blue squares, 30 BL Lacs
in orange circles, and 8 BCUs in green triangles). Side panels show the
projected histograms on each axis.

Figure 3 shows the cumulative distributions of the $0.1-316$ GeV gamma-ray
flare duty cycles ($f_{\rm fl}$) and flare energy fractions ($b^{\gamma}_{\rm
fl}$) of 145 blazars (106 FSRQs, 31 BL Lacs, and 8 BCUs). The $f_{\rm fl}$ and
$b^{\gamma}_{\rm fl}$ 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 test to check if the
cumulative distributions of either $f_{\rm fl}$ or $b^{\gamma}_{\rm fl}$ 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 Lacs
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 Lacs at a significance level of 10 %.

Figure 4 shows a strong correlation between the flare duty cycles ($f_{\rm
fl}$) and the flare energy fractions ($b^{\gamma}_{\rm fl}$) of 141 blazars
(103 FSRQs, 30 BL Lacs, and 8 BCUs); here, non-flaring sources are excluded at
the $6\sigma$ threshold. We find an almost linear relation between the two
quantities for $f_{\rm fl}\lesssim 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_{\rm fl}\propto f_{\rm fl}^{-\beta_{f}}$ for $f_{\rm fl}\gtrsim 0.01$
and $dN/db_{\rm fl}\propto(b_{\rm fl}^{\gamma})^{-\beta_{b}}$ for
$b^{\gamma}_{\rm fl}\gtrsim 0.03$, with $\beta_{f}~{}\sim 1.7$ and
$\beta_{b}~{}\sim 1.4$. In other words, 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 ${\beta}_{f}=1.7$ for $0.006-1.0$ flare duty cycles
and ${\beta}_{b}=1.4$ for $0.02-1.0$ flare energy fractions.

So far, we have presented results for $f_{\rm fl}$ and $b^{\gamma}_{\rm fl}$
obtained using the $6\sigma$ 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_{\rm fl}$ and $b^{\gamma}_{\rm fl}$ distributions, shift to
higher values as expected. The duty cycle distribution tend to saturate only
for $s\geq 8$. Still, for $s\geq 6$, we find small differences in the
cumulative distributions of the duty cycle for $f_{\rm fl}\gtrsim 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}=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_{0}$
at least for the values we consider.

Figure 5:  Cumulative histograms of the $0.1-316$ GeV gamma-ray flare duty
cycles (a) and flare energy fractions (b) of 145 blazars for different integer
values of $s$ ranging from 2 to 9. In the case of $s=6$ cumulative histograms
of $p_{0}=0.01,0.10$ are also plotted.

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$.

Figure 6: Scatter plot between $0.1-316$ GeV gamma-ray flare duty cycle and
quiescent gamma-ray luminosity for 137 blazars of 103 FSRQs (blue squares), 30
BL Lacs (orange circles), and 4 BCUs (green triangles).

##  4 Implications for high-energy neutrino emission

In order to discuss implications for neutrino emission, we need to know the
physical relationship between neutrinos and gamma rays. In $pp$ scenarios for
either gamma-ray or neutrino data, neutrinos and gamma rays can be simply
connected via the hadronic production process itself (Murase et al., 2013).
However, blazar gamma-ray emission is usually explained by primary leptonic
emission (e.g., inverse Compton scattering) and/or secondary cascade emission
from photohadronic ($p\gamma$) interactions, thus making the multimessenger
connection much more model dependent. It is therefore useful to parameterize
the multimessenger relationship in a generic manner using the luminosity
weighting factor $\gamma$ (e.g., Yuan et al., 2020),

$\displaystyle L^{\rm fl}_{\nu}=L^{q}_{\nu}\left(\frac{L^{\rm
fl}_{\gamma}}{L^{q}_{\gamma}}\right)^{\gamma},$ (4)

where $L^{\rm fl}_{\nu}$ is the flare neutrino luminosity, $L^{q}_{\nu}$ is
the quiescent neutrino luminosity, and $L^{\rm fl}_{\gamma}$ is the gamma-ray
luminosity in the flaring state. Theoretical models typically predict
$\gamma\sim 1.0-2.0$ (e.g., Murase et al., 2018). Note that the differential
gamma-ray energy flux, $F_{E_{\gamma}}$, is related to the bolometric gamma-
ray luminosity as

$\displaystyle L_{\gamma}=4\pi d_{L}^{2}\int dE_{\gamma}\,F_{E_{\gamma}}.$ (5)

There are different ways to estimate the non-flaring (quiescent) neutrino
flux, $F^{q}_{\nu}$. In what follows, we explore two cases where $F^{q}_{\nu}$
is benchmarked with the quiescent gamma-ray flux (scenario 1) or quiescent
X-ray flux (scenario 2). For each scenario, we consider the impacts of the
neutrino spectrum. For illustration purposes, we present in Figure 7 a scatter
plot of the quiescent gamma-ray and X-ray fluxes for 138 blazars of 101 FSRQs,
31 BL Lacs, and 6 BCUs. Note that the majority of sources has
$F_{\gamma}^{q}>F_{\rm X}^{q}$, which relates also to their broadband spectral
types.

Figure 7: Scatter plot between the quiescent gamma-ray and X-ray fluxes for
138 blazars of 101 FSRQs (blue squares), 31 BL Lacs (orange circles), and 6
BCUs (green triangles). The isolated upper left FSRQ is 3C 273.

###  4.1 Gamma-ray proxy (scenario 1)

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 non-trivial 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

$\displaystyle E_{{\nu}_{\mu}}F^{\rm
fl}_{E_{{\nu}_{\mu}}}=E_{{\nu}_{\mu}}F^{q}_{E_{{\nu}_{\mu}}}\left(\frac{F^{\rm
fl}_{\gamma}}{F^{q}_{\gamma}}\right)^{\gamma}=A_{\gamma}\,\frac{E_{\gamma}F^{q}_{E_{\gamma}}}{3}\left(\frac{F^{\rm
fl}_{\gamma}}{F^{q}_{\gamma}}\right)^{\gamma},$ (6)

where $A_{\gamma}$ is a normalization parameter, and the factor of 3 comes
from neutrino flavour mixing in vacuum
${\nu}_{e}:{\nu}_{\mu}:{\nu}_{\tau}~{}{\simeq}~{}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 weakly 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 GeV 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^{q}_{E_{\gamma}}$ 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, 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_{\gamma}=1$ as an upper limit for the normalisation parameter. A
stronger upper limit on $A_{\gamma}$ can be obtained from the non-detection 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_{\gamma}=1.0$ and ${\gamma}=1.5$ as a function of ${\sin}({\delta})$, where
$\delta$ is the source declination, neglecting the muon neutrino flux equal to
or less than its uncertainty. Results for one-week and 10-year bins are shown
in the left and right panels, respectively777In Appendix F, we also present
the estimated muon neutrino flare fluxes using the gamma-ray photon fluxes
instead of the gamma-ray energy fluxes.. The neutrino flare fluxes for a
10-year bin are the averages of the weekly neutrino flare fluxes, with zero
fluxes assumed for the non-flaring 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_{\gamma}$ for different
values of the index $\gamma$ (gray lines), as demonstrated in Figure 9. As
expected, the weaker the relation between neutrino and gamma-ray energy fluxes
is (i.e., lower $\gamma$ values), the higher $A_{\gamma}$ 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-year bin are lower
than those with a one-week bin because of the averaging. Meanwhile, the
IceCube 90% sensitivity for a 10-year bin is much lower than that for one-week
bin. Hence, the limits placed on $A_{\gamma}$ from the 10-year binning
analysis are more stringent than those set by the weekly binned analysis.

Figure 8: Estimated muon neutrino flare fluxes of 136 blazars (102 FSRQs, 26
BL Lacs, and 8 BCUs) for a one-week bin (a) and a 10-year bin (b) of scenario
1 with $A_{\gamma}=1.0$ and ${\gamma}=1.5$ as a function of ${\sin}({\delta})$
where ${\delta}$ is the source declination. The IceCube 90 % sensitivities
(red solid lines) for a point source with an $E_{\nu}^{-2}$ neutrino spectrum
for one week and 10 years are taken from Vandenbroucke et al. (2019) and
Aartsen et al. (2020b), respectively. Figure 9: Upper limits on the
normalization parameter $A_{\gamma}$ of the analyzed blazars in gray as a
function of the indices $\gamma$ based on the non-detection of neutrino
events, compared to the IceCube 90 % sensitivity for scenario 1. The left and
right panels show the normalization parameters with one-week binning and
10-year binning, respectively. When sources whose flares are all less than the
sensitivity flux, the normalization parameters of $A_{\gamma}$ are set to 1.0.
The normalization parameter of the source is adjusted so that the maximum
neutrino flux of the source is equal to the sensitivity flux if the flux of at
least one flare exceeds the sensitivity. Several sources with the lowest
normalization parameters are plotted with red, blue, brown, and purple lines.

###  4.2 X-ray proxy (scenario 2)

This case is motivated by modeling studies of the TXS 0506+056 multi-messenger
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^{q}_{E_{X}}$, and the muon neutrino flare flux is written as

$\displaystyle E_{{\nu}_{\mu}}F^{\rm
fl}_{{\nu}_{\mu}}=E_{{\nu}_{\mu}}F^{q}_{E_{{\nu}_{\mu}}}\left(\frac{F^{\rm
fl}_{\gamma}}{F^{q}_{\gamma}}\right)^{\gamma}=A_{X}\,\frac{E_{X}F^{q}_{E_{X}}}{3}\left(\frac{F^{\rm
fl}_{\gamma}}{F^{q}_{\gamma}}\right)^{\gamma}$ (7)

where $A_{X}$ is a normalization parameter. For $\gamma=1.0-2.0$, we can set
an upper limit on $A_{X}$ so that our predictions are consistent with the non-
detection 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 years of _Swift_ -XRT data Giommi et al.
(2019). The 0.5 keV, 1.0 keV, 1.5 keV, 3.0 keV, 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^{q}_{E_{X}}$, is then defined as the
average of the quiescent fluxes in these five energies. Figure D.1 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 Lacs, 6
BCUs), ignoring 16 sources without flaring states found at the $6\sigma$
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 ${\gamma}=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 one or two 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 $\gamma$ based on the non-detection of
neutrino events, compared to the IceCube 90 % sensitivity for one-week binning
(left) and 10-year binning (right). The red, blue, green, and orange lines
show sources with the lowest normalization parameters (see inset legend).

Figure 10: Estimated muon neutrino flare fluxes of 129 blazars (97 FSRQs, 26
BL Lacs, and 6 BCUs) for a one-week bin (a) and a 10-year bin (b) of scenario
2 with $A_{X}=1.0$ and ${\gamma}=1.5$ as a function of the ${\sin}({\delta})$
with declination ${\delta}$. The 90 % IceCube sensitivities are the same as in
Figure 8. Figure 11: Upper limits on the normalization parameter $A_{X}$ of
the analyzed blazars in gray as a function of the indices $\gamma$ based on
the non-detection of neutrino events, compared to the IceCube 90% sensitivity
for scenario 2. Left and right panels show the normalization parameters with
one week binning and 10 year binning, respectively. For sources whose flares
are all less than the sensitivity flux, the normalization parameter $A_{X}$ is
set to 1.0. The normalization parameter of the source is adjusted so that the
maximum neutrino flux of the source is equal to the sensitivity flux if the
flux of at least one flare exceeds the sensitivity. Several sources with the
lowest normalization parameters are plotted with red, blue, green, and orange
lines.

###  4.3 Impacts of realistic spectra and upper limits

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-year binning to test this variant of
scenarios 1 and 2 respectively, and to place upper limits on the normalization
parameters, $A_{\gamma}$ and $A_{X}$. To calculate the expected number of
neutrinos, we assume the IC86 point source effective area given in Aartsen et
al. (2017).

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 two lower than that
of an $E_{\nu}^{-2}$ neutrino spectrum. Therefore the respective upper limits
of the more realistic model are a factor of two 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 than
with a different spectral shape for a given total neutrino energy flux.
Regardless of the neutrino spectral template, the upper limits on $A^{\rm
up}_{\gamma/X}$ 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). The assumed upper limit is 2.44 neutrinos in the livetime of IceCube. The IC86 point source effective area has been assumed. |  |  | LMBB2b | $E_{\nu}^{-2}$
---|---|---|---|---
Source | Scenario | $\gamma$ |  $N_{\nu_{\mu}+\bar{\nu}_{\mu}}(A_{\gamma/X}=1$) |  $A_{\gamma/X}^{\rm up}$ (90% CL) |  $N_{\nu_{\mu}+\bar{\nu}_{\mu}}(A_{\gamma/X}=1$ ) |  $A_{\gamma/X}^{\rm 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\times 10^{-3}$
3C 454.3 | 1 | 2.0 | 2000 | $10^{-3}$ | 6000 | $4\times 10^{-4}$
3C 454.3 | 2 | 1.0 | 10 | 0.2 | 30 | 0.08
3C 454.3 | 2 | 1.5 | 100 | $2\times 10^{-2}$ | 300 | $8\times 10^{-3}$
3C 454.3 | 2 | 2.0 | 1000 | $2\times 10^{-3}$ | 3000 | $8\times 10^{-4}$

###  4.4 Contribution of blazar flares to the all-sky neutrino flux

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., 2017; 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
declination of each blazar in our sample.

As shown in Figures 8 and 10, some neutrino flare fluxes with $A_{X}=1.0$ and
$A_{\gamma}=1.0$ exceed the IceCube 90 % sensitivity. The normalization
parameters of $A_{\gamma}$ 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 catalogue 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{\pi}$ and then multiplied by two factors for the completeness correction.
First, the correction factor of all 4LAC blazars and the blazars of our sub-
sample is obtained by the fraction of the total gamma-ray energy fluxes in 0.1
– 316 GeV (raised to the power of ${\gamma}$) of 4LAC blazars to the sample
blazars:

${\sum_{i=1}^{N_{c}}({\rm EnergyFlux}_{i})^{\gamma}}/{\sum_{j=1}^{N_{s}}({\rm
EnergyFlux}_{j})^{\gamma}}$ (8)

where $N_{c}$ is the number of blazars in the 4LAC catalog, $N_{s}$ is the
number of blazars in our sample, ${\rm EnergyFlux}_{i,j}$ are time-average
gamma-ray energy fluxes in 0.1 – 316 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 $\gamma$, 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.

Table 2: Correction factors for determining the contribution of blazar flares to the all-sky neutrino flux. | 4LAC/Sample 888The correction factor of all the 4LAC FSRQs/BL Lacs | Total/4LAC 999The correction factor of the whole population of blazars
---|---|---
$\gamma$ |  FSRQs |  BL Lacs |
1.0 | 1.76 | 4.24 | 3.87
1.5 | 1.22 | 2.21 | 1.09
2.0 | 1.06 | 1.60 | 1.01

Figure 12 shows the total fluxes of muon neutrinos estimated from our bright
gamma-ray FSRQs and BL Lacs, 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 Lacs of scenarios 1 and 2 are plotted as a function of
$\gamma$, compared to a reference flux of the isotropic diffuse muon neutrino
flux (Aartsen et al., 2016b) 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 $\gamma$, the sum of the
neutrino energy fluxes of the sample blazars gradually increases with
$\gamma$. On the other hand, the correction factors sharply decrease with
$\gamma$ up to ${\gamma}~{}{\simeq}~{}1.5$ and stay nearly constant for
$\gamma~{}{\gtrsim}~{}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 $\gamma$ up to
${\gamma}~{}{\simeq}~{}1.5$ and turn to increase gradually with $\gamma$ for
$\gamma~{}{\gtrsim}~{}1.5$ as shown in scenario 1 of Fig. 12. On the other
hand, the total neutrino energy fluxes are maintained nearly constant for
$\gamma~{}{\gtrsim}~{}1.5$ as shown in scenario 2 of Fig. 12, since the gamma-
ray normalization factors $A_{\gamma}$ are strongly suppressed in tens of
blazars shown in Fig. 9 (right). The upper limits of scenario 1 are in
agreement with those of Hooper et al. (2019) for $\gamma\gtrsim 1.5$ even
though the analyses and assumptions are different. We predict weaker upper
limits though for smaller values of $\gamma$. This is because the second
correction factor should be larger for small values of $\gamma$ (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 $\gamma$.

Figure 12: The total muon neutrino fluxes of gamma-ray flaring FSRQs and BL
Lacs as a function of $\gamma$. The total fluxes of muon neutrinos are
estimated from our samples, multiplied by correction factors of the
contribution of the 4LAC blazars and not-yet-detected blazars that are too
gamma-ray faint to be included in the 4LAC catalog. The results of scenarios 1
and 2 are shown in blue solid squares and orange solid circles, compared to a
reference flux of all-sky muon neutrinos (gray dashed line) (Aartsen et al.,
2016b). The three upper limits of Hooper et al. (2019) derived for different
assumptions about the source evolution ($1/D_{L}^{2}$ , Flat, $\gamma$-Ray),
and reduced to $1/3$ for muon neutrinos, are overplotted for comparison.

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

Figure 13: The upper limits of $E_{{\nu}_{\mu}}^{2}{\Phi}_{{\nu}_{\mu}}$ of a
10-year binning as a function of the flare significance $s$ above the
quiescent level for FSRQs (open blue squares) and BL Lacs (open orange
circles) of scenario 1 and FSRQs (solid blue squares) and BL Lacs (solid
orange circles) of scenario 2 with $\gamma=1.5$. The dashed horizontal line
shows a reference flux of the isotropic diffuse muon neutrinos (Aartsen et
al., 2016b).

##  5 Summary and Conclusions

We analyzed 145 bright gamma-ray blazars among the 4LAC blazars. Their flare
duty cycles and energy fractions represent power-law-like distributions,
correlating strongly with each other. We found a significant difference
between blazar sub-classes for the flare duty cycles at the 5 % significant
level. The flare duty cycles and energy fractions also do not correlate with
the quiescent gamma-ray luminosities. Using monthly-binned light curves of the
2FGL catalog, Ackermann et al. (2011) showed that bright FSRQs and BL Lacs
have flare duty cycles of about $0.05-0.10$. Our flare duty cycles of the
weekly-averaged time bin light curves show much broader distributions for both
blazar sub-classes, ranging from 0.0 to 0.6. Abdollahi et al. (2017) present
the second catalog of flaring gamma-ray sources (2FAV), whose Table 2 shows
the number of weekly-bin flares detected for each source for the first 387
weeks of Fermi observations. According to Table 2 in Abdollahi et al. (2017),
the flare duty cycles of the FAVA analysis appear to suppress less than
${\sim}0.2$, while our flare duty cycles extend to ${\sim}0.6$. As described
in appendix E, the photon flux distribution from the FAVA analysis is likely
to systematically underestimate the flaring contribution in the total energy
output of blazars. For the blazars analized in our sample, we find the power-
law indices $\alpha$ of the photon flux distributions are less than 2.5. This
suggests that high-energy neutrinos of blazars might be produced mainly during
the flare phase.

We estimated muon neutrino flare fluxes by using a scaling law of
$L_{\nu}^{\rm fl}=L_{\nu}^{q}(L_{\gamma}^{\rm fl}/L_{\gamma}^{q})^{\gamma}$.
The quiescent neutrino flux is benchmarked with two proxies: the quiescent
${\gamma}$-ray flux and the quiescent X-ray flux with the normalization
parameters of $A_{\gamma}$ and $A_{X}$, whose upper limits are 1.0. By
comparing the estimated muon neutrino flare fluxes to the IceCube 90%
sensitivities, $A_{\gamma}$ and $A_{X}$ are restricted to values much lower
than 1.0 for several tens of blazars. As shown in Figure 10, the predicted
neutrino flaring flux for TXS 0506+056 is below the sensitivity in scenario 2,
which is motivated by this source. The expected number of muon neutrinos from
TXS 0506+056 is $\sim 1$ as presented in Table 1, which is below the 90%
sensitivity corresponding to $\sim 2.4$ events.

The origin of all-sky neutrinos observed in IceCube is one of the most
important puzzles in high-energy neutrino astrophysics. We found that
scenarios 1 and 2 suggest that no more than ${\sim}50$ % and ${\sim}14$ % of
the all-sky neutrino flux can originate from gamma-ray flares of FSRQs and BL
Lacs, respectively. A more realistic neutrino spectrum than the usual
$E_{\nu}^{-2}$ power law yields upper limits that are a factor of two less
constraining. The upper limits are consistent with those obtained the previous
literature despite different methods and assumptions.

## Acknowledgments

We thank Chengchao Yuan for providing the data of the completeness factor.
K.M. is supported by the NSF Grant No. AST-1908689, No. AST-2108466 and No.
AST-2108467, and KAKENHI No. 20H01901 and No. 20H05852. M.P. acknowledges
support from the MERAC Fondation through the project THRILL and from the
Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “2nd call
for H.F.R.I. Research Projects to support Faculty members and Researchers”
through the grant number 3013 (UNTRAPHOB).

## References

  * Aartsen et al. (2013a) Aartsen, M., et al. 2013a, Phys.Rev.Lett., 111, 021103, doi: 10.1103/PhysRevLett.111.021103
  * Aartsen et al. (2013b) —. 2013b, Science, 342, 1242856, doi: 10.1126/science.1242856
  * Aartsen et al. (2014) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014, ApJ, 796, 109, doi: 10.1088/0004-637X/796/2/109
  * Aartsen et al. (2015a) Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2015a, Phys. Rev. Lett., 115, 081102, doi: 10.1103/PhysRevLett.115.081102
  * Aartsen et al. (2015b) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2015b, ApJ, 807, 46, doi: 10.1088/0004-637X/807/1/46
  * Aartsen et al. (2016a) Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2016a, ApJ, 833, 3, doi: 10.3847/0004-637X/833/1/3
  * Aartsen et al. (2016b) —. 2016b, ApJ, 833, 3, doi: 10.3847/0004-637X/833/1/3
  * Aartsen et al. (2017) —. 2017, ApJ, 835, 45, doi: 10.3847/1538-4357/835/1/45
  * Aartsen et al. (2017) Aartsen, M. G., et al. 2017, Astrophys. J., 835, 151, doi: 10.3847/1538-4357/835/2/151
  * Aartsen et al. (2020a) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020a, Phys. Rev. Lett., 125, 121104, doi: 10.1103/PhysRevLett.125.121104
  * Aartsen et al. (2020b) —. 2020b, Phys. Rev. Lett., 124, 051103, doi: 10.1103/PhysRevLett.124.051103
  * Abbasi et al. (2022) Abbasi, R., et al. 2022.  https://arxiv.org/abs/2207.04946
  * Abdollahi et al. (2017) Abdollahi, S., Ackermann, M., Ajello, M., et al. 2017, ApJ, 846, 34, doi: 10.3847/1538-4357/aa8092
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Ackermann et al. (2011) Ackermann, M., Ajello, M., Allafort, A., et al. 2011, ApJ, 743, 171, doi: 10.1088/0004-637X/743/2/171
  * Ahlers & Halzen (2017) Ahlers, M., & Halzen, F. 2017, Progress of Theoretical and Experimental Physics, 2017, doi: 10.1093/ptep/ptx021
  * Ajello et al. (2020) Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105, doi: 10.3847/1538-4357/ab791e
  * Albert et al. (2018) Albert, A., André, M., Anghinolfi, M., et al. 2018, ApJ, 868, L20, doi: 10.3847/2041-8213/aaeecf
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Atwood et al. (2009) Atwood, W. B., et al. 2009, Astrophys. J., 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Buson et al. (2022) Buson, S., Tramacere, A., Pfeiffer, L., et al. 2022, arXiv e-prints, arXiv:2207.06314.  https://arxiv.org/abs/2207.06314
  * Capel et al. (2020) Capel, F., Mortlock, D. J., & Finley, C. 2020, Phys. Rev. D, 101, 123017, doi: 10.1103/PhysRevD.101.123017
  * Costamante (2020) Costamante, L. 2020, MNRAS, 491, 2771, doi: 10.1093/mnras/stz3018
  * Dermer et al. (2014) Dermer, C. D., Murase, K., & Inoue, Y. 2014, JHEAp, 3-4, 29, doi: 10.1016/j.jheap.2014.09.001
  * Gao et al. (2017) Gao, S., Pohl, M., & Winter, W. 2017, Astrophys. J., 843, 109, doi: 10.3847/1538-4357/aa7754
  * Giommi et al. (2019) Giommi, P., Brandt, C. H., Barres de Almeida, U., et al. 2019, A&A, 631, A116, doi: 10.1051/0004-6361/201935646
  * Hooper et al. (2019) Hooper, D., Linden, T., & Vieregg, A. 2019, J. Cosmology Astropart. Phys, 2019, 012, doi: 10.1088/1475-7516/2019/02/012
  * Hovatta et al. (2021) Hovatta, T., et al. 2021, Astron. Astrophys., 650, A83, doi: 10.1051/0004-6361/202039481
  * IceCube Collaboration et al. (2018a) IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018a, Science, 361, eaat1378, doi: 10.1126/science.aat1378
  * IceCube Collaboration et al. (2018b) —. 2018b, Science, 361, 147, doi: 10.1126/science.aat2890
  * Kadler et al. (2016) Kadler, M., et al. 2016, Nature Phys., 12, 807, doi: 10.1038/NPHYS3715
  * Keivani et al. (2018) Keivani, A., et al. 2018, Astrophys. J., 864, 84, doi: 10.3847/1538-4357/aad59a
  * Mészáros (2017) Mészáros, P. 2017, Annual Review of Nuclear and Particle Science, 67, 45, doi: 10.1146/annurev-nucl-101916-123304
  * Meyer et al. (2019) Meyer, M., Scargle, J. D., & Blandford, R. D. 2019, ApJ, 877, 39, doi: 10.3847/1538-4357/ab1651
  * Murase (2017) Murase, K. 2017, Active Galactic Nuclei as High-Energy Neutrino Sources, ed. T. Gaisser & A. Karle, 15–31, doi: 10.1142/9789814759410_0002
  * Murase et al. (2013) Murase, K., Ahlers, M., & Lacki, B. C. 2013, Phys. Rev. D, 88, 121301, doi: 10.1103/PhysRevD.88.121301
  * Murase et al. (2016) Murase, K., Guetta, D., & Ahlers, M. 2016, Phys. Rev. Lett., 116, 071101, doi: 10.1103/PhysRevLett.116.071101
  * Murase et al. (2014) Murase, K., Inoue, Y., & Dermer, C. D. 2014, Phys. Rev. D, 90, 023007, doi: 10.1103/PhysRevD.90.023007
  * Murase et al. (2018) Murase, K., Oikonomou, F., & Petropoulou, M. 2018, ApJ, 865, 124, doi: 10.3847/1538-4357/aada00
  * Murase & Stecker (2022) Murase, K., & Stecker, F. W. 2022, arXiv e-prints, arXiv:2202.03381.  https://arxiv.org/abs/2202.03381
  * Murase & Waxman (2016) Murase, K., & Waxman, E. 2016, Phys. Rev. D, 94, 103006, doi: 10.1103/PhysRevD.94.103006
  * Nalewajko (2013) Nalewajko, K. 2013, MNRAS, 430, 1324, doi: 10.1093/mnras/sts711
  * Neronov et al. (2017) Neronov, A., Semikoz, D. V., & Ptitsyna, K. 2017, A&A, 603, A135, doi: 10.1051/0004-6361/201630098
  * Oberto et al. (2019) Oberto, A., Allen, M., Loup, C., & Lesteven, S. 2019, in Astronomical Society of the Pacific Conference Series, Vol. 521, Astronomical Data Analysis Software and Systems XXVI, ed. M. Molinaro, K. Shortridge, & F. Pasian, 96
  * Oikonomou et al. (2019) Oikonomou, F., Murase, K., Padovani, P., Resconi, E., & Mészáros, P. 2019, Mon. Not. Roy. Astron. Soc., 489, 4347, doi: 10.1093/mnras/stz2246
  * Oikonomou et al. (2021) Oikonomou, F., Petropoulou, M., Murase, K., et al. 2021, JCAP, 10, 082, doi: 10.1088/1475-7516/2021/10/082
  * Padovani et al. (2015) Padovani, P., Petropoulou, M., Giommi, P., & Resconi, E. 2015, MNRAS, 452, 1877, doi: 10.1093/mnras/stv1467
  * Paliya et al. (2020) Paliya, V. S., Böttcher, M., Olmo-García, A., et al. 2020, Astrophys. J., 902, 29, doi: 10.3847/1538-4357/abb46e
  * Palladino et al. (2019) Palladino, A., Rodrigues, X., Gao, S., & Winter, W. 2019, ApJ, 871, 41, doi: 10.3847/1538-4357/aaf507
  * Petropoulou et al. (2016) Petropoulou, M., Coenders, S., & Dimitrakoudis, S. 2016, Astropart. Phys., 80, 115, doi: 10.1016/j.astropartphys.2016.04.001
  * Petropoulou et al. (2015) Petropoulou, M., Dimitrakoudis, S., Padovani, P., Mastichiadis, A., & Resconi, E. 2015, MNRAS, 448, 2412, doi: 10.1093/mnras/stv179
  * Petropoulou et al. (2020) Petropoulou, M., Oikonomou, F., Mastichiadis, A., et al. 2020, Astrophys. J., 899, 113, doi: 10.3847/1538-4357/aba8a0
  * Petropoulou et al. (2020) Petropoulou, M., Murase, K., Santander, M., et al. 2020, ApJ, 891, 115, doi: 10.3847/1538-4357/ab76d0
  * Plavin et al. (2021) Plavin, A. V., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2021, ApJ, 908, 157, doi: 10.3847/1538-4357/abceb8
  * Resconi et al. (2009) Resconi, E., Franco, D., Gross, A., Costamante, L., & Flaccomio, E. 2009, A&A, 502, 499, doi: 10.1051/0004-6361/200911770
  * Rodrigues et al. (2021) Rodrigues, X., Garrappa, S., Gao, S., et al. 2021, Astrophys. J., 912, 54, doi: 10.3847/1538-4357/abe87b
  * Sahakyan et al. (2022) Sahakyan, N., Giommi, P., Padovani, P., et al. 2022.  https://arxiv.org/abs/2204.05060
  * Scargle et al. (2013) Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ApJ, 764, 167, doi: 10.1088/0004-637X/764/2/167
  * Smith et al. (2021) Smith, D., Hooper, D., & Vieregg, A. 2021, J. Cosmology Astropart. Phys, 2021, 031, doi: 10.1088/1475-7516/2021/03/031
  * Stathopoulos et al. (2022) Stathopoulos, S. I., Petropoulou, M., Giommi, P., et al. 2022, MNRAS, 510, 4063, doi: 10.1093/mnras/stab3404
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * Vandenbroucke et al. (2019) Vandenbroucke, J., Pizzuto, A., & Meagher, K. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 1026
  * Wood et al. (2017) Wood, M., Caputo, R., Charles, E., et al. 2017, in International Cosmic Ray Conference, Vol. 301, 35th International Cosmic Ray Conference (ICRC2017), 824.  https://arxiv.org/abs/1707.09551
  * Yuan et al. (2020) Yuan, C., Murase, K., & Mészáros, P. 2020, ApJ, 890, 25, doi: 10.3847/1538-4357/ab65ea
  * Zhou et al. (2021) Zhou, B., Kamionkowski, M., & Liang, Y.-f. 2021, Phys. Rev. D, 103, 123018, doi: 10.1103/PhysRevD.103.123018

##  Appendix A List of Our Sources

A list of the sources analyzed in this work is given in Table
LABEL:tab:source_list, which is available in its entirety in machine-readable
form. A partial list is presented here.

Table A.1: List of our sources 4FGL Name | Object Name | RAJ2000 | DEJ2000 | Optical Class | SED Class | $z$
---|---|---|---|---|---|---
4FGL J0001.5+2113 | TXS 2358+209 | 0.3815 | 21.2183 | FSRQ | ISP | 1.106
4FGL J0017.5-0514 | PMN J0017-0512 | 4.3949 | -5.2347 | FSRQ | LSP | 0.227
4FGL J0038.2-2459 | PKS 0035-252 | 9.5652 | -24.9899 | FSRQ | LSP | 1.196
4FGL J0102.8+5824 | TXS 0059+581 | 15.701 | 58.4092 | FSRQ | LSP | 0.644
4FGL J0108.6+0134 | 4C +01.02 | 17.1695 | 1.5819 | FSRQ | LSP | 2.099
4FGL J0112.8+3208 | 4C +31.03 | 18.2227 | 32.1399 | FSRQ | LSP | 0.603
4FGL J0132.7-1654 | PKS 0130-17 | 23.176 | -16.9103 | FSRQ | LSP | 1.02
4FGL J0133.1-5201 | PKS 0131-522 | 23.2938 | -52.0202 | FSRQ | LSP | 0.925
4FGL J0210.7-5101 | PKS 0208-512 | 32.6946 | -51.0218 | FSRQ | LSP | 1.003
$\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$
Notes. Optical and SED classes are from Fourth LAT AGN Catalog (4LAC) Ajello
et al. (2020) . Redshifts are from the 4LAC and SIMBAD Astronomical Database
Oberto et al. (2019). A portion is shown here. (This table is available in its
entirety in machine-readable form.)

##  Appendix B Examples of 1-316 GeV gamma-ray light curves

Figure B.1 presents examples of weekly-binned _Fermi_ LAT light curves in the
higher energy range of 1-316 GeV (black points). The Bayesian blocks
representations and the quiescent flux levels, which are derived by using the
method described in section 3.1.1, are also shown in solid orange lines and
dashed cyan lines.

Figure B.1:  Indicative examples of weekly-binned blazar gamma-ray light
curves (black points) in the 1-316 GeV energy range. Grey triangles indicate
95% upper limits. See the caption of Figure 2 in details.

##  Appendix C Gamma-ray quiescent flux levels with time binning

As mentioned in Section 3, even though there is no widely accepted agreement
on how to determine the quiescent flux, it is reasonable to assume that it
remains low and constant. When the gamma-ray flux is assumed to be low and
constant, the quiescent gamma-ray flux does not increase with coarser time
binning unless the flaring state is in effect. However, the coarser time
binning should cause an increase in the quiescent gamma-ray flux, accounting
for the flaring fluxes. Figure C.1 presents the mean gamma-ray quiescent
fluxes of 106 FSRQs, 31 BL Lacs, and 8 BCUs with standard deviations as a
function of time bins. The quiescent fluxes remain stable below $\sim$20 weeks
and gradually increase above $\sim$20 weeks, indicating that they can
reasonably be determined to be relatively low and constant using a 1-week time
bin.

Figure C.1: The mean gamma-ray quiescent fluxes of 106 FSRQs (blue squares),
31 BL Lacs (orange circles), and 8 BCUs (green triangles) with the standard
deviations as a function of time bins.

##  Appendix D X-ray quiescent flux level

Examples of 1.0 keV _Swift_ X-ray light curves are shown in Fig. D.1 (black
points). The Bayesian blocks representations are overplotted with solid orange
lines. The quiescent flux levels are also presented with dashed cyan lines and
are estimated in the same way as for the gamma-ray quiescent flux levels – see
Section 3.1.1.

Figure D.1:  Examples of the 1.0 keV quiescent fluxes (dashed cyan lines) of
PKS 0208-512 and PKS 1510-089 X-ray light curves (black points) with the
Bayesian block representations (orange solid lines).

##  Appendix E Comparison with photon flux distribution obtained by the FAVA
analysis

The photon flux distribution has been used to study the time variation to hold
some of the important clues to the origin and nature of their variability.
Murase et al. (2018) used the flux distributions with a convolution of a power
law with a Poisson distribution, selecting 6 BL Lacs at intermediate redshifts
from the FAVA analysis. Their analysis indicates that the power-law index
$\alpha$ ranges from 1.7 to 3.0. However, it has been known that the FAVA
analysis is not ideal for obtaining detailed information on gamma-ray light
curves. Smaller values of $\alpha$ imply the larger output of neutrinos during
the flaring phase, so it is important to determine $f_{\rm fl}$ and $b_{\rm
fl}^{\gamma}$ more precisely.

In the standard leptonic scenario for the blazar gamma-ray emission, in
common, the leptonic models predict a neutrino luminosity $L_{\nu}$ is in
proportion to a gamma-ray luminosity $L_{\gamma}$ to the power of $\gamma\sim
1.5-2.0$ (see Murase & Waxman (2016), and references therein), giving

$L_{\nu}^{2}\frac{dN}{dL_{\nu}}\propto
L_{\nu}^{\frac{\gamma-\alpha+1}{\gamma}},$ (E1)

which implies that the flaring contribution can be dominant for the larger
$\gamma$ and smaller $\alpha$, e.g., $\gamma\geq 1.5$ and $\alpha\leq 2.5$
(Murase et al., 2018). Thus, the index $\alpha$ is a good indicator of the
flaring contribution to the neutrino output. In this section, we compare our
results of the flux distributions with those from the FAVA analysis.

Our results for several indicative sources are presented in Table E.1. The
source flux distribution per week is modeled as a convolution of a power law
with a Poisson distribution of the source photons plus the background photons.
While the number of the source photons is related to the source flux with the
effective area and the exposure time, the effective area multiplied by the
exposure time changes one-week bin by one-week bin. The number of the
background photons changes as well. Hence, we include their systematic
uncertainties to the fitting errors as the sum in quadrature of the errors.
Figure E.1 shows an example of the gamma-ray flux distributions fitted with a
power-law function convolved with a Poisson distribution including the
background photons, where the fitting is carried out with an unbinned maximum
likelihood method. Our results on $\alpha$, ranging from 1.2 to 2.4, are
systematically smaller than those of FAVA analysis, which range from 1.7 to
3.0 (Murase et al., 2018). Hence, the photon flux distribution from the FAVA
analysis tends to underestimate the flaring contribution to the neutrino
output of a blazar, and our results would better justify the “flare dominance”
in neutrino emission, considered by Murase et al. (2018).

Figure E.1:  The gamma-ray flux distributions of TXS 0506+056, PKS 0426-380,
OJ 287, and PKS 0235+164 in 0.1–316 GeV fitted with a power-law function
convolved with a Poisson distribution including the background photons. See
text for details.  Table E.1: Gamma-ray flare duty cycle ($f_{\rm fl}$), flare
energy fraction ($b^{\gamma}_{\rm fl}$), and best-fit power-law index of the
flux distribution ($\alpha$) for a few indicative sources from our sample.

Name | $f_{\rm fl}$ | $b^{\gamma}_{\rm fl}$ | $f^{\rm LE}_{\rm fl}$ | $b^{\gamma{\rm LE}}_{\rm fl}$ | $f^{\rm HE}_{\rm fl}$ | $b^{\gamma{\rm HE}}_{\rm fl}$ | $\alpha$ | $\alpha^{\rm HE}$
---|---|---|---|---|---|---|---|---
TXS 0506+056 | $0.014\pm 0.002$ | $0.055\pm 0.011$ | $0.009\pm 0.002$ | $0.039\pm 0.009$ | $0.038\pm 0.004$ | $0.132\pm 0.016$ | $2.05_{-0.31}^{+0.27}$ | $2.82_{-0.16}^{+0.15}$
PKS 0426-380 | $0.164\pm 0.008$ | $0.364\pm 0.020$ | $0.136\pm 0.008$ | $0.322\pm 0.019$ | $0.152\pm 0.008$ | $0.343\pm 0.020$ | $1.23_{-0.18}^{+0.28}$ | $2.18_{-0.20}^{+0.43}$
OJ 287 | $0.031\pm 0.004$ | $0.115\pm 0.015$ | $0.029\pm 0.004$ | $0.110\pm 0.015$ | $0.023\pm 0.003$ | $0.113\pm 0.019$ | $2.13_{-0.16}^{+0.12}$ | $2.84_{-0.19}^{+0.18}$
PKS 0235+164 | $0.059\pm 0.005$ | $0.183\pm 0.017$ | $0.051\pm 0.005$ | $0.165\pm 0.017$ | $0.070\pm 0.006$ | $0.212\pm 0.019$ | $1.58_{-0.15}^{+0.13}$ | $1.87_{-0.09}^{+0.13}$
PKS 0301-243 | $0.014\pm 0.002$ | $0.076\pm 0.016$ | $0.010\pm 0.002$ | $0.067\pm 0.016$ | $0.009\pm 0.002$ | $0.058\pm 0.017$ | $2.43_{-0.22}^{+0.21}$ | $3.15_{-0.16}^{+0.43}$
S5 0716+71 | $0.290\pm 0.011$ | $0.507\pm 0.021$ | $0.221\pm 0.010$ | $0.415\pm 0.020$ | $0.224\pm 0.010$ | $0.463\pm 0.022$ | $1.43_{-0.13}^{+0.19}$ | $2.15_{-0.12}^{+0.13}$
S4 0954+65 | $0.024\pm 0.003$ | $0.114\pm 0.018$ | $0.021\pm 0.003$ | $0.104\pm 0.018$ | $0.021\pm 0.003$ | $0.138\pm 0.024$ | $2.37_{-0.17}^{+0.13}$ | $2.51_{-0.17}^{+0.22}$

##  Appendix F Estimated muon neutrino flare fluxes by using the gamma-ray
photon fluxes.

Figure F.1 shows the estimated muon neutrino flare fluxes from the gamma-ray
photon fluxes for a one-week bin (a) and a 10-year bin (b) of scenario 1 with
$A_{\gamma}=1.0$ and ${\gamma}=1.5$ as a function of ${\sin}({\delta})$. As
seen by comparing Fig. F.1 with Fig. 8, the estimated muon neutrino fluxes of
the gamma-ray photon fluxes are almost similar to those of the gamma-ray
energy fluxes.

Figure F.1: Estimated muon neutrino flare fluxes of 136 blazars (101 FSRQs, 28
BL Lacs, and 7 BCUs) from the gamma-ray photon fluxes for a one-week bin (a)
and a 10-year bin (b) of scenario 1 with $A_{\gamma}=1.0$ and ${\gamma}=1.5$
as a function of ${\sin}({\delta})$. The 90 % IceCube sensitivities are the
same as in Figure 8.