High-Energy Neutrinos from Blazar Flares and Implications of TXS 0506+056

Motivated by the observation of a $>290$ TeV muon neutrino by IceCube, coincident with a $\sim$6 month-long $\gamma$-ray flare of the blazar TXS 0506+056, and an archival search which revealed $13 \pm 5$ further, lower-energy neutrinos in the direction of the source in 2014-2015, we discuss the likely contribution of blazars to the diffuse high-energy neutrino intensity, the implications for neutrino emission from TXS 0506+056 based on multi-wavelength observations of the source, and a multi-zone model that allows for sufficient neutrino emission so as to reconcile the multi-wavelength cascade constraints with the neutrino emission seen by IceCube in the direction of TXS 0506+056.


Introduction
The IceCube Collaboration recently reported the observation of a > 290 TeV muon neutrino, IceCube-170922A, coincident with a ∼6 month-long γ-ray flare of the blazar TXS 0506+056 [1], at redshift z ≈ 0.3365 [2]. The neutrino detection prompted multi-wavelength electromagnetic follow-up of the source, and the blazar flare was observed with several instruments at energies up to > 100 GeV [3]. The correlation of the neutrino with the flare of TXS 0506+056 is inconsistent with arising by chance at the ∼ 3σ level. An archival search revealed 13 ± 5 further, lower-energy neutrinos in the direction of TXS 0506+056 during a 6-month period in 2014-2015 [4]. These events were not accompanied by a γ-ray flare. Such an accumulation of events is inconsistent with arising from a background fluctuation at the 3.5σ level. Motivated by these observations in [5], we considered the implications of the possible neutrino-blazar flare association. Here, we summarise and update some of these results, in light of Ice-Cube and other related analysis updates on the topic. In section 2.1 we review existing constraints to blazars as dominant sources of the diffuse neutrino intensity seen by IceCube. In section 2.2 we determine the duty cycle of γray flares of several Fermi-detected blazars and compute the resulting neutrino enhancement factor. In section 2.3 we present the constraints imposed to single-zone models of neutrino emission from the contemporaneous broadband spectral energy distribution (SED) of the source, and from X-ray and γ-ray observations of putative neutrino sources in general. Section 3 presents a multi-zone model of neutrino emission based on cosmic-ray induced neutral beams, which allows for enhanced neutrino production with respect to standard one-zone models, without violating the multi-wavelength constraints from the SED of TXS 0506+056. Conclusions are presented in section 4.
2 Blazar contribution to the diffuse neutrino flux

Clustering constraints
The absence of high-energy multiplets in the IceCube data can be used to constrain the number density of sources contributing to the diffuse neutrino background [6,7] (see also earlier work by [8][9][10]). The IceCube diffuse flux, E 2 ν Φ ν , gives an estimate of the neutrino production rate which is at the level of ∼ (c · t H · n eff · ε ν L ave ε ν )/∆Ω for standard candle sources with effective number density, n eff and time-averaged luminosity ε ν L ave ε ν [11]. Here, t H is the Hubble time, c the speed of light, and ∆Ω, the solid angle covered by the detector. Here, and throughout, primed quantities correspond to comoving frame quantities. Unprimed scripted, ε X , corresponds to cosmic-rest frame energy, and capital unscripted, E X , to observer frame energy.
As shown in [6], one can express the number of highenergy neutrino doublets, N m≥2 , as, where, d n=1 , is the distance within which the number of expected events from a single source is equal to 1, n eff 0 = n eff [z = 0] the local source number density, and q L , is a luminosity dependent function that depends on the redshift evolution of the source population and approaches unity for low-luminosities. For example, q L = 0.94 and q L = 2.0 for redshift evolution of the form n s (z) ∝ (1 + z) 0 and n s (z) ∝ (1 + z) 3 respectively, for luminosity corresponding to d n=1 /(c/H 0 ) = 0.1, where H 0 is the Hubble constant. Using eq. 1, and expressing d n=1 , as [6], where F lim ∼ 3 × 10 −10 GeV cm −2 s −1 [12] is the Ice-Cube point-source sensitivity (90% CL) for an E −2 neutrino spectrum, one can derive a luminosity dependent upper limit on n eff , in the absence of doublets in the IceCube data, i.e. imposing N m≥2 < 1, where F lim = 10 −9.5 F lim,−9.5 GeV cm −2 s −1 . Using eq. 3 we can derive an upper limit to the contribution of a source population with number density n 0,eff to the IceCube allflavor neutrino flux, where the function ξ z , has been introduced to parametrise the redshift evolution of the neutrino emissivity of the source population; ξ z = 0.7 for the γ-ray luminosity density evolution of BL Lacs, ξ z = 8 for that of FSRQs [13], and ξ z = 3 for the X-ray luminosity density evolution of AGNs [14]. The above limits give conservative constraints on flaring neutrino sources, as shown in [5]. In general, the multiplet limits apply to transient sources, and therefore in the present context to blazar flares. The rate density is constrained as, where we substituted N fl ≈ T obs /∆T fl , with N fl , the number of flaring periods, and T obs , the observation duration, and φ lim = 0.1 φ lim,−1 GeV cm −2 is the muon neutrino fluence sensitivity, φ lim ∼ 0.04 GeV cm −2 , which can be calculated by the publicly available effective area. The limit of equation 2.1 applies to the 2017, as well as the 2014-15 flares. Even in the latter case, where a large number of low energy neutrinos were observed, the number of > 50 TeV neutrino multiplets, N m≥2 ≤ 1 and the limit of equation holds. For rare transients (no repeating bursts in the observation time), we have ρ eff 0 = n eff 0 /∆T fl , with ∆T fl , the interval between flares. The average number of sources in flaring state is n eff 0 ∆T dur /∆T fl . Based on the above arguments, we conclude that the absence of doublets and higher order multiplets in the Ice-Cube data, disfavours BL Lacs as the dominant source of the diffuse IceCube flux. However, the constraints could be relaxed by several reasons [5]. For example, rapidly evolving FSRQs (including MeV blazars) could make a substantial contribution only in the PeV range (see also, e.g., [15,16]).
The total blazar contribution to the diffuse neutrino background has further been constrained in other analyses of event clustering and autocorrelation [17,18]. In addition, the absence of extremely-high energy neutrinos (>5 PeV) in diffuse searches so far [19], constrains optimistic models of blazar neutrino emission [e.g. 20,21]. The reason is that typically blazar neutrino models peak at > PeV energy, as the target photon density is larger at lower energies in blazars. Finally, the contribution of blazars to the IceCube flux has been constrained by cross-correlation and stacking analyses [18,22]. The latest update of this analysis constrains the contribution of γ-ray bright blazar emission to ≤ 27% of the IceCube flux.

Contribution from blazar flares: Duty cycle and enhancement factor
In order to assess the fraction of neutrinos that could be produced during flaring states we studied the γ-ray lightcurves of a sample of blazars from the Fermi Large Area Telescope, FAVA catalogue [23]. In order to investigate the fraction of time spent in a flaring state we introduce f fl , the flare duty factor, defined as, where N tot is the total number of time bins, L th is the threshold luminosity beyond which we consider the source to be flaring, and dN/dL is the distribution of luminosity states. In addition, we quantify the fraction of energy emitted in the flaring state, b fl , defined as, where the average luminosity is given by L ave = (1/N tot ) dL L(dN/dL) and the average flaring luminosity is L fl = (b fl / f fl )L ave . In the flare-dominated limit, where, b fl ≈ 1, L ave approaches L ave ≈ f fl L fl . We determined b fl and f fl for a sample of FAVA sources by using the weeklybinned publicly available photon counts for each source.
We thus obtain the photon-count distribution, dN/dN γ , which is proportional to the γ-ray luminosity distribution, dN/dL γ , of the source in this time and energy bin, as long as the photon index doesn't change. For the rest of this discussion we assume dN/dN γ ∼ dN/dL γ .  Figure 1. Histogram of the number of photons, N γ , detected per week by the FAVA analysis in the direction of TXS 0506+056, OJ 287, PKS 0301-243, and S4 0954+65 in the high-energy bin (800 MeV−300 GeV). The photon distribution is modelled as a power law with spectral index α convolved with a Poissonian distribution (solid red line). Table 1. Flare duty factor, f fl , and fraction of energy released during flares, b fl , for a sample of sources as derived from the FAVA analysis. We report values for the low-energy, LE, 100 − 800 MeV, and high-energy, HE, 800 MeV−300 GeV, bins.
The duty factors quoted are for flux enhancement ≥ 5σ according to the FAVA definition. The γ-ray luminosity of the sources, L γ (in units of erg s −1 ), is derived from the 3FGL in the 1 − 100 GeV energy range. Rounded values of the power-law index, α, are also shown.
Name , and obtain f fl ≈ 0.02 − 0.1 for TXS 0506+056. The corresponding fraction of emitted photons is b fl ∼ 0.1 1, implying that the bulk of the γ-ray emission of the studied sources comes from quiescent periods. We further investigated the frequency and duty cycle of flares by fitting the photon count distribution of each FAVA source with a power-law, dN/dN γ ∼ dN/dL γ ∝ L −α γ . The number of detected photons per time interval is then given by a convolution of this power law with a Poissonian distribution. Example fits for TXS 0506+056 and some of the other sources analysed are shown in Fig. 1. For TXS 0506+056, α 3. For the other sources analysed we found α ∼ 2 − 4. The above finding that the duty factor of flares and corresponding fraction of emission during flaring is subdominant might lead to the conclusion that the majority of neutrinos are produced during quiescent states in blazars. However, even though for the sources that we studied highly-significant flaring states are achieved for a small fraction of time, during which a non-dominant part of the observed γ-ray emission is released this may not be the case for the neutrino emission for at least two reasons, which we discuss below in turn.
Firstly, neutrinos can be preferentially produced during γ-ray flares as typically, in models where the majority of γ-rays are leptonic in origin, L ν ∝ L γ rad , with γ ∼ 1.5 − 2.0. This was demonstrated, for example, in [6,24,25] and references therein. Therefore, which implies that neutrino production from flaring states can dominate the total neutrino emission of such sources even if flares don't dominate the radiative output. This is, for example, the case for sources with γ ∼ 2.0 and α ≤ 3. As shown in fig. 1 the FAVA data of TXS 0506+056 and some of the other sources we analysed are consistent with these sources being in the latter category. Second, neutrino emission can be significantly enhanced during flaring, if the proton spectrum changes to a harder spectral index. In order to demonstrate this, we consider as an example a source where protons are produced with a power-law spectrum with index s l = 2.3 during steady state emission, and a proton spectrum which becomes harder during the flare, with index, s fl = 1.8. For such a scenario we find that the enhancement factor c[ε ν ] is given by, If we further assume that the neutrino energy is ε ν = 0.05ε flp,max = 0.1 PeV, and the minimum proton energy during steady state is ε fl l = 10 GeV, we find c[ε ν ] ∼ 30 f fl pγ / f l pγ , which further demonstrates that it is possible that neutrino emission during electromagnetic flares dominates the total neutrino output of the source.
We conclude that the contribution of flaring blazars like those observed from TXS 0506+056 to the diffuse background can be limited as, Note that the general contribution cannot exceed the limit of equation 4.
A further limit is imposed on the maximum neutrino luminosity from the requirement that the synchrotron cascade flux produced by electron-positron pairs, injected by Bethe-Heitler interactions, should not exceed the observed X-ray flux.
where, ε BH syn ≈ 6 keVB 0.5 G (ε p /6 PeV) 2 (20/δ), with δ the Doppler factor of the relativistic motion of the emitting region, B the magnetic field strength as measured by a comoving observer, and Y IC the Compton dominance parameter which is at most 1 for TXS 0506+056 [27]. For the flaring spectrum of TXS 0506+056 from the analysis of [27] this corresponds to ε γ L X ε γ ≤ 3 × 10 44 erg/s. This imposes a constraint on the maximum muon neutrino luminosity at the level of, ε ν L 0.1−1 PeV ε νµ ε γ L X ε γ /3 ∼ 10 44 erg s −1 , where the factor of 3 accounts for going from all-flavour to single-flavour neutrino luminosity. Eq. 18 is plotted as the mid-grey exclusion region in fig. 2 assuming Y IC = 1 which gives the most optimistic estimate for the maximum neutrino luminosity of TXS 0506+056.
In addition, a synchrotron cascade flux component from electron-positron pairs injected from photo-meson production and electron-positron pair production from hadronic γ rays is unavoidable in the single-zone model. The minimum cascade flux from these processes is, where ε pγ syn 60 MeV (B /0.3 G))(ε p /6 PeV) 2 (20/δ). The constraint of eq. 19 is shown as a further exclusion region in fig. 2, where we conservatively assume that ε γ L ε γ | ε pγ syn ∼60 MeV ≤ 5×10 45 ergs −1 . As evident, for the 2017 flare of TXS 0506+056, the X-ray cascade constraint of eq. 18 is stronger. Fig. 3 shows the cascade constraints for the archival, 2014-15-neutrino flare. Here, we plot exclusion regions for ε γ L ε γ | ε pγ syn ∼60 MeV ≤ 2 × 10 45 ergs −1 , and ε γ L ε γ | ε pγ syn ∼10 keV ≤ 2 × 10 46 ergs −1 , based on the analyses of [30,31]. In this case the MeV data provide a stronger constraint. In deriving constraints in the MeV region for fig. 3 and 3 we have assumed that there is no additional, Bethe-Heitler bump in the SED (see e.g fig. 1 of [30]). If such a component did exist, our upper limit for the 2014-15 flare would be too conservative.
Equations 18 and 19 show that the luminosity of the synchrotron cascade is comparable to the neutrino luminosity in single-zone models. Thus, the 2017 and 2014-15 flares reported by the IceCube Collaboration should be accompanied by X-ray emission with which should be detectable by X-ray monitoring instruments such as Swift and MAXI. . Schematic representation of the cosmic-ray beam model for high-energy neutrino production. While the neutrino emission is highly beamed the associated cascade emission in the X-ray range is isotropised.

Multi-zone model and cosmic-ray induced neutral beams
The electromagnetic cascade emission discussed in the previous section is a consequence of energy conservation, and hence expected not only in single-zone models, but also in multi-zone models. It is, therefore, crucial to take into account cascades inside the source.
In this section we discuss the cosmic-ray induced neutral beam model, which was previously studied in [32] and [33] as a possible way to explain the 2017 and 2014-15 flare of TXS 0506+056 in a "common" framework. It has the appealing feature that it can avoid the cascade constraints, as isotropisation of high energy electrons and positrons is expected to take place if such conditions are found in the source environment.
For the cosmic-ray induced neutral beam model to be in operation, nuclei must be accelerated in the emitting region of the blazar. Subsequently, the following steps take place: 1. Neutrons are produced via the photo-disintegration of nuclei in the cosmic-ray acceleration region. This region can be the vicinity of the central black hole.
2. Nuclei and protons remain confined and eventually cool via adiabatic losses while neutrons escape the cosmic-ray acceleration zone.
3. The neutrons keep to interact with external radiation fields (e.g., scattered accretion disk emission and non-thermal emission from the sheath) or dense clouds that could exist at larger radii producing neutrinos.
4. The relativistic pairs produced in nγ or np interactions get isotropised by the magnetic fields in the jet or the surrounding medium. The emitted synchrotron radiation is not boosted, thus suppressing the electromagnetic cascade that is otherwise expected.
A schematic representation of the model is shown in fig. 4. The electromagnetic cascade is doubly suppressed since the electron-positron pairs can be isotropised in the larger scale jet, which also causes a spread in time, and in addition, neutrons don't undergo Bethe-Heitler interactions. The expected neutrino flux in this model is, ε ν L ε ν ≈ 3.8 × 10 46 erg/s 2K 1 + K f nγ/np ε n L ε n 10 47 erg/s , (21) where, f nγ/np is the neutrino production efficiency in the interactions of neutrons with photons/hadrons respectively, K = 1 for np and K = 2 for nγ interactions. This model can therefore explain the energetics of both the 2017 and 2014-15 neutrino flares of TXS 0506+056 as long as the neutrino production efficiency is > 0.1 and ∼ 1, respectively. The neutron luminosity quoted in eq. 21 can be produced via photo-nuclear interactions in the blazar zone as detailed in [5].

Conclusions
In light of the recent finding of a high-energy neutrino in the direction of the flaring TXS 0506+056, we presented the summary of constraints on blazars as sources of the diffuse neutrino flux seen by IceCube, and specifically the constraints on TXS 0506+056, and by extension, single powerful blazars as neutrino sources. Our results reiterate that it is difficult to reconcile blazars with all existing neutrino observations as the dominant source population. We quantified the maximum neutrino flux that can be produced by individual blazars, and showed that this is limited by the requirement that the cascade flux produced by photons co-produced in interactions of hadrons inside the source, should not exceed multi-wavelength contemporaneous observations. We presented a multi-zone model based on the concept of cosmic-ray induced neutral beams, which can escape the X-ray and γ-ray constraints that otherwise limit the maximum attainable neutrino luminosity from TXS 0506+056.