Blazar Flares as an Origin of High-energy Cosmic Neutrinos?

The blazar contribution to the INB has been constrained by different types of analyses: (i) diffuse searches for extremely high energy (EHE) neutrinos (Aartsen et al. 2016, 2017a), (ii) event clustering and autocorrelation analyses (Aartsen et al. 2015, 2017b; Murase & Waxman 2016), and (iii) stacking and cross-correlation analyses (Aartsen et al. 2017c, 2017b).

Neutrino spectra have been predicted by most theoretical models to be hard and peaking at energies beyond 1 PeV. The hardness of the spectrum is related to the fact that the target photon density in blazars is higher at lower energies. Upon normalization to the IceCube flux at ∼1 PeV, most model-predicted fluxes at 10 PeV are found to be ${E}_{\nu }^{2}{{\rm{\Phi }}}_{\nu }\gtrsim (3\mbox{--}5)\times {10}^{-8}\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}$ for all flavors (Murase 2017, and references therein). With the 9 yr diffuse analysis, the IceCube Collaboration reported an upper limit on the INB, ${E}_{\nu }^{2}{{\rm{\Phi }}}_{\nu }\lesssim 1\times {10}^{-8}\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}$ (90% CL) at 5–10 PeV (Aartsen et al. 2017a), excluding some of the optimistic physical models for blazar neutrino emission (see also Neronov et al. 2017).

Another type of less model-dependent constraint is obtained from the absence of sources of high-energy multiplets  (Ahlers & Halzen 2014; Murase & Waxman 2016; see also Lipari 2008; Silvestri & Barwick 2010; Murase et al. 2012a, for earlier works). Let the number of the sources with multiplets be denoted as ${N}_{m\geqslant k}$. Then, constraints can be placed by requiring ${N}_{m\geqslant k}\lt 1$ at sufficiently high energies (e.g., >50 TeV for muons). Although the energy-dependent effective area should be taken into account for detailed calculations as in Murase & Waxman (2016), the basic results can be understood by using a limit from the nondetection of point sources. The 8 yr point-source sensitivity (90% CL) for an E⁻² neutrino spectrum is F_lim ∼ (5–6) × 10⁻¹⁰ $\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ (Aartsen et al. 2017b). For such a flat energy spectrum with a time-averaged luminosity of ${\varepsilon }_{\nu }{L}_{{\varepsilon }_{{\nu }_{\mu }}}^{\mathrm{ave}}\sim {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, the number density of the sources is constrained as

$$\begin{eqnarray}\begin{array}{rc}\,{n}_{0}^{\mathrm{eff}} & \lesssim \,1.9\times {10}^{-10}\,{\mathrm{Mpc}}^{-3}\,{\left(\displaystyle \frac{{\varepsilon }_{\nu }{L}_{{\varepsilon }_{{\nu }_{\mu }}}^{\mathrm{ave}}}{{10}^{44}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{-3/2}\\ & \,\times \,{({b}_{m}{q}_{L}/6.6)}^{-1}{F}_{\mathrm{lim},-9.2}^{3/2}\left(\displaystyle \frac{2\pi }{{\rm{\Delta }}{\rm{\Omega }}}\right),\end{array}\end{eqnarray} \tag{ 1 }$$

where q_L ∼ 1−3 is a luminosity-dependent correction factor determined by the redshift evolution and ΔΩ is the solid angle covered by the detector. This limit depends on spectral templates of the sources, and the differential sensitivity can be worse than the integrated sensitivity by $\sim \mathrm{ln}(10)$. Thus, for harder spectra explaining neutrinos only in the PeV range (see Figures 2 and 4 of Murase & Waxman 2016), the upper limit in Equation (1) can be relaxed by ∼3. Note that ${b}_{m}\geqslant 1$ is a factor that represents details of the analysis. For example, in the case of the multiplet analysis, we have b_m ≈ 6.6 for $m\geqslant 2$ and b_m ≈ 1.6 for $m\geqslant 3$.

Using Equation (1), we derive an upper limit on the contribution of a source population with ${n}_{0}^{\mathrm{eff}}$ to the INB:

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\nu }^{2}{{\rm{\Phi }}}_{\nu } & \approx & \displaystyle \frac{{\xi }_{z}{{ct}}_{H}}{4\pi }3({\varepsilon }_{\nu }{L}_{{\varepsilon }_{{\nu }_{\mu }}}^{\mathrm{ave}}){n}_{0}^{\mathrm{eff}}\\ & \lesssim & 6.9\times {10}^{-9}\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}\,\left(\displaystyle \frac{{\xi }_{z}}{0.7}\right){\left(\displaystyle \frac{6.6}{{b}_{m}{q}_{L}}\right)}^{2/3}\\ & & \times {\left(\displaystyle \frac{{n}_{0}^{\mathrm{eff}}}{{10}^{-7}{\mathrm{Mpc}}^{-3}}\right)}^{1/3}{F}_{\mathrm{lim},-9.2}{\left(\displaystyle \frac{2\pi }{{\rm{\Delta }}{\rm{\Omega }}}\right)}^{2/3},\end{array}\end{eqnarray} \tag{ 2 }$$

where t_H is the Hubble time and ξ_z represents the redshift evolution of the neutrino luminosity density of the sources: ξ_z ∼ 0.7 for the γ-ray luminosity density evolution of BL Lacs, ξ_z ∼ 8 for that of FSRQs, and ξ_z ∼ 3 for the X-ray luminosity density evolution of AGNs (Ajello et al. 2014; Ueda et al. 2014). Note that if one uses the number density evolution (corresponding to the equal luminosity weight), we have ξ_z ∼ 0.2 for all BL Lacs with ${n}_{s}\propto {(1+z)}^{-3.5}$ and ξ_z ∼ 0.1 for HSP objects with n_s ∝(1 + z)⁻⁶, respectively (Ajello et al. 2014). In the extreme case where the faintest BL Lacs are equally neutrino emitters, we have ${n}_{0}^{\mathrm{eff}}={n}_{0}^{\mathrm{tot}}\sim (1\mbox{--}3)\times {10}^{-7}\,{\mathrm{Mpc}}^{-3}$ (Ajello et al. 2014) and ξ_z ∼ 0.2, implying that the contribution to the INB is ≲10% at 0.1 PeV. Although the limit could be further relaxed by a factor of two by integrating the number density down to the faintest source tail, blazars are unlikely to be dominant in the INB for such weak redshift evolution. One should keep in mind that the effective number density, which depends on the neutrino luminosity function (i.e., dn_s/dL_ν), should always be smaller than the total number density, i.e., ${n}_{0}^{\mathrm{eff}}\lt {n}_{0}^{\mathrm{tot}}$ (Murase & Waxman 2016). This is because physical models typically predict ${L}_{\nu }\propto {L}_{\gamma }-{L}_{\gamma }^{2}$ (Murase et al. 2014; Petropoulou et al. 2015; Tavecchio & Ghisellini 2015; Murase & Waxman 2016). In the fiducial case of the leptonic scenario for BL Lacs, the effective number density is ${n}_{0}^{\mathrm{eff}}\sim {10}^{-9}\mbox{--}{10}^{-8}\,{\mathrm{Mpc}}^{-3}$, which gives ≲(5–10)(ξ_z/0.7)% of the INB. For FSRQs, we have ${n}_{0}^{\mathrm{eff}}\sim {10}^{-12}\mbox{--}{10}^{-11}\,{\mathrm{Mpc}}^{-3}$, leading to ≲(6–10)(ξ_z/8)% in the 0.1 PeV range. Note that the neutrino luminosity density evolution would be stronger than the γ-ray one if the γ-ray luminosity more strongly weighs on the neutrino luminosity, but such cases are constrained by the stacking limits. In this sense these constraints are complementary to each other.

Another limit can be placed by the autocorrelation analysis on the small-scale anisotropy (Aartsen et al. 2015, 2017b; Ando et al. 2017). With the measured INB and the latest anisotropy limit (Aartsen et al. 2017b), the upper limit on the Poisson angular spectrum is estimated to be ${E}_{\nu }^{4}{C}_{P}\lt 4\,\times {10}^{-19}{\mathrm{GeV}}^{2}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-2}\,{\mathrm{sr}}^{-1}$. Then, using the known formula of C_P for standard candle sources, we obtain ${n}_{0}^{\mathrm{eff}}$ ≲ 1.1 ×${10}^{-9}$ ${\mathrm{Mpc}}^{-3}$ ${({\varepsilon }_{\nu }{L}_{{\varepsilon }_{{\nu }_{\mu }}}^{\mathrm{eff}}/{10}^{44}\mathrm{erg}{{\rm{s}}}^{-1})}^{-3/2}$ ${q}_{L}^{-1}{(2\pi /{\rm{\Delta }}{\rm{\Omega }})}^{1/4}$, which gives a comparable limit.

Blazars are highly variable objects, and as a result, their luminosity density could be dominated by the flaring states characterized by a "flaring" neutrino luminosity ${L}_{\nu }^{\mathrm{fl}}$. In this "flare-dominated" case, the average neutrino luminosity can be written as ${L}_{\nu }^{\mathrm{ave}}\approx {f}_{\mathrm{fl}}{L}_{\nu }^{\mathrm{fl}}$, where f_fl is the duty factor of flares, which will be discussed in the next subsection. The differential neutrino luminosity density is then written as ${\varepsilon }_{\nu }{Q}_{{\varepsilon }_{\nu }}\,=({\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }}^{\mathrm{ave}}){n}_{0}^{\mathrm{eff}}\approx ({\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }}^{\mathrm{fl}})({f}_{\mathrm{fl}}{n}_{0}^{\mathrm{eff}})$. For transients, including flaring sources, the atmospheric background can be reduced owing to the shorter time window, thereby improving, in general, the fluence sensitivity. The power of such a time-dependent search was demonstrated in IceCube Collaboration (2018) (although the excess besides IceCube-170922A was not significant in the time-integrated search). Thus, Equation (1) for the time-averaged emission can still be regarded as a conservative limit, and the constraint given by Equation (2) is still applicable even if blazars are highly variable.

As noted above, stacking and cross-correlation analyses can provide tighter constraints especially for physically motivated models. In particular, for the 2LAC catalog consisting of 862 blazars, the blazar contribution to the INB is restricted to be $\leqslant (19\mbox{--}27)$% (Aartsen et al. 2017c). Note that theory typically predicts ${L}_{\nu }\propto {L}_{\gamma }-{L}_{\gamma }^{2}$ (Murase et al. 2014; Petropoulou et al. 2015; Tavecchio & Ghisellini 2015), and a flux weighting with F_ν ∝ F_γ leads to ≲7% (Aartsen et al. 2017c). For HSP BL Lacs, the preliminary result gives $\leqslant (4.5\mbox{--}5.7)$% (Aartsen et al. 2017b) in the equal flux weight (${F}_{\nu }\propto \mathrm{const}.$) assumption. The stacking limits are powerful and more meaningful when the physical luminosity weight is strong (in which the neutrinos mostly come from the resolved blazars; see Murase et al. 2014). Thus, γ-ray-bright blazars that significantly contribute to the extragalactic γ-ray background are disfavored as the dominant (∼100%) origin of IceCube's neutrinos (see also Wang & Li 2016; Ando et al. 2017; Neronov et al. 2017; Palladino & Vissani 2017; Zhang & Li 2017).

All the constraints discussed so far can be relaxed by a factor of a few under different assumptions. (i) The diffuse EHE limits can be avoided if the cosmic-ray (CR) spectrum is sufficiently soft or the maximum CR energy is lower than ∼10–100 PeV (far below UHECR energies), as considered in Murase et al. (2014) and Dermer et al. (2014). This is because the neutrino production efficiency increases with energy, and the resulting neutrino spectra are hard for a power-law CR spectrum with s_cr ∼ 2−2.6. (ii) The multiplet limit is sensitive to ξ_z. Weakly evolving sources such as BL Lacs are strongly constrained. On the other hand, rapidly evolving sources such as FSRQs could give a significant contribution to the INB if ${n}_{0}^{\mathrm{eff}}\sim {n}_{0}^{\mathrm{tot}}\sim {10}^{-9}\,{\mathrm{Mpc}}^{-3}$. However, this is contrary to the fiducial theoretical prediction, ${n}_{0}^{\mathrm{eff}}\ll {n}_{0}^{\mathrm{tot}}$. Dermer et al. (2014) proposed such a model, in which flaring blazars significantly contribute to the INB only in the PeV range, but this model does not explain the UHECRs (see also Murase 2017). (iii) The stacking limits do not apply to γ-ray-dark blazars. For example, a subset of FSRQs with a spectral energy distribution (SED) peak in the MeV range are dim in the Fermi LAT band, and such MeV blazars could significantly contribute to the INB as hidden CR accelerators (Murase et al. 2016).

The combination of all constraints indicates that the blazar contribution to the INB is likely to be subdominant at least in the 0.1 PeV range. This is even more so the case for the medium-energy component in the 10–100 TeV range, which requires models prohibiting the escape of γ-rays (Murase et al. 2016). On the other hand, at present, the contribution could be more significant in the PeV or higher-energy range. We note that fiducial models (normalized to the UHECR flux), presented in Murase et al. (2014), are consistent with the above constraints and give ∼2%–10% of the INB in the 0.1 PeV range (and more at higher energies).

Bright flaring sources are detectable in neutrinos whether the blazars are dominant or subdominant in the extragalactic neutrino sky (Dermer et al. 2014; Murase & Waxman 2016; Guépin & Kotera 2017). Nevertheless, it is natural to estimate the contribution of blazar flares like the ones from TXS 0506+056 to the INB and discuss the detectability of similar flaring events.

The flaring state lasts only for a fraction of the observation time. For a given time binning, one can measure the number of detected particles (e.g., photons) per bin, which is proportional to luminosity. We introduce the flaring state when the number of photons in a bin exceeds a certain threshold corresponding to the luminosity L^th. Then one can construct the distribution of the number of time bins with luminosity, dN/dL. The fraction of time spent in the flaring state (i.e., the duty factor) is given by

$$\begin{eqnarray}&&{f}_{\mathrm{fl}}=\displaystyle \frac{1}{{N}_{\mathrm{tot}}}{\int }_{{L}^{\mathrm{th}}}{dL}\,\displaystyle \frac{{dN}}{{dL}},\end{eqnarray} \tag{ 3 }$$

where ${N}_{\mathrm{tot}}$ is the total number of time bins. The fraction of energy emitted in the flaring state is

$$\begin{eqnarray}&&{b}_{\mathrm{fl}}=\displaystyle \frac{1}{{L}^{\mathrm{ave}}{N}_{\mathrm{tot}}}{\int }_{{L}^{\mathrm{th}}}{dL}\,L\displaystyle \frac{{dN}}{{dL}},\end{eqnarray} \tag{ 4 }$$

where the average luminosity is given by ${L}^{\mathrm{ave}}\,=(1/{N}_{\mathrm{tot}})\int {dL}\,L({dN}/{dL})$ and the average flaring luminosity is rewritten as L^fl = (b_fl/f_fl)L^ave. Note that in the flare-dominated limit, b_fl ≈ 1, we have ${L}^{\mathrm{ave}}\approx {f}_{\mathrm{fl}}{L}^{\mathrm{fl}}$.

Using the public Fermi All-sky Variability Analysis (FAVA) data by Fermi LAT (Abdollahi et al. 2017), one can obtain the luminosity distribution, dN/dL_γ, of a source, under the assumption that the spectral shape does not change during flaring states (since only photon counts are available in the FAVA analysis). In the latter case, dN/dL_γ ∝ dN/dN_γ. Henceforth, we use dN/dN_γ and dN/dL_γ interchangeably. An example for TXS 0506+056 is shown in Figure 1, in which the photon distribution is modeled as a power law with slope α above a minimum number of photons per bin N_γ,0 that corresponds to a "quiescent" flux. The number of detected photons per time interval is then given by a convolution of this power law with a Poissonian distribution. We applied the same method to a selection of BL Lacs at intermediate redshifts from FAVA and find that they are well described by a power law ${dN}/{{dL}}_{\gamma }\propto {L}_{\gamma }^{-\alpha }$ with α ∼ 2−4. One can also calculate f_fl and b_fl for various blazars using Equations (3) and (4) (see Table 1). The exact values depend on the definition of the flaring state (see, e.g., Resconi et al. 2009), but our main conclusions are not expected to change, if flares are defined differently. We found that the duty factor lies in the range of 0.3%–10% for $\geqslant 5\sigma$ flares (as per the FAVA definition) and obtained f_fl ≈ 0.02–0.1 for TXS 0506+056. The corresponding fraction of emitted photons is ${b}_{\mathrm{fl}}\sim 0.1\ll 1$, implying that the γ-ray emission is dominated by steady emission. Although a detailed statistical study is beyond the purpose of the present work, we underline the need for a systematic study of the properties of flaring, γ-ray-bright blazars. For the purposes of the present work, we treat TXS 0506+056 as a "characteristic" test case.

Zoom In Zoom Out Reset image size

Table 1.  Flare Duty Factors (f_fl) and Flare Energy Fractions (b_fl) for TXS 0506+056, OJ 287, PKS 0426–380, PKS 0301–243, S5 0716+071, and S4 0954+065 as Derived from the FAVA Analysis (Abdollahi et al. 2017)

Name	${L}_{\gamma }$	${f}_{\mathrm{fl}}^{\mathrm{LE}}$	${f}_{\mathrm{fl}}^{\mathrm{HE}}$	${b}_{\mathrm{fl}}^{\mathrm{LE}}$	${b}_{\mathrm{fl}}^{\mathrm{HE}}$	α
TXS 0506+056	1046.3	0.1	0.03	0.1	0.1	3.0
OJ 287	1046.1	0.02	0.01	0.04	0.1	2.9
PKS 0426–380	1048	0.1	0.1	0.2	0.2	1.7
PKS 0301–243	1046	0.04	0.05	0.1	0.3	2.5
S5 0716+071	1046.7	0.07	0.08	0.1	0.2	1.7
S4 0954+065	1045.5	0.04	0.03	0.07	0.3	2.5

Note. Values are reported for the low-energy (LE: 100–800 MeV) and high-energy (HE: 800 MeV–300 GeV) bins. The duty factors quoted are for flaring periods when the flux is enhanced by $\geqslant 5\sigma$ according to the FAVA definition. The "integrated" γ-ray luminosity of the sources, L_γ (in units of erg s⁻¹), is derived from the 3FGL in the 1–100 GeV energy range. Rounded values of the power-law index (α) are also shown.

Download table as:  ASCIITypeset image

During the flare of TXS 0506+056 in the period 2017 September 15–27, the 1–100 GeV "differential" γ-ray flux (${\varepsilon }_{\gamma }{L}_{{\varepsilon }_{\gamma }}^{\mathrm{fl}}\sim 2\times {10}^{46}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$) was about six times higher than the average in the 3FGL catalog (${\varepsilon }_{\gamma }{L}_{{\varepsilon }_{\gamma }}^{\mathrm{ave}}\sim 4\,\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$; Acero et al. 2015; Tanaka et al. 2017).⁷ Thus, this flare can be regarded as one of the brightest flares of this object. Note that the γ-ray photon index in the 0.1–300 GeV range had a similar value during the September flare, β = 2.0, to that in the 3FGL.

Naively, high-energy neutrino emission is dominated by the nonflaring contribution, if one assumes L_ν ∝ L_γ. However, in the standard leptonic scenario for the blazar γ-ray emission, flares can dominate the neutrino output of a blazar. Under different assumptions, the leptonic models predict ${L}_{\nu }\propto {L}_{\gamma }^{\gamma }$ with γ ∼ 1.5–2 (see Murase & Waxman 2016, and references therein), giving

$$\begin{eqnarray}&&{L}_{\nu }^{2}\displaystyle \frac{{dN}}{{{dL}}_{\nu }}\propto {L}_{\nu }^{1-\tfrac{\alpha -1}{\gamma }},\end{eqnarray} \tag{ 5 }$$

which implies that the flaring contribution can be dominant, e.g., for BL Lacs with γ ∼ 2 and $\alpha \leqslant 3$. Such a situation might be the case for the multi-messenger flare associated with IceCube-170922A. In addition, the neutrino flare emission can be pronounced if the CR spectrum is harder during the high state. To demonstrate this, let us consider a toy model where the low-state and high-state CR spectra are described as ${\varepsilon }_{p}{L}_{{\varepsilon }_{p}}^{l}\propto {\varepsilon }_{p}^{2-{s}_{l}}$ (for s_l > 2) and ${\varepsilon }_{p}{L}_{{\varepsilon }_{p}}^{\mathrm{fl}}\propto {\varepsilon }_{p}^{2-{s}_{\mathrm{fl}}}$ (for s_fl < 2), respectively. Assuming that the maximum and minimum energies are ${\varepsilon }_{p}^{\max }$ and ${\varepsilon }_{p}^{\min }$, the flux enhancement at ε_ν is

$$\begin{eqnarray}&&c[{\varepsilon }_{\nu }]\equiv \displaystyle \frac{{\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }}^{\mathrm{fl}}}{{\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }}^{l}}\approx \displaystyle \frac{(2-{s}_{\mathrm{fl}}){f}_{p\gamma }^{\mathrm{fl}}}{({s}_{l}-2){f}_{p\gamma }^{l}}{\left(\displaystyle \frac{20{\varepsilon }_{\nu }}{{\varepsilon }_{p,\max }^{\mathrm{fl}}}\right)}^{2-{s}_{\mathrm{fl}}}{\left(\displaystyle \frac{20{\varepsilon }_{\nu }}{{\varepsilon }_{p,\min }^{l}}\right)}^{{s}_{l}-2},\end{eqnarray} \tag{ 6 }$$

where f_pγ is the effective photomeson production optical depth that we will discuss below. If we adopt indicative values, namely, s_l ∼ 2.3, s_fl ∼ 1.8, ${\varepsilon }_{\nu }\sim 0.05{\varepsilon }_{p,\max }^{\mathrm{fl}}\sim 0.1$ PeV, and ${\varepsilon }_{p,\min }^{l}\sim 10\,\mathrm{GeV}$, we find $c[{\varepsilon }_{\nu }]\sim 30{f}_{p\gamma }^{\mathrm{fl}}/{f}_{p\gamma }^{l}$. Unless the physical conditions during flares change so radically that ${f}_{p\gamma }^{\mathrm{fl}}/{f}_{p\gamma }^{l}\ll 1$, flares may dominate the neutrino output of a blazar. Such a situation could be realized in the 2014–2015 neutrino flare. For either case, in the flare-dominated regime, the time-averaged neutrino luminosity can be written as ${\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }}^{\mathrm{ave}}\approx {f}_{\mathrm{fl}}({\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }}^{\mathrm{fl}})$.

From the above considerations, neutrinos can be copiously produced during the high states. In particular, the contribution of blazar flares like those observed from TXS 0506+056 to the INB is constrained as

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\nu }^{2}{{\rm{\Phi }}}_{\nu }\lesssim 3.8\times {10}^{-10}\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}\,\left(\displaystyle \frac{2\pi }{{\rm{\Delta }}{\rm{\Omega }}}\right)\left(\displaystyle \frac{{\xi }_{z}}{0.7}\right)\\ \,\,\times \,{\left(\displaystyle \frac{0.05}{{f}_{\mathrm{fl}}}\right)}^{1/2}{\left(\displaystyle \frac{{10}^{46}\mathrm{erg}{{\rm{s}}}^{-1}}{{\varepsilon }_{\nu }{L}_{{\varepsilon }_{{\nu }_{\mu }}}^{\mathrm{fl}}}\right)}^{1/2}\left(\displaystyle \frac{6.6}{{b}_{m}{q}_{L}}\right){F}_{\mathrm{lim},-9.2}^{3/2},\,\,\,\end{array}\end{eqnarray} \tag{ 7 }$$

where Equation (2) is used. They can contribute up to a few percent of the INB, allowing uncertainties in the redshift evolution and neutrino spectrum. Equation (7) also implies that dimmer blazar flares could make a larger contribution, but the total neutrino intensity must not exceed the upper limit given by Equation (2).

Finally, we estimate how many bright flares can be detected with multi-messenger (neutrino and γ-ray) coincidence searches in near future. For variable sources, such as TXS 0506+056, we can use the muon neutrino fluence sensitivity, ${\phi }_{\mathrm{lim}}\,\sim 0.04\,\mathrm{GeV}\,{\mathrm{cm}}^{-2}$ (for a spectrum that is broad around 0.1 PeV), which can also be calculated from the public effective area. Due to the γ-ray monitoring of blazars with Fermi LAT in GeV energies and HAWC in TeV energies, the detection rate of flaring blazars in neutrinos can be estimated to be (see Equation (4.14) of Murase et al. 2012a)

$$\begin{eqnarray}\begin{array}{rcl}{\dot{N}}_{\mathrm{blazar}} & \approx & \displaystyle \frac{{\rm{\Delta }}{\rm{\Omega }}}{3}{\rho }_{0}^{\mathrm{eff}}{\left(\displaystyle \frac{{\varepsilon }_{\nu }{L}_{{\varepsilon }_{{\nu }_{\mu }}}^{\mathrm{fl}}{t}_{\mathrm{dur}}}{4\pi {\phi }_{\mathrm{lim}}}\right)}^{3/2}\\ & \lesssim & 0.9\,{\mathrm{yr}}^{-1}\,{\left(\displaystyle \frac{{f}_{\mathrm{fl}}}{0.05}\right)}^{-1/2}{t}_{\mathrm{dur},7}^{1/2}{\left(\displaystyle \frac{{F}_{\mathrm{lim},-9.2}}{{\phi }_{\mathrm{lim},-1.4}}\right)}^{3/2}.\end{array}\end{eqnarray} \tag{ 8 }$$

Here ${\rho }_{0}^{\mathrm{eff}}={f}_{\mathrm{fl}}{n}_{0}^{\mathrm{eff}}/{t}_{\mathrm{dur}}$ is the effective rate density of blazar flares, and the above equation is valid when the IceCube observation time is longer than t_dur/f_fl. This estimate is consistent not only with existing observational constraints but also with the fact that no other flares besides that of TXS 0506+056 have been identified with a high significance. The prospects of detecting neutrinos from short-duration blazar flares are less favorable, because of ${\dot{N}}_{\mathrm{blazar}}\propto {t}_{\mathrm{dur}}^{1/2}$ (see also Petropoulou et al. 2016; Guépin & Kotera 2017). If the association with TXS 0506+056 is physical, according to the standard leptonic scenario including FSRQs (Dermer et al. 2014; Murase et al. 2014), we predict that neutrinos associated with FSRQ flares should also be identified in the future. It also suggests that dedicated time-dependent neutrino searches (Turley et al. 2016, 2018) are important to test these predictions.

Keivani et al. (2018) provided a detailed study of the TXS 0506+056 flare, using the multi-messenger data that have been obtained quasi-simultaneously with IceCube-170922A. The authors found a viable model for both high-energy neutrinos and γ-rays only in the leptonic scenario, where γ-ray emission is attributed to the inverse Compton (IC) mechanism. The X-ray and γ-ray light curves were variable on a day timescale, thus implying a comoving size of the emission region of $l^{\prime} \approx \delta {{ct}}_{\mathrm{var}}/(1+z)\simeq 4.5\times {10}^{16}\,\mathrm{cm}\,(\delta /20){t}_{\mathrm{var},5}$, where typical values of the Doppler factor are δ ∼ 10–30. The observed SED suggests that the Compton dominance parameter is around unity, suggesting a magnetic field of B' ∼ 0.1–1 G(Keivani et al. 2018), which is consistent with population-based estimates for BL Lacs (Celotti & Ghisellini 2008; Murase et al. 2014).

Neutrinos and hadronic γ-rays are coproduced by photomeson production, which is characterized by its effective optical depth, f_pγ. Let us consider a relativistically moving blob and a target photon spectrum, ${n}_{{\varepsilon }_{t}^{{\prime} }}$ (where ${\varepsilon }_{t}^{{\prime} }\approx {\varepsilon }_{t}/\delta$ is the target photon energy in the comoving frame). Approximating the spectrum by ${\varepsilon }_{t}^{{\prime} }{n}_{{\varepsilon }_{t}^{{\prime} }}={n}_{0}^{{\prime} }{({\varepsilon }_{t}^{{\prime} }/{\varepsilon }_{0}^{{\prime} })}^{1-\beta }$ with β > 1, where ${\varepsilon }_{0}^{{\prime} }$ is the reference energy, f_pγ is given by (e.g., Murase et al. 2016)

$$\begin{eqnarray}&&{f}_{p\gamma }[{\varepsilon }_{p}]\approx {\eta }_{p\gamma }[\beta ]{\hat{\sigma }}_{p\gamma }l^{\prime} {n}_{0}^{{\prime} }{({\varepsilon }_{p}^{{\prime} }/{\tilde{\varepsilon }}_{p\gamma 0}^{{\prime} })}^{\beta -1},\end{eqnarray} \tag{ 9 }$$

where ${\eta }_{p\gamma }[\beta ]\approx 2/(1+\beta )$, ${\hat{\sigma }}_{p\gamma }\sim 0.7\times {10}^{-28}\,{\mathrm{cm}}^{2}$ is the attenuation cross section, ${\bar{\varepsilon }}_{{\rm{\Delta }}}\sim 0.3\,\mathrm{GeV}$, and ${\tilde{\varepsilon }}_{p\gamma 0}^{{\prime} }\,=0.5{m}_{p}{c}^{2}{\bar{\varepsilon }}_{{\rm{\Delta }}}/{\varepsilon }_{0}^{{\prime} }$. This estimate is valid when the meson production is dominated by the Δ-resonance and direct pion production. For target photons with observed energy E_t = ε_t/(1 + z), the characteristic energy of protons producing neutrinos of energy ε_ν is ${\varepsilon }_{p}\approx 20{\varepsilon }_{\nu }\approx 0.5{\delta }^{2}{m}_{p}{c}^{2}{\bar{\varepsilon }}_{{\rm{\Delta }}}{{\varepsilon }_{t}}^{-1}$. This results in ${\varepsilon }_{t}\sim 8\,\mathrm{keV}\,{(\delta /20)}^{2}{({\varepsilon }_{\nu }/300\mathrm{TeV})}^{-1}$, corresponding to UV photons or X-rays for neutrino energies ranging from ∼3 PeV to ∼30 TeV. For example, the high-energy synchrotron component of the SED during the TXS 0506+056 flare is well described by β = 2.8 in the optical-to-X-ray range (above the lower-energy peak at ε_syn ∼1 eV; Keivani et al. 2018). In this case, for a CR spectrum with s_cr = 2, the predicted neutrino spectrum is so hard (see also Equation (9)) that it would contradict the nondetection of >10 PeV neutrinos during the flare, unless the CR proton spectrum cuts off at 10–100 PeV (Keivani et al. 2018).

The same target photons lead to the Bethe–Heitler pair production, to which the effective optical depth is

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{BH}}[{\varepsilon }_{p}] & \approx & {\eta }_{p\gamma }[\beta ]{\hat{\sigma }}_{\mathrm{BH}}l^{\prime} {n}_{0}^{{\prime} }{({\varepsilon }_{p}^{{\prime} }/{\tilde{\varepsilon }}_{\mathrm{BH}0}^{{\prime} })}^{\beta -1}\\ & = & g[\beta ]{f}_{p\gamma }[{\varepsilon }_{p}],\end{array}\end{eqnarray} \tag{ 10 }$$

where ${\hat{\sigma }}_{\mathrm{BH}}\sim 0.8\times {10}^{-30}\,{\mathrm{cm}}^{2}$, ${\bar{\varepsilon }}_{\mathrm{BH}}\sim 10(2{m}_{e}{c}^{2})\sim 10$ MeV(Chodorowski et al. 1992), ${\tilde{\varepsilon }}_{\mathrm{BH}0}^{{\prime} }=0.5{m}_{p}{c}^{2}{\bar{\varepsilon }}_{\mathrm{BH}}/{\varepsilon }_{0}^{{\prime} }$, and$g[\beta ]\sim 0.011{(30)}^{\beta -1}$.

The same photons also prevent γ-rays from escaping the source. The γγ optical depth can be written in terms of ε_ν (see Equation (8) in Murase et al. 2016) as

$$\begin{eqnarray}&&{\tau }_{\gamma \gamma }[{\varepsilon }_{\gamma \gamma -p\gamma }]\sim {10}^{3}{f}_{p\gamma }[20\,{\varepsilon }_{\nu }],\end{eqnarray} \tag{ 11 }$$

where ${\varepsilon }_{\gamma \gamma -p\gamma }\sim 10\,\mathrm{GeV}\,({\varepsilon }_{\nu }/300\,\mathrm{TeV})$. The fact that 10–100 GeV γ-rays were observed during the flare suggests that the neutrino production in the same emission region has to be inefficient (e.g., Waxman & Bahcall 1997; Levinson 2006; Dermer et al. 2007; Murase et al. 2016; Petropoulou et al. 2017). Imposing τ_γγ < 1 at 100 GeV leads to ${f}_{p\gamma }\lt {10}^{-3}{({\varepsilon }_{p}/60\mathrm{PeV})}^{\beta -1}$. Note that the neutrino energy is related to the proton energy as ${\varepsilon }_{\nu }\approx 0.05{\varepsilon }_{p}$ = $0.3\,\mathrm{PeV}({\varepsilon }_{p}/6\,\mathrm{PeV})$.

For an isotropic-equivalent proton luminosity, ${\varepsilon }_{p}{L}_{{\varepsilon }_{p}}$, the differential neutrino luminosity is then given by

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }} & \approx & \displaystyle \frac{3}{8}{f}_{p\gamma }({\varepsilon }_{p}{L}_{{\varepsilon }_{p}})\\ & \simeq & 1.2\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{f}_{p\gamma ,-4}\left(\displaystyle \frac{{\varepsilon }_{p}{L}_{{\varepsilon }_{p}}}{{10}^{49.5}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right),\end{array}\end{eqnarray} \tag{ 12 }$$

which is consistent with the results of Keivani et al. (2018). The remaining fraction (i.e., 5/8) of energy should be carried by pionic γ-rays with ${\varepsilon }_{\gamma }^{{\prime} }\approx 0.1{\varepsilon }_{p}^{{\prime} }$ and secondary electrons and positrons with ${\gamma }_{e}^{{\prime} }\approx 0.05{\varepsilon }_{p}^{{\prime} }/({m}_{e}{c}^{2})\simeq 2.9\times {10}^{7}\,({\varepsilon }_{p}/6\,\mathrm{PeV})(20/\delta )$. The TeV–PeV γ-rays are attenuated inside the source, which also generate the pairs. The Bethe–Heitler process also injects high-energy pairs with ${\gamma }_{e}^{{\prime} }\approx 5\times {10}^{-4}{\varepsilon }_{p}^{{\prime} }/({m}_{e}{c}^{2})\,\simeq 2.9\times {10}^{5}\,({\varepsilon }_{p}/6\,\mathrm{PeV})(20/\delta )$ (e.g., Mastichiadis & Kirk 1995); even more energetic pairs can be produced by interactions happening far from the energy threshold of the process (e.g., Kelner & Aharonian 2008). These highly relativistic pairs quickly lose their energies via synchrotron and IC cooling. The cooling Lorentz factor is ${\gamma }_{c}^{{\prime} }\approx 2300\,B{{\prime} }_{-0.5}^{-2}{{l}_{17}^{{\prime} }}^{-1}{(1+{Y}_{\mathrm{IC}})}^{-1}$, implying that the resulting cascade spectrum lies in the fast-cooling regime. In the case of TXS 0506+056, the synchrotron peak is comparable to the IC peak, and the Compton Y parameter (Y_IC) is at most unity (Keivani et al. 2018).

The synchrotron emission from pairs injected via the Bethe–Heitler process is not always negligible in blazars, as demonstrated by Petropoulou & Mastichiadis (2015). It turns out to be important also for TXS 0506+056 during its high state (Keivani et al. 2018). The minimum synchrotron cascade flux associated with the neutrino flux at ${\varepsilon }_{\nu }$ is

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{\gamma }{L}_{{\varepsilon }_{\gamma }}{| }_{{\varepsilon }_{\mathrm{syn}}^{\mathrm{BH}}} & \approx & \displaystyle \frac{1}{2(1+{Y}_{\mathrm{IC}})}g[\beta ]{f}_{p\gamma }({\varepsilon }_{p}{L}_{{\varepsilon }_{p}})\\ & \approx & \displaystyle \frac{4g[\beta ]}{3(1+{Y}_{\mathrm{IC}})}{\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }},\end{array}\end{eqnarray} \tag{ 13 }$$

where ${\varepsilon }_{\mathrm{syn}}^{\mathrm{BH}}\simeq 6\,\mathrm{keV}B{{\prime} }_{\,-0.5}^{}{({\varepsilon }_{p}/6\mathrm{PeV})}^{2}(20/\delta )$ is the characteristic frequency of synchrotron emission by pairs from protons with ${\varepsilon }_{p}\approx 20{\varepsilon }_{\nu }$. Because of the broad distribution of pairs injected by the Bethe–Heitler processes, even if the protons are monoenergetic, (Dimitrakoudis et al. 2012), the expected synchrotron spectrum will be extending over several decades in energy (e.g., Petropoulou & Mastichiadis 2015; Petropoulou et al. 2015). Note that for sufficiently high energy pairs we expect ${Y}_{\mathrm{IC}}\ll 1$ owing to the Klein–Nishina suppression.

Similarly, for synchrotron emission from pairs injected via photomeson production and two-photon annihilation for pionic γ-rays, the synchrotron cascade flux is

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{\gamma }{L}_{{\varepsilon }_{\gamma }}{| }_{{\varepsilon }_{\mathrm{syn}}^{p\gamma }} & \approx & \displaystyle \frac{1}{2(1+{Y}_{\mathrm{IC}})}\displaystyle \frac{5}{8}{f}_{p\gamma }({\varepsilon }_{p}{L}_{{\varepsilon }_{p}})\\ & \approx & \displaystyle \frac{5}{6(1+{Y}_{\mathrm{IC}})}{\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }},\end{array}\end{eqnarray} \tag{ 14 }$$

where ${\varepsilon }_{\mathrm{syn}}^{p\gamma }\simeq 60\,\mathrm{MeV}\,B{{\prime} }_{-0.5}{({\varepsilon }_{p}/6\mathrm{PeV})}^{2}(20/\delta )$ and the contribution of pionic γ-rays is included assuming that they are converted into pairs inside the source.

In addition to the synchrotron cascade components considered above, the IC emission and subsequent regeneration processes can affect the pair-injection spectrum. Although the exact spectral shape of a cascade photon spectrum depends on details of the pair injection and possible contributions from muon and meson radiation, the resulting energy spectrum becomes approximately flat, that is, it can be expressed as ${E}_{\gamma }{F}_{{E}_{\gamma }}\propto {E}_{\gamma }^{2-\beta }$ with β ∼ 1.5–2 varying in the X-ray and γ-ray range.

For the TXS 0506+056 flare coincident with IceCube-170922A, Swift and NuSTAR measured X-rays quasi-simultaneously. A more recent sophisticated analysis gave the X-ray flux, ${E}_{\gamma }{F}_{{E}_{\gamma }}^{X}\approx 0.8\times {10}^{-12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ at E_γ ≈ 2−3 keV (Keivani et al. 2018). This leads to tight limits on the high-energy neutrino flux from the TXS 0506+056 flare. Combining Equation (13) with the observed X-ray flux, the neutrino flux in the 0.1–1 PeV range is constrained as ${E}_{\nu }{F}_{{E}_{\nu }}^{0.1-1\,\mathrm{PeV}}\lesssim {E}_{\gamma }{F}_{{E}_{\gamma }}^{X}\sim {10}^{-12}$ $\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ for all flavors, or we have ${\varepsilon }_{\nu }{L}_{{\varepsilon }_{{\nu }_{\mu }}}^{0.1-1\,\mathrm{PeV}}\lesssim {\varepsilon }_{\gamma }{L}_{{\varepsilon }_{\gamma }}^{X}/3\sim {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, where a factor of 3 comes from ν_e:ν_μ:ν_τ ≈ 1:1:1. This is fully consistent with the detailed numerical results presented in Keivani et al. (2018), and the upper-limit neutrino luminosity is much lower than the luminosity suggested for the 2017 and 2014–2015 flares of TXS 0506+056, ${\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }}(\sim 3{\varepsilon }_{\nu }{L}_{{\varepsilon }_{{\nu }_{\mu }}})\sim {10}^{46}\mbox{--}{10}^{47}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. Note that the neutrino luminosity around the peak at higher energies (see Figure 4 of Keivani et al. 2018) is slightly higher owing to g[β] and details of the cascade spectrum.

Equations (13) and (14) show that the luminosity of the synchrotron cascade components is comparable to the neutrino luminosity, as long as ${Y}_{\mathrm{IC}}\ll 1$. Thus, the cascade bound on the neutrino flux is unavoidable as long as the photomeson production occurs in a compact region such as the blazar zone. If the canonical picture of blazars based on the single-zone modeling is correct, these results allow us to predict that the 2017 and 2014–2015 neutrino flares reported by the IceCube Collaboration (IceCube Collaboration 2018) should be accompanied by X-ray emission with ${E}_{\gamma }{F}_{{E}_{\gamma }}^{X}\sim {E}_{\nu }{F}_{{E}_{\nu }}^{0.1\mbox{--}1\,\mathrm{PeV}}\sim (3\mbox{--}30)\times {10}^{-11}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$, which should be detectable by X-ray sky monitors such as Swift and MAXI.

Finally, we comment on high-energy neutrino emission through hadronuclear processes. For a relativistic jet of blazars, pp interactions such as the neutrino production mechanism cannot be globally efficient at least in the blazar zone; otherwise, the jet becomes too heavy to remain relativistic. However, for flares, one could potentially avoid this problem by invoking some mechanism to entrain dense clouds in the jets. Several possibilities have been discussed, which include clouds in the broad-line region, stellar winds, and the tidal stripping of the stellar envelope (e.g., Bednarek & Protheroe 1997; Bosch-Ramon et al. 2012; Perucho et al. 2017). In addition, massive stars could supply matter via supernova explosions or mass eruptions. Although various possibilities could potentially be considered (see Figure 2), there are several issues in the explanation for the month-scale neutrino flares from TXS 0506+056. First, there has been much evidence that the clouds in the broad-line region are mostly located around the equatorial plane, not on the jet axis (e.g., Wills & Browne 1986; Runnoe et al. 2013; Shen & Ho 2014). The profiles of the broad-line emission often show a disk-like rotation (e.g., Eracleous et al. 2009; Luo et al. 2013), and the vertical velocity component can be naturally explained by the dusty disk wind (e.g., Czerny & Hryniewicz 2011; Czerny et al. 2015). On the other hand, there are a large number of stars around the AGN core, and some of them could be on the jet axis. No matter what the origin is, the cloud mass required to achieve a high column density of $\gtrsim {10}^{24}\,{\mathrm{cm}}^{-2}$ is ${M}_{c}\gtrsim 0.4\,{M}_{\odot }\,{l}_{c,16}^{2}$ for cloud size l_c (note that the observations suggest that the comoving size of the blazar zone is $l^{\prime} \sim {10}^{16}\mbox{--}{10}^{17}$ cm). With a typical jet luminosity of BL Lacs, it is difficult for a massive cloud to become highly relativistic. If the CR generation occurs via collisions between the jet and a dense cloud at rest, the characteristic timescale is given by $(1+z){l}_{c}/{\delta }^{2}c$ $\simeq \,1.1\times {10}^{4}$ ${\rm{s}}\,{l}_{c,17}{(20/\delta )}^{2}$, which is shorter than the observed variability and duration, although the evaporation and expansion subsequent to the jet–cloud interaction may eventually lead to a longer timescale. Third, the existence of a powerful jet implies that such a cloud can be largely ionized, and strong photoelectric absorption of X-rays is unlikely to happen except for the innermost region of the cloud. For a cloud density of ${n}_{c}\simeq 2.8\times {10}^{6}\,{\mathrm{cm}}^{-3}({M}_{c}/10\,{M}_{\odot }){l}_{c,17}^{-3}$, the ionization radius is estimated to be $\sim 8\,\times {10}^{17}\,\mathrm{cm}\,{({\varepsilon }_{\gamma }{L}_{{\varepsilon }_{\gamma }}/{10}^{45}\mathrm{erg}{{\rm{s}}}^{-1})}^{1/3}\,{({\varepsilon }_{\gamma }/100\mathrm{eV})}^{-1/3}{n}_{c,7}^{-2/3}\,{{ \mathcal T }}_{5}^{1/4}$, where ${ \mathcal T }$ is the cloud temperature. Except for extreme cases with N_H > Σ_T⁻¹, we expect that X-rays escape and the cascade limit is important even for pp models, and we have

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{\gamma }{L}_{{\varepsilon }_{\gamma }}{| }_{{\varepsilon }_{\mathrm{syn}}^{{pp}}} & \approx & \displaystyle \frac{1}{2(1+{Y}_{\mathrm{IC}})}\displaystyle \frac{1}{2}{f}_{{pp}}({\varepsilon }_{p}{L}_{{\varepsilon }_{p}})\\ & \approx & \displaystyle \frac{1}{2(1+{Y}_{\mathrm{IC}})}{\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }}.\end{array}\end{eqnarray} \tag{ 15 }$$

For a spectrum of CRs with s_cr ∼ 2, the above relation holds at energies from ${\varepsilon }_{\mathrm{syn}}^{c}\simeq 0.39\,\mathrm{eV}\,B{{\prime} }_{-0.5}^{-3}l{{\prime} }_{17}^{-2}{(1+{Y}_{\mathrm{IC}})}^{-2}(\delta /20)$ to ${\varepsilon }_{\mathrm{syn}}^{{pp},\max }\simeq 60\,\mathrm{MeV}\,{B}_{-0.5}^{{\prime} }{({\varepsilon }_{p}/6\mathrm{PeV})}^{2}(20/\delta )$.

Zoom In Zoom Out Reset image size

An electromagnetic cascade is a consequence of energy conservation, so the cascade limit discussed in the previous subsection exists for not only single-zone but also multizone models.

One of the possibilities is that neutrinos are mainly produced around the base of the jet or even in the vicinity of a supermassive black hole (i.e., AGN core models; e.g., Stecker et al. 1991; Bednarek & Protheroe 1999; Becker & Biermann 2009; Stecker 2013; Kimura et al. 2015). This includes high-energy neutrino production caused by pp interactions with dense clouds as discussed above. Because of the higher compactness of the system, broadband cascade emission should accompany the neutrino emission. In such an inner region, magnetic fields should also be much stronger, which could affect the resulting cascade spectrum. Importantly, not only protons but also electrons are accelerated, and both emission components should be taken into account. Detailed studies are beyond the scope of this work (A. Mastichiadis & M. Petropoulou 2018 in preparation).

Here we consider the CR-induced neutral beam model (see Figure 2), which can avoid the cascade constraints. We consider interactions between beamed CRs escaping from the blazar zone and an external radiation field. Although we do not specify the origin of the external photons, this setup is analogous to the one considered in Murase et al. (2014) and Dermer et al. (2014). In this sense, this two-zone model can be regarded as a natural extension of the standard leptonic scenario. In particular, we assume that escaping CRs are neutrons that can be produced via the photodisintegration of nuclei in the blazar zone. The γ-ray signatures produced by CR-induced neutral beams were previously studied in Murase (2012) and Dermer et al. (2012) (see also Essey et al. 2010, 2011; Murase et al. 2012b, for related discussion about intergalactic cascades). As we show below, the cascade emission can be largely diminished via the isotropization of relativistic electrons and positrons.

Following Dermer et al. (2014), let us assume that CRs are accelerated via the second-order Fermi acceleration mechanism. The maximum energy accelerated in the blazar zone can be ε_cr/Z ∼ 1−10 PeV, where ε_cr is the CR ion energy and Z is the nuclear charge. The CR acceleration zone can be the γ-ray emission site or inner regions of the blazar zone, and disintegrated nuclei are accompanied by not only protons but also neutrons (e.g., Murase & Beacom 2010; Rodrigues et al. 2018). The protons may lose energy via adiabatic losses during the confinement in the blazar zone, while neutrons can escape. The neutron luminosity is given by ${\varepsilon }_{n}{L}_{{\varepsilon }_{n}}\,\approx (1/2){f}_{A\gamma }({\varepsilon }_{\mathrm{cr}}{L}_{{\varepsilon }_{\mathrm{cr}}})$, where f_Aγ is the effective optical depth to the photodisintegration process (Murase & Beacom 2010) and ${\varepsilon }_{\mathrm{cr}}{L}_{{\varepsilon }_{\mathrm{cr}}}$ is the CR ion luminosity. For neutron production in the γ-ray emission region, we have ${f}_{A\gamma }\sim 0.1{({\varepsilon }_{\gamma }{L}_{{\varepsilon }_{\gamma }}/{10}^{46}\mathrm{erg}{{\rm{s}}}^{-1})(\delta /20)}^{-3}l{{\prime} }_{17}^{-1}{({\varepsilon }_{\gamma }/1\mathrm{eV})}^{-1}$ ${({\varepsilon }_{\mathrm{cr}}^{{\prime} }/{\tilde{\varepsilon }}_{A\gamma ,\mathrm{syn}}^{{\prime} })}^{\beta -1}$, where ${\tilde{\varepsilon }}_{A\gamma ,\mathrm{syn}}^{{\prime} }=0.5{m}_{A}{c}^{2}{\bar{\varepsilon }}_{\mathrm{GDR}}/{\varepsilon }_{\mathrm{syn}}^{{\prime} }$ and ${\bar{\varepsilon }}_{\mathrm{GDR}}\sim 20\mbox{--}30\,\mathrm{MeV}$. This also implies that the neutron emission may predominantly come from smaller dissipation radii at which efficient photodisintegration (i.e., ${f}_{A\gamma }\gtrsim 1$) occurs.

In single-zone models, the neutrino flares of TXS 0506+056 require unpleasantly large CR luminosities (Keivani et al. 2018). This problem still exists at some level even in multizone models, although it can be alleviated in the CR beam model in the sense that the meson production efficiency is enhanced by additional target photons or nucleons. One should keep in mind that observations and modeling of radio galaxies (based on larger-scale jets than the blazar zone) have shown that the absolute jet power averaged over the lifetime of the AGN jet is ${P}_{j}\lesssim {10}^{45}\mbox{--}{10}^{46}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ for Fanaroff–Riley I galaxies (Cavagnolo et al. 2010; Godfrey & Shabala 2013), which are believed to be off-axis counterparts of BL Lacs. For the supermassive black hole mass M_BH, the Eddington luminosity⁸ is ${L}_{\mathrm{Edd}}\simeq 1.3\times {10}^{47}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,({M}_{\mathrm{BH}}/{10}^{9}\,{M}_{\odot })$. Ghisellini et al. (2014) showed that the absolute jet power of blazars may exceed the accretion luminosity, and our study implies that the flaring jet power is larger than the time-averaged one by ${b}_{\mathrm{fl}}/{f}_{\mathrm{fl}}\sim 3\mbox{--}10$. The isotropic-equivalent CR luminosity during the flaring phase can then be written as ${L}_{\mathrm{cr}}^{\mathrm{fl}}$ $\approx \,(2/{\theta }_{\mathrm{beam}}^{2})$ ${\epsilon }_{\mathrm{cr}}{P}_{j}({b}_{\mathrm{fl}}/{f}_{\mathrm{fl}})$ $\simeq 6.0\times {10}^{49}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ ${({\theta }_{\mathrm{beam}}/0.05)}^{-2}$ $({\epsilon }_{\mathrm{cr}}/0.2)$ $({b}_{\mathrm{fl}}{f}_{\mathrm{fl}}^{-1}/10)$ $({P}_{j}/0.3{L}_{\mathrm{Edd}})({M}_{\mathrm{BH}}/{10}^{9}\,{M}_{\odot })$, where _cr is the energy fraction carried by CR ions and θ_bream is the opening angle of the CR beam. The neutron luminosity during the flaring phase results in ${L}_{n}^{\mathrm{fl}}$ ≃ 3.0 × ${10}^{49}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ ${f}_{A\gamma }$ ${({\theta }_{\mathrm{beam}}/0.05)}^{-2}$ $({\epsilon }_{\mathrm{cr}}/0.2)({b}_{\mathrm{fl}}{f}_{\mathrm{fl}}^{-1}/10)$ $({P}_{j}/0.3{L}_{\mathrm{Edd}})({M}_{\mathrm{BH}}/{10}^{9}\,{M}_{\odot })$.

The neutrons that leave the CR acceleration zone propagate along the jet and may interact with external radiation fields that could exist on larger scales or perhaps a dense cloud. For LSP and ISP objects like TXS 0506+056, it is possible to invoke such a setup. For example, if the jet is structured, nonthermal photons can be provided by the sheath region. Moreover, a fraction of UV and X-ray emission from the accretion disk can be scattered by clumps of matter that may be present at outer radii. In addition, there could be high-velocity clumps such as the broad-line region, although they are usually seen in FSRQs. Note that the neutrons with γ_n ∼ 10⁷–10⁸ can travel ∼0.1–1 kpc.

Interestingly, the detailed modeling of the SED of TXS 0506+056 (Keivani et al. 2018) already suggested that such an external radiation field is necessary to explain the X-ray and γ-ray spectrum. If this is the case, it is natural for escaping CRs to keep interacting with the ambient photons, leading to the production of more neutrinos.

As a toy model, we assume that the decelerated jet or slower jet of the sheath region provides soft photons with a luminosity of ${L}_{\mathrm{ext}}\sim 3\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and the characteristic energy at ${\varepsilon }_{\mathrm{ext}}\sim 10\,\mathrm{eV}$, over a length scale of ${R}_{\mathrm{ext}}\,\sim 3\times {10}^{19}$ cm. The external radiation energy density is ${U}_{\mathrm{ext}}\approx 3{L}_{\mathrm{ext}}/(4\pi {R}_{\mathrm{ext}}^{2}c)\sim 3\times {10}^{-5}\,\mathrm{erg}\,{\mathrm{cm}}^{-3}$, which is consistent with the parameters used in Keivani et al. (2018). Noting ${\hat{\sigma }}_{n\gamma }\approx {\hat{\sigma }}_{p\gamma }$, the photomeson production efficiency is ${f}_{n\gamma }$ $\approx \,[{\eta }_{n\gamma }{\hat{\sigma }}_{n\gamma }3{L}_{\mathrm{ext}}/(4\pi {R}_{\mathrm{ext}}c{\varepsilon }_{\mathrm{ext}})]$ ${({\varepsilon }_{n}/{\tilde{\varepsilon }}_{n\gamma ,\mathrm{ext}})}^{\beta -1}$ $\sim 3\times {10}^{-3}\,{\eta }_{n\gamma }$ ${L}_{\mathrm{ext},45.5}{R}_{\mathrm{ext},19.5}^{-1}$ ${({\varepsilon }_{\mathrm{ext}}/10\mathrm{eV})}^{-1}$ ${({\varepsilon }_{n}/{\tilde{\varepsilon }}_{n\gamma ,\mathrm{ext}})}^{\beta -1}$ (where η_nγ = η_pγ and ${\tilde{\varepsilon }}_{n\gamma ,\mathrm{ext}}$ = 0.5 ${m}_{n}{c}^{2}{\bar{\varepsilon }}_{{\rm{\Delta }}}/{\varepsilon }_{\mathrm{ext}}$ $\approx \,{\tilde{\varepsilon }}_{p\gamma ,\mathrm{ext}}$), which can be larger by a factor of $\sim {R}_{\mathrm{ext}}/({\rm{\Gamma }}l^{\prime} )$ compared to f_pγ in the blazar zone (see Equation (9)). Alternatively, instead of nγ interactions, one could expect np interactions between the neutral beam and a dense cloud in the line of sight. Although the existence of such a cloud is speculative, with a clump like a giant molecular cloud with a mass of ${M}_{c}\sim ({10}^{5}\mbox{--}{10}^{6})\,{M}_{\odot }$ and a size of ${l}_{c}\sim {10}^{19}$ cm, the effective np optical depth can be as high as ${f}_{{np}}\,\simeq 0.011\,({M}_{{\rm{c}}}/{10}^{5}\,{M}_{\odot }){l}_{c,19}^{-2}$.

The all-flavor neutrino luminosity is

$$\begin{eqnarray}&&{\varepsilon }_{\nu }{L}_{{\varepsilon }_{\nu }}\simeq 3.8\times {10}^{46}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,\left(\displaystyle \frac{2K}{1+K}\right)\left(\displaystyle \frac{{f}_{n\gamma /{np}}{\varepsilon }_{n}{L}_{{\varepsilon }_{n}}}{{10}^{47}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right),\end{eqnarray} \tag{ 16 }$$

where K = 1 and K = 2 for np and nγ interactions, respectively. In the nγ case, the neutrino energy corresponding to the neutron energy with ${\tilde{\varepsilon }}_{n\gamma ,\mathrm{ext}}$ is

$$\begin{eqnarray}&&{\varepsilon }_{\nu }\approx 0.05{\tilde{\varepsilon }}_{n\gamma ,\mathrm{ext}}\sim 0.7\,\mathrm{PeV}\,{({\varepsilon }_{\mathrm{ext}}/10\mathrm{eV})}^{-1},\end{eqnarray} \tag{ 17 }$$

which could explain both of the 2017 and 2014–2015 neutrino flares of TXS 0506+056, e.g., with f_Aγ ∼ 0.1 and f_Aγ ≳ 1, respectively. The duration of neutrino emission is comparable to the "lifetime" of the CR-induced beam, which is being determined by the duration of particle energization in the CR acceleration zone corresponding to the observed t_dur.

The key point of the CR-induced neutral beam model considered here is that the cascade signature except for γ-rays in the very high energy range can be largely diminished because of the isotropization of the relativistic pairs in the larger-scale jet or other magnetized environments. Let us assume that the magnetic field in the main scale of neutrino production is as small as B_ext ∼ 0.1–10 mG (typical for large-scale jets). The deflection of pairs occurs before they cool via synchrotron and IC losses. Following Murase (2012), the deflection angle during the radiation cooling time is given by

$$\begin{eqnarray}&&{\theta }_{\mathrm{def}}\approx \sqrt{\displaystyle \frac{2}{3}}\displaystyle \frac{{{ct}}_{\mathrm{syn}}}{{r}_{L}}\simeq 3.5\,{\gamma }_{e,9}^{-2}{B}_{\mathrm{ext},-2.5}^{-1},\end{eqnarray} \tag{ 18 }$$

where t_syn is the synchrotron cooling time, r_L is the Larmor radius, and the Lorentz factor of pairs in the black hole rest frame is ${\gamma }_{e}\approx 0.05{\varepsilon }_{n}/({m}_{e}{c}^{2})\simeq 5.8\times {10}^{8}\,({\varepsilon }_{n}/6\,\mathrm{PeV})$. Only synchrotron cooling is considered, because the IC cooling is suppressed for such high-energy electrons and positrons. The above equation implies that the pairs lose their energy after they become isotropized. Note that the X-ray emission is doubly suppressed, because the deflection is larger than the jet opening angle, i.e., ${\theta }_{\mathrm{def}}\gg {\theta }_{j}$, and the time spread in the cascade emission is longer than the intrinsic flare duration, i.e., ${R}_{\mathrm{ext}}/c\gg {t}_{\mathrm{dur}}$. Also, the CR neutrons do not produce extra pairs via the Bethe–Heitler process. In addition, in the np scenario, X-rays could also be isotropized via the Compton scattering. Thus, X-rays from the neutral beam can be highly suppressed if the CR energy is lower than UHE energies, and the associated hadronic cascade signatures can be masked by the GeV γ-ray emission from the blazar zone. Note that the UHE neutron beam, if it exists, can also produce beamed cascade emission at $\gtrsim 0.1$ TeV energies. The resulting synchrotron pair-echo emission (Dermer et al. 2012; Murase 2012) is expected to have some time delay. We also remark that the very high energy flare seen by MAGIC was detected about 10 days later after IceCube-170922A.