# Strong constraints on hadronic models of blazar activity from Fermi and IceCube stacking analysis

Andrii Neronov, Dmitri V. Semikoz, Ksenia Ptitsyna

► To cite this version:

Andrii Neronov, Dmitri V. Semikoz, Ksenia Ptitsyna. Strong constraints on hadronic models of blazar activity from Fermi and IceCube stacking analysis. Astronomy & Astrophysics - A&A, 2017, 603, pp.A135. 10.1051/0004-6361/201630098. hal-01584681

HAL Id: hal-01584681

https://hal.science/hal-01584681v1

Submitted on 27 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Strong constraints on hadronic models of blazar activity from Fermi and IceCube stacking analysis

Andrii Neronov<fnref n="1" />, Dmitri V. Semikoz<fnref n="2" />,<fnref n="3" />, and Ksenia Ptitsyna<fnref n="4" />

Received 19 November 2016 / Accepted 4 June 2017

## ABSTRACT

*Context.* High-energy emission from blazars is produced by electrons that are either accelerated directly, under the assumption of leptonic models of blazar activity, or produced in interactions of accelerated protons with matter and radiation fields, under the assumption of hadronic models. The hadronic models predict that γ-ray emission is accompanied by neutrino emission with comparable total energy.

*Aims.* We derive constraints on the hadronic models of activity of blazars imposed by non-detection of neutrino flux from a population of γ-ray emitting blazars.

*Methods.* We stacked the γ-ray and muon neutrino flux from 749 blazars situated in the declination strip above -5°

*Results.* Non-detection of neutrino flux from the stacked blazar sample rules out the proton-induced cascade models in which the high-energy emission is powered by interactions of a shock-accelerated proton beam in the active galactic nucleus (AGN) jet with the ambient matter or radiation field of the black hole accretion disk. The result also remains valid for the case of interactions in the scattered radiation field in the broad line region. The IceCube constraint could be avoided if the spectrum of accelerated protons is sharply peaking in the ultra-high-energy cosmic ray range, as in the models of acceleration in the magnetic reconnection regions or in the vacuum gaps of black hole magnetospheres. Models based on these acceleration mechanisms are only consistent with the data if characteristic energies of accelerated protons are higher than $10^{19}$ eV. The constraint could also be avoided if the hadronic emission component of γ-ray flux is largely sub-dominant compared to the leptonic component and/or if it only appears during flaring activity episodes, which provide negligible contribution to the time-averaged source flux.

**Key words.** galaxies: active – gamma rays: galaxies – neutrinos

## 1. Introduction

Supermassive black holes in some 10% of active galactic nuclei (AGNs) accelerate particles and produce jets that could be occasionally aligned along the line of sight (Urry & Padovani 1995). In this case the AGN appears as a blazar that is detectable in the γ-ray band. The details of the mechanisms of jet generation and particle acceleration in AGN are not clear. High-energy electrons producing the observed synchrotron and inverse Compton emission can either be accelerated directly, as in leptonic models of AGN activity, or can be generated in interactions of high-energy protons and nuclei, as in hadronic models.

Different types of hadronic models correspond to different possible mechanisms of acceleration and interactions of high-energy protons and nuclei. The proton-induced cascade model(s) ascribe the γ-ray emission to the inverse Compton emission from electromagnetic cascade initiated by the interactions of high-energy protons during their propagation through the radiation field of the AGN central engine and jet (Mannheim & Biermann 1989, 1992; Mannheim 1993). Otherwise, the γ-ray component of the spectrum could be dominated by the synchrotron radiation from the highest energy protons, as suggested

in the proton synchrotron models (Mücke & Protheroe 2001; Aharonian 2002; Mücke et al. 2003).

The proton induced cascade (PIC) models are subdivided onto two sub-types, depending on the nature of interactions of the high-energy protons. In one case, the cascade can be initiated by interactions of high-energy protons with dense radiation field created by the AGN accretion disk or the jet (Mannheim & Biermann 1989, 1992; Begelman et al. 1990; Mannheim 1993; Halzen & Zas 1997; Neronov & Semikoz 2002; Kalashev et al. 2015). Otherwise, the cascade can be initiated by interactions of protons with low-energy protons from the ambient medium, for example from the accretion flow or from the interstellar environment of the black hole (Eichler 1979; Neronov et al. 2008; Neronov & Ribordy 2009a).

The details of cascade development depend on the spectrum of primary protons. The most common assumption is that the acceleration site is a relativistic shock in the AGN jet. In this case, the proton spectrum is expected to be a cut-off power-law $dN_p/dE \propto E^{-\Gamma} \exp(-E/E_{\text{cut}})$ with the slope $\Gamma \approx 2$ and cut-off energy $E_{\text{cut}}$, which could reach the ultra-high-energy cosmic ray (UHECR) energy range. This statement is somewhat uncertain because of the uncertainty of the efficiency of geometry of magnetic field geometry in

<footnote n="1">¹ Astronomy Department, University of Geneva, Ch. d'Ecogia 16, 1290 Versoix, Switzerland
e-mail: andrii.neronov@unige.ch</footnote>

<footnote n="2">² APC, Universite Paris Diderot, CNRS/IN2P3, CEA/IRFU, Observatoire de Paris, Sorbonne Paris Cite, 119, 75205 Paris, France</footnote>

<footnote n="3">³ National Research Nuclear University MEPHI (Moscow Engineering Physics Institute), Kashirskoe highway 31, 115409 Moscow, Russia</footnote>

<footnote n="4">⁴ Institute for Nuclear Research of the Russian Academy of Sciences, 60th October Anniversary Prospect 7a, 117312 Moscow, Russia</footnote>
the AGN jet shocks, such as turbulence properties and orientation of the ordered field component with respect to the shock normal (Bednarz & Ostrowski 1998; Lemoine et al. 2006; Pelletier et al. 2009). An alternative possibility for particle acceleration is acceleration in magnetic reconnection regions or events in the AGN accretion flow (Lesch & Birk 1997; Romanova & Lovelace 1992) or in the jet, and acceleration in the vacuum gaps in the black hole magnetosphere (Beskin et al. 1992; Hirotani & Okamoto 1998; Levinson 2000; Neronov et al. 2005, 2009; Neronov & Aharonian 2007; Aleksić et al. 2014; Hirotani & Pu 2016; Broderick & Tchekhovskoy 2015; Ptitsyna & Neronov 2016). Contrary to the shock acceleration, these types of acceleration mechanisms produce high-energy particle spectra sharply peaked at a characteristic energy so that the effective $\Gamma$ is $\Gamma \ll 2$.

A straightforward difference between leptonic and hadronic models of blazar activity is the absence or presence of neutrino emission accompanying the $\gamma$-ray emission. No significant neutrino flux is expected in the leptonic models. To the contrary, the specific of hadronic cascades is that the overall neutrino energy flux is always comparable (at least within an order of magnitude) with the electromagnetic flux from the source.

IceCube neutrino telescope has discovered an astrophysical neutrino signal in the energy range in which AGN and, in particular, blazars are expected to produce neutrino flux (Aartsen et al. 2013, 2014b, 2015b,a,c; IceCube Collaboration 2013; IceCube Collaboration et al. 2016b). Search for neutrino signal from individual blazars or blazar samples did not provide positive detections up to now (Tchernin et al. 2013; IceCube Collaboration et al. 2016a). Also the stacking analysis did not provide solid evidence for neutrino flux from blazars (Aartsen et al. 2014a; Glüsenkamp 2016; Padovani et al. 2016). The absence of multiplet events in the IceCube muon neutrino data gives additional constraints on the density of neutrino sources, which also constrain blazars as possible neutrino sources (Murase & Waxman 2016). However, time coincidence of a neutrino arrival with a blazar outburst was claimed to support the hypothesis of the blazar origin of the astrophysical neutrino signal (Kadler et al. 2016).

Non-detection of the signal from the $\gamma$-ray brightest blazars was used by Tchernin et al. (2013) to derive limits on the parameters of the PIC hadronic models of blazars. The data of IC-40 detector limit $E_{cut}$ and $\Gamma$ to be roughly

$$E_{\text{cut}} \ge 10^{18} \text{ eV, } \Gamma \le 2. \qquad (1)$$

In what follows we combine the updated IceCube neutrino data (IceCube Collaboration et al. 2016b) with the data of Fermi $\gamma$-ray telescope to improve previously derived bounds on $\Gamma$ and $E_{cut}$. We use stacking, rather than individual source analysis, and consider cumulative $\gamma$-ray and neutrino spectra of a large number of northern sky blazars.

Our results severely limit the class of PIC models consistent with the combined $\gamma$-ray and neutrino data. Most of the hadronic component-dominated models of time-averaged flux from blazars, based on the assumption of shock acceleration of protons interacting either with low-energy protons or with the UV radiation field, appear to be ruled out. Models that avoid the $\gamma$-ray + neutrino constraint are those predicting the accelerated particle spectra sharply peaked in the UHECR range. However, in these models the parameters of particle acceleration mechanisms also have to be tuned to provide the highest possible efficiencies. The limits imposed by the IceCube data could be avoided in the “mixed” hadronic plus leptonic models of blazar

emission. For such models, our analysis severely limits the fraction of electromagnetic power that could come from hadronic interactions. The hadronic and leptonic components could be also “mixed in time”; for instance the hadronic emission might be responsible for flaring activity, while leptonic emission dominates the steady-state flux. Such models are also not directly constrained by our analysis. The contribution of such sources to IceCube flux is subdominant.

## 2. Stacking analysis approach to $\gamma$-ray and neutrino spectra

Within the PIC models, the overall energy flux of neutrinos and photons from the source is comparable, at least within an order of magnitude, i.e.

$$F_{\gamma} \sim F_{\gamma}. \qquad (2)$$

The neutrino signal from individual blazars is detectable in background-free regime in the energy band above several hundred TeV (IceCube Collaboration et al. 2016b). Suppose that the upper limit on neutrino flux from an individual source is $F_{v,lim}$. Non-detection of neutrinos from any of $N$ stacked sources provides an upper limit on “typical” neutrino flux from a source, i.e.

$$F_{\gamma} \le \frac{F_{v,lim}}{N}. \qquad (3)$$

The upper limits on the flux of individual sources in the Northern Hemisphere reported by IceCube (IceCube Collaboration et al. 2016a) range between $F_{v,lim} \approx 10^{-9} \text{ GeV/(cm}^2\text{s)}$ and $3 \times 10^{-9} \text{ GeV/(cm}^2\text{s)}$, for the declinations ranging from 0 to 90 degrees. The best sensitivity is achieved in the energy range around 100 TeV, where the atmospheric neutrino background becomes low enough so that the signal with the flux level $F_{v,lim}$ is detectable in nearly background free regime. The limiting flux level for individual sources is comparable to the typical flux level of blazars detected by Fermi-LAT telescope. Thus, source-by-source analysis could only mildly constrain the limits of the PIC type models by imposing a requirement that the bulk of the neutrino flux is not emitted in the 100 TeV band. This limits the properties of the parent proton spectrum, as reported by Tchernin et al. (2013).

A size of the population of GeV $\gamma$-ray detected blazars is $\sim 10^3$. Although the flux of most of the Fermi-LAT detected blazars is typically one or two orders of magnitude lower than the fluxes of the brightest blazars, the cumulative $\gamma$-ray flux of all the detected blazars is still an order of magnitude larger than that of an individual bright blazar. Thus, the expected cumulative neutrino signal from the population of the $\gamma$-ray detected blazars is a factor of 10 larger than the flux from any individual bright blazar. If the typical energies of neutrinos from blazars are much higher than 100 TeV, as suggested by the IceCube constraints derived from the analysis of individual bright blazars (Tchernin et al. 2013), the signal from individual blazars is not detectable, but an order of magnitude stronger cumulative neutrino signal from $\gamma$-ray loud blazar population might still be detectable at the highest energies.

## 3. Data analysis

### 3.1. Fermi-LAT

For our stacking analysis we consider publicly available data of Fermi-LAT telescope collected during the time period between
<figure><img src="image_2.png" /><figcaption>Fig. 1. IceCube muon neutrino events (white ellipses) with direction uncertainty less than 4 degrees, overlaid over Fermi-LAT countmap pf the northern sky in the energy range above 1 TeV, smoothed with 3 degree Gaussian. Green crosses show blazars selected for the stacking analysis.</figcaption></figure>

August 2008 and June 2016. We filter the LAT event list with the help of `gtselect` - `gtmktime` tools following the recommendations of the Fermi-LAT team<fnref n="1" /> to select only events belonging to the CLEAN sub-class. For each selected source we extract the spectrum using the aperture photometry method estimating the exposure with the help of `gtexposure` tool. We sum the source counts, the background counts and the exposures in the direction of all the selected sources to produce a cumulative source spectrum. For each source, the source signal is extracted from a circle of the radius 1° around the source position. The background is estimated from a 1° circle displaced by 2° from the source position.

The list of blazars selected for the stacking analysis includes $N = 749$ flat spectrum radio quasars (FSRQ) and BL Lacerta (BL Lac) type objects at declinations $Dec > -5°$ listed in the third Fermi-LAT source catalogue (Acero et al. 2015). This choice of the declination range is determined by sky region from which the IceCube astrophysical neutrino signal is collected (IceCube Collaboration et al. 2016b).

Figure 2 shows the resulting cumulative spectrum of blazars in this part of the sky. One could see that it is nearly identical to the cumulative spectrum of high Galactic latitude sources found in the Ackermann et al. (2015). This is not surprising, given the fact that blazars constitute the dominant extragalactic source population.

## 3.2. IceCube

We use the list of 29 muon neutrino events with energy proxies above 200 TeV reported by IceCube Collaboration et al. (2016b). Three out of the 29 events have large statistical uncertainty of direction reconstruction (larger than 3°). The chance coincidence probability for such events to have one of the 749 blazars within their error ellipse is of the order of one. We remove these events from the analysis. Rejection of these events reduces the effective exposure of IceCube data set by $3/29 \approx 10\%$.

The 90% error ellipses of the 26 neutrino events retained for the analysis do not contain γ-ray detected blazars, except for one blazar, OP 313, which is at the border of the error ellipse of the muon neutrino event 36. This is consistent with a chance coincidence expectation. The event number 36 has the

<figure><img src="image_3.png" /><figcaption>Fig. 2. Cumulative γ-ray (blue data points) spectrum and neutrino flux upper limit (red) for the Northern Hemisphere blazars. The blue shaded band shows the spectrum of extragalactic sources resolved by Fermi telescope from Ackermann et al. (2015). The black hatched bow-tie shows the IceCube astrophysical neutrino flux in the muon neutrino channel from IceCube Collaboration et al. (2016b). The grey dashed line shows a model neutrino spectrum for a PIC model from Tchernin et al. (2013).</figcaption></figure>

energy proxy $E = 200$ TeV and only 0.45 “signalness” value (IceCube Collaboration et al. 2016b), i.e. is more likely to be part of the atmospheric neutrino background. A 90% upper limit on the number of muon neutrino events with energy proxies above 200 TeV from blazars is, assuming the Poisson statistics of the signal,

$$N_{\text{lim}} = 4. \qquad (4)$$

The spectrum of neutrino emission from blazars is, in general, unknown. To derive an upper limit on the neutrino flux from blazars for arbitrary spectral shape, for different slopes $\Gamma$ we calculate the maximal possible normalisation $\kappa$ of a power-law neutrino flux

$$\frac{dN_{\nu}}{dE} = \kappa \left(\frac{E}{E_{*}}\right)^{-\Gamma}, \qquad (5)$$

where $E_*$ is normalisation energy fixed here to $E_* = 1$ PeV. We scan over the slopes $\Gamma$ to find an “envelope” curve of the

<footnote n="1">¹ http://fermi.gsfc.nasa.gov/ssc/data/analysis/</footnote>
different maximal possible flux power laws, as described by Tchernin et al. (2013). Realistic neutrino emission spectra can typically be well approximated by power laws in the energy range where the IceCube sensitivity is highest (about PeV for the considered IceCube data set), unless the spectrum has a high- or low-energy cut-off exactly in the IceCube sensitivity range. This implies that envelope curve of the upper limits on the power-law type spectra can be also used for the realistic spectra; the spectra consistent with the data can at most “touch” the envelope curve from below.

To derive the envelope curve, the IceCube exposure for $\nu_\mu$ (or $\bar{\nu}_\mu$) in the energy range above 400 TeV is well approximated by a power law

$$T_{\text{exp}}A_{\text{eff}} \approx TA_{*}(E_{\nu}/E_{*})^{p} \quad (6)$$

with the normalisation $TA_* \approx (7/2) \times 10^{14} \text{ cm}^2 \text{ s}$, averaged over the solid angle $\Omega = 2\pi(1 - \cos(95^\circ))$ and counting only muon neutrinos, at a reference energy $E_* = 1$ PeV and the slope $p = 0.34$ (IceCube Collaboration et al. 2016b).

The expected number of muon neutrinos and anti-neutrinos in an energy range $E_{\text{min}} < E_{\nu} < E_{\text{max}}$ for a given flux normalisation $\kappa$ is

$$ \begin{aligned} N_{\nu_{\mu}} &= \frac{1}{3}\kappa \int \Omega T_{\text{exp}}A_{\text{eff}}(E_{\nu}) \left(\frac{E_{\nu}}{E_{*}}\right)^{-\Gamma} dE_{\nu} \\ &= \frac{\kappa \Omega T A_{*} E_{*}}{3(p - \Gamma + 1)} \left(\left[\frac{E_{\text{max}}}{E_{*}}\right]^{p-\Gamma+1} - \left[\frac{E_{\text{min}}}{E_{*}}\right]^{p-\Gamma+1}\right), \end{aligned} \quad (7) $$

where the factor 1/3 accounts for the fact that the muon neutrinos constitute 1/3 of the signal (adopting standard assumptions about production mechanism and mixing).

Muons produced in the charged current interactions outside the IceCube detector have initial energies $E_\mu \approx (1-y_{cc})E_\nu$, where $y_{cc}$ is the average inelasticity of the charged current interactions (Gandhi et al. 1996). The energies of most detected neutrino-induced muons are much lower than the initial muon energy because of the energy loss in the rock $dE_\mu/dx = -(a+bE_\mu)$ (Chirkin & Rhode 2004), so that the distribution of final energies of muons originating from monoenergetic neutrinos with energy $E_\nu$ steadily undergoing charged current interactions all over a large distance through the rock/ice is $dN_\mu/dE \propto E^{-1}$ in the energy range $E < E_0$, down to the energy $E_{\text{crit}} = a/b \approx 1$ TeV. This suggests the probability density function for the muon energy (Neronov & Ribordy 2009b)

$$ \frac{d p(E_{\mu}, E_{\nu})}{d E} = \left( \ln \left( 1 + \frac{(1 - y_{cc}) E_{\nu}}{E_{\text{crit}}} \right) \right)^{-1} \frac{1}{(E_{\mu} + E_{\text{crit}})}. \quad (8) $$

The distribution of the energies of the muon events is then

$$ \begin{aligned} \frac{dN_{\mu}}{dE_{\mu}} &= \int_{\frac{E_{\mu}}{1+y_{cc}}}^{\infty} \frac{dp(E_{\mu}, E_{\nu})}{dE} T_{\text{exp}}A_{\text{eff}}(E_{\nu}) \Omega \frac{\kappa}{3} \left(\frac{E_{\nu}}{E_{*}}\right)^{-\Gamma} dE_{\nu} \\ &\approx \frac{\Omega T A_{*}}{3 \ln\left(\frac{E_{\mu}}{E_{\text{crit}}}\right) (p - \Gamma + 1)} \left(\frac{E_{\mu}}{E_{*}}\right)^{p-\Gamma}. \end{aligned} \quad (9) $$

Integrating the muon distribution in the energy range ($E_{\text{min}}, E_{\text{max}}$) one finds

$$ N_{\mu} \approx \frac{\kappa \Omega T A_{*} E_{*}}{3(\Gamma - p - 1)} \times \left( \left[ \frac{(E_{\text{min}}/E_{*})^{p-\Gamma+1}}{\ln(E_{\text{min}}/E_{\text{crit}})} \right] - \left[ \frac{(E_{\text{max}}/E_{*})^{p-\Gamma+1}}{\ln(E_{\text{max}}/E_{\text{crit}})} \right] \right), \quad (10) $$

which differs from the neutrino number by a factor $(\ln(E_\mu/E_{\text{crit}}))^{-1} \approx 0.2$.

Normalising $N_\mu = N_{\text{lim}}$ one finds $\kappa$ for different values of $\Gamma$. The result is shown by the red straight lines in Fig. 2. The red thick curve shows the envelope of all maximal allowed power laws. In $\Gamma < p + 1 \approx 1.3$, the signal statistics is dominated by the highest energy events and the maximal normalisation of the power law depends on the high-energy cut-off in the spectrum. We have fixed $E_{\text{max}}$ to $10^{18}$ eV in our calculation.

# 4. Results and discussion

## 4.1. Inconsistency of the data with proton induced cascade model with shock-accelerated protons in the jet

One could see from Fig. 2 that the upper limit on the neutrino flux from blazars in the PeV energy range is two orders of magnitude below the expected flux level. The mismatch between the $\gamma$-ray flux level and neutrino flux upper limits rule our hadronic models in which neutrino spectra are expected to produce most of the power in the PeV energy range.

The pion production in $\gamma\nu$ interactions occurs only above an energy threshold

$$ E_{p,\text{thr}} = 10^{16} \left[ \frac{\epsilon}{10 \text{ eV}} \right]^{-1}, \quad (11) $$

where $\epsilon$ is the characteristic energy of photons from the radiative environment of the AGN. The energies of neutrinos originating from the charged pion decays are typically below $\le 10\%$ of the parent proton energy,

$$ E_{\nu} \sim 10^{15} \left[ \frac{\epsilon}{10 \text{ eV}} \right]^{-1}. \quad (12) $$

Pion decays transfer nearly equal power to the neutrino and electromagnetic emission. Development of electromagnetic cascade in the source transfers the electromagnetic power to the GeV-TeV energy band. The horizontal dashed line in Fig. 2 shows an estimate of the electromagnetic power from the blazar population. Contrary to the electromagnetic power, the neutrino emission power stays in the energy band (12).

This is the case for the proton-induced cascade models in which the shock-accelerated protons interact with the radiation field of accretion disk. Conventional geometrically thin and optically thick accretion disks in AGN have temperatures reaching $10^4$ K in the innermost portions of the disk close to the last stable orbit (Shakura & Sunyaev 1976; Frank et al. 1992). Protons accelerated near the black hole or in the innermost portion of the AGN jet interact with the direct UV radiation from the disk or with the disk radiation scattered in the broad line region. This inevitably produces neutrinos with energies in the PeV range. An example of neutrino spectrum calculated for proton spectrum with $\Gamma = 2$ and $E_{\text{cut}} = 10^{17}$ eV interacting with the soft photons from the accretion disk with the spectrum peaking at $\epsilon = 15$ eV (from Tchernin et al. 2013) is shown by the grey dashed line in Fig. 2. A realistic neutrino spectrum produced by an entire population of blazars might be more complicated than this model because both $\Gamma$ and $E_{\text{cut}}$ are distributed in a certain range. Still a model with fixed $E_{\text{cut}}$ and $\Gamma$ could serve as a “reference” model for the entire population, where the reference parameter values are the distribution averages. This makes sense if the parameter distributions are concentrated around certain values, for example $\Gamma$ is always close to 2 in the models of shock acceleration.

The IceCube bound could be avoided if the bulk of neutrino power is emitted in an energy band different from 0.1–10 PeV.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of the IceCube upper limit on the neutrino flux with predictions of PIC models of blazars.</figcaption></figure>

This is the case if the soft photon targets for *py* interactions are in the infrared or microwave range. This type of interaction is considered in the model of Essey et al. (2010, 2011). The characteristic energy of the cosmic microwave background photons is $\epsilon \approx 10^{-3}$ eV. The magenta dotted line in Fig. 3 shows a representative model of neutrino spectrum from interactions of protons with the spectrum $E^{-2}$ with high-energy cut-off at $10^{20}$ eV during their propagation through the CMB and extragalactic background light radiation fields (Essey et al. 2011). Normalising the model flux on the average $\gamma$-ray flux from blazars, one finds that the model spectrum is consistent with the IceCube constraint.

Other types of hadronic models that are severely constrained by combined IceCube and Fermi data are those in which shock accelerated protons interact with the ambient medium protons and nuclei. In this case the pion production threshold is in the 100 MeV range and the neutrino spectrum in the energy range above the threshold approximately repeats the proton spectrum. If the proton spectrum is a power law with a slope close to 2, the neutrino spectrum in the PeV energy range is also a power law with a slope close to 2 and with the energy flux comparable to that of the $\gamma$-ray flux within a factor of two. The mismatch between the maximal possible normalisation of an $E^{-2}$ type neutrino spectrum and the dashed horizontal line in Figs. 2 and 3 is two orders of magnitude, which means that the model is ruled out. Consistency of the model with the data can only be restored if the acceleration process typically produces proton spectrum with a slope that is softer than 2.4. In this case the spectra of $\gamma$-rays and neutrinos originating from pion decays are also soft power laws and the neutrino flux in the PeV range are within the IceCube upper limit. This type of hadronic model is, however, constrained via source-by-source analysis, as reported for example by Tchernin et al. (2013).

Blazars are highly variable sources exhibiting periods of activity intermittent with the quiescence periods. Differences in the chains of production reactions of the PeV neutrinos and GeV photons could, in principle, lead to significant differences in the timing properties of the $\gamma$-ray and neutrino signals. In particular, higher energy neutrinos that are produced at the very beginning of development of the PIC might come before the lower energy photons emitted at the end of the cascade development. This mismatch between the timing properties of neutrino and $\gamma$-ray signals hampers the search of the neutrino emission from blazars through a search of the neutrino counterparts of the $\gamma$-ray

signal and the $\gamma$-ray counterparts of the neutrino signal on a source-by-source basis. However, it does not affect the results of the time-averaged stacking analysis presented above because the timescale of the analysis (eight years of Fermi data and six years of IceCube data) is much longer than the timescale of the $\gamma$-ray flares (minutes-to-months). All the neutrino signals associated with the $\gamma$-ray flares of individual blazars should be contained within the multi-year exposure with the possible exception of flares that occurred at the very beginning of the exposure time.

## 4.2. Constraints on models with sharply peaked proton spectra

The combined IceCube and Fermi constraint on the proton-induced cascade model can be avoided if the high-energy protons are not produced by the shock acceleration process. If the protons are injected by acceleration taking place in large-scale electric fields as is the case for the field in magnetic reconnection regions or in vacuum gaps in the black hole magnetosphere, the spectrum of protons is sharply peaked at a particular energy rather than having a power-law shape. If the characteristic energy of the acceleration process is large enough, the peak energy of the neutrino spectrum is determined by the characteristic proton energy, rather than by the threshold of the *py* reaction. The blue dash-dotted line in Fig. 3 shows the neutrino spectrum from interactions of protons with energies $3 \times 10^{19}$ GeV interacting with $\epsilon = 10$ eV photons. The low-energy part of the model spectrum is tangent to the envelope of the IceCube upper limits. This means that the models of PIC in the AGN accretion disk radiation fields are consistent with the limits if the proton spectrum has a sharp low-energy cut-off,

$$E_p \ge 2 \times 10^{19} \text{ eV (PIC in UV radiation field).} \quad (13)$$

An alternative possibility for avoiding the IceCube plus Fermi constraint is to consider models in which the accretion on the black hole forms a radiatively inefficient accretion flow (RIAF). These types of models are believed to be applicable to low luminosity radio galaxies and BL Lac type objects. In this case the low-energy radiation from the accretion flow is synchrotron radiation from electrons heated to relativistic temperatures by collisions with protons. Synchrotron radiation from a RIAF peaks in the infrared range according to Narayan et al. (1998) $\epsilon \le 0.1$ eV, and the neutrino spectrum peaks in the energy range above $10^{17}$ eV. This is illustrated by the red solid curve in Fig. 3, which is the neutrino spectrum produced in interactions of $10^{19}$ eV protons with 0.1 eV photons. The model spectrum is consistent with the IceCube data, so that models with proton spectrum with low-energy cut-off at

$$E_p \ge 0.6 \times 10^{19} \text{ eV (PIC in IR radiation field)} \quad (14)$$

are also consistent with the data.

The combined IceCube and Fermi constraint can be also avoided in the PIC model, in which high-energy protons interact with low-energy protons from the accretion flow if the high-energy proton spectrum is sharply peaked at high energies. The green dashed line in Fig. 3 shows the calculation of the model neutrino spectrum from interactions of protons with energy $2 \times 10^{20}$ eV with the low-energy protons. One could see that the low-energy part of the neutrino spectrum is tangent to the IceCube upper bound envelope curve. This means that models based on *pp* interactions are constrained to produce extremely high-energy protons,

$$E_p \ge 1 \times 10^{20} \text{ eV (PIC initiated by } pp \text{ interactions).} \quad (15)$$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Energies of protons accelerated in magnetospheric vacuum gaps near a black hole of mass $M = 5 \times 10^9 M_{\odot}$ (left column) or $5 \times 10^8 M_{\odot}$ (right column) surrounded by RIAF with synchrotron emission spectrum peaking in the infrared at $\epsilon = 0.1$ eV (top row) to $10^{-4}$ eV (bottom row). Red dashed lines show the dimensionless gap height $h$. Green solid curves show proton energies. Dark grey parts of the diagram show the parameter range excluded by the IceCube plus Fermi data set. Green shaded areas correspond to the allowed range of parameters.</figcaption></figure>

### 4.3. Proton energies in black hole magnetospheric gap models

Energies of protons accelerated in the gap of the height $H$ in the magnetosphere of a black hole of the mass $M$ are limited

by the finite extent of the gap, which is defined by the onset of electron-positron pair production on the soft photon background field present in the magnetosphere. The gap height depends on the luminosity $L$ and size $R$ of the soft photon field, the characteristic soft photon energy $\epsilon$, and the rate of
electron/proton acceleration in the gap, which is determined by the specific angular momentum of the black hole a (Beskin et al. 1992; Hirotani & Okamoto 1998; Levinson 2000; Neronov et al. 2005, 2009; Neronov & Aharonian 2007; Aleksić et al. 2014; Hirotani & Pu 2016; Broderick & Tchekhovskoy 2015; Ptitsyna & Neronov 2016).

Figure 4 shows the attainable proton energy as a function of source luminosity L and magnetic field B, calculated within the framework discussed by Ptitsyna & Neronov (2016). The assumption of the model is that the black hole accretes in the RIAF mode in which the soft photon field is produced via synchrotron emission from electrons heated to the temperatures 10–100 MeV by the protons. The RIAF synchrotron radiation spectrum typically peaks in the infrared range as opposed to the UV-dominated spectrum of an optically thick, geometrically thin accretion disk. The two columns of the figure correspond to two different black hole masses. The four rows correspond to various soft photon fields in the AGN central engine. In all the cases the infrared/microwave soft photon source is supposed to be distributed over a region of the size $10R_{Schw}$, where $R_{Schw}$ is the Schwarzschild radius of the black hole.

The process of the pair production starts to limit the gap height when the luminosity reaches certain (magnetic field dependent) value around $10^{40}$ erg/s. At lower luminosities of the RIAF the density of the soft photon field is not sufficient for the pair production within the extent of the black hole magnetosphere on the distance scale $R \sim R_{Schw}$ about the Schwarzschild radius of the black hole. In this case the limiting proton energy is estimated as

$$E_p \sim eBR_{Schw} \approx 10^{19} \left[ \frac{B}{100 G} \right] \left[ \frac{M}{3 \times 10^9 M_{\odot}} \right] \text{eV}, \quad (16)$$

where M is the black hole mass. The limit $E_p > 6 \times 10^{18}$ eV converts into a lower bound on magnetic field strength in the innermost part of the RIAF

$$B > 100 \left[ \frac{M}{3 \times 10^9 M_{\odot}} \right]^{-1} \text{G.} \quad (17)$$

Proton energies also drop below $6 \times 10^{18}$ eV at high values of B. In this case the proton synchrotron loss limits the energies of protons, as in the models of Mücke & Protheroe (2001), Aharonian (2002), Mücke et al. (2003). This type of model is not directly constrained by a combination of the Fermi-LAT and IceCube data discussed above.

The pair production initiated by electrons accelerated in the gap also gets suppressed at high luminosities because the strong inverse Compton loss rate does not allow electrons to accelerate to the energies needed for the pair production. Suppression of the pair production at low and high luminosities allows acceleration of protons to the energies in excess of $10^{19}$ eV. However, in this case protons start to produce pairs themselves when accelerated to energies higher than the pair production threshold. This also limits the proton energies in the case of high luminosity RIAF as is clear from Fig. 4.

The green shaded regions in the various panels of Fig. 4 show the ranges of L, B parameter space in which proton energies reach $>6 \times 10^{18}$ GeV. Acceleration in the magnetospheric vacuum gaps near black holes surrounded by RIAF with such parameters would produce neutrino and electromagnetic emission consistent with IceCube and Fermi-LAT data. The IceCube plus Fermi-LAT constraints could be satisfied in only a very limited range of parameter space.

## 5. Conclusions

We have shown that a combination of IceCube and Fermi-LAT data rules out certain types of hadronic models of emission from blazars. The models that are inconsistent with the data are those in which the observed γ-ray emission is produced by a particle cascade initiated by shock accelerated protons in the UV radiation field of the AGN central engine. This result is consistent with analysis of Aartsen et al. (2017), which appeared after this paper was completed. Our analysis confirms the main result of that paper, i.e. that the astrophysical neutrino flux is generated by sources different from the GeV γ-ray loud blazars (BL Lacs and flat spectrum radio quasars).

Hadronic models of steady-state blazar emission consistent with the data are those in which high-energy protons spectra are sharply peaked in the UHECR range. An example of this type of model is that of proton acceleration in the vacuum gaps of black hole magnetospheres. We have shown that the IceCube and Fermi-LAT data constrain the parameter space of such models (luminosity and magnetic field in the RIAF surrounding the black hole). Models that are consistent with the data predict a neutrino flux originating from UHECR production in the blazars. This suggests that the model is testable via observations of neutrino flux in the energy range higher than that accessible with IceCube, $E_{\gamma} \sim 0.1-1$ eV. An increase of IceCube exposure or an exploration of this energy range with dedicated detectors optimised for the 0.1–1 eV range, such as CHANT (Neronov et al. 2017) or ARA (Ara Collaboration et al. 2012) could be used to test the model.

Our analysis also does not rule out the possibility of hadronic flares of blazars. It is still possible that the steady-state emission that dominates the time-averaged flux is determined by leptonic processes, while episodes of proton acceleration followed by PIC are responsible for the flaring activity. Such a possibility could be constrained via the search of neutrino counterparts of the γ-ray flares (or γ-ray counterparts of neutrino emission episodes of individual sources, but not with the time-averaged/stacking analysis. In any case, the hadronic flare contribution to the IceCube time-averaged flux should be small.

*Acknowledgements.* We would like to thank S. Shoening for useful discussions on details of IceCube results. The work of K.P. is supported by the Russian Science Foundation grant 14-12-01340.

## References

Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013, Phys. Rev. Lett., 111, 021103

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014a, ApJ, 796, 109

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014b, Phys. Rev. Lett., 113, 101101

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2015a, ApJ, 809, 98

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2015b, Phys. Rev. D, 91, 022001

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2015c, Phys. Rev. Lett., 115, 081102

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017, ApJ, 835, 45

Aero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23

Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 799, 86

Aharonian, F. A. 2002, MNRAS, 332, 215

Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014, Science, 346, 1080

Ara Collaboration, Allison, P., Auffenberg, J., et al. 2012, Astropart. Phys., 35, 457

Bednarz, J., & Ostrowski, M. 1998, Phys. Rev. Lett., 80, 3911

Begelman, M. C., Rudak, B., & Sikora, M. 1990, ApJ, 362, 38

Beskin, V. S., Istomin, Y. N., & Parev, V. I. 1992, Sov. Astron., 36, 642

Broderick, A. E., & Tchekhovskoy, A. 2015, ApJ, 809, 97

Chirkin, D., & Rhode, W. 2004, High Energy Physics – Phenomenology [arXiv:hep-ph/0407075]
Eichler, D. 1979, ApJ, 232, 106

Essey, W., Kalashev, O. E., Kusenko, A., & Beacom, J. F. 2010, Phys. Rev. Lett., 104, 141102

Essey, W., Kalashev, O., Kusenko, A., & Beacom, J. F. 2011, ApJ, 731, 51

Frank, J., King, A., & Raine, D. 1992, Accretion power in astrophysics

Gandhi, R., Quigg, C., Hall Reno, M., & Sarcevic, I. 1996, Astropart. Phys., 5, 81

Glüsenkamp, T. 2016, in EPJ Web Conf., 121, 05006

Halzen, F., & Zas, E. 1997, ApJ, 488, 669

Hirotani, K., & Okamoto, I. 1998, ApJ, 497, 563

Hirotani, K., & Pu, H.-Y. 2016, ApJ, 818, 50

IceCube Collaboration 2013, Science, 342, 1242856

IceCube Collaboration, Aartsen, M. G., Abraham, K., et al. 2016a, ApJ, 835, 151

IceCube Collaboration, Aartsen, M. G., Abraham, K., et al. 2016b, ApJ, 833, 3

Kadler, M., Krauß, F., Mannheim, K., et al. 2016, Nat. Phys., 12, 807

Kalashev, O., Semikoz, D., & Tkachev, I. 2015, J. Exp. Theor. Phys., 120, 541

Lemoine, M., Pelletier, G., & Revenu, B. 2006, ApJ, 645, L129

Lesch, H., & Birk, G. T. 1997, A&A, 324, 461

Levinson, A. 2000, Phys. Rev. Lett., 85, 912

Mannheim, K. 1993, A&A, 269, 67

Mannheim, K., & Biermann, P. L. 1989, A&A, 221, 211

Mannheim, K., & Biermann, P. L. 1992, A&A, 253, L21

Mücke, A., & Protheroe, R. J. 2001, Astropart. Phys., 15, 121

Mücke, A., Protheroe, R. J., Engel, R., Rachen, J. P., & Stanev, T. 2003, Astropart. Phys., 18, 593

Murase, K., & Waxman, E. 2016, Phys. Rev. D, 94, 103006

Narayan, R., Mahadevan, R., & Quataert, E. 1998, in Theory of Black Hole Accretion Disks, eds. M. A. Abramowicz, G. Björnsson, & J. E. Pringle, 148

Neronov, A., & Aharonian, F. A. 2007, ApJ, 671, 85

Neronov, A., & Ribordy, M. 2009a, Phys. Rev. D, 80, 083008

Neronov, A., & Ribordy, M. 2009b, Phys. Rev. D, 79, 043013

Neronov, A. Y., & Semikoz, D. V. 2002, Phys. Rev. D, 66, 123003

Neronov, A., Tinyakov, P., & Tkachev, I. 2005, Sov. J. Exp. Theor. Phys., 100, 656

Neronov, A., Semikoz, D., & Sibiryakov, S. 2008, MNRAS, 391, 949

Neronov, A. Y., Semikoz, D. V., & Tkachev, I. I. 2009, New J. Phys., 11, 065015

Neronov, A., Semikoz, D. V., Anchordoqui, L. A., Adams, J. H., & Olinto, A. V. 2017, Phys. Rev. D, 95, 023004

Padovani, P., Resconi, E., Giommi, P., Arsioli, B., & Chang, Y. L. 2016, MNRAS, 457, 3582

Pelletier, G., Lemoine, M., & Marcowith, A. 2009, MNRAS, 393, 587

Ptitsyna, K., & Neronov, A. 2016, A&A, 593, A8

Romanova, M. M., & Lovelace, R. V. E. 1992, A&A, 262, 26

Shakura, N. I., & Sunyaev, R. A. 1976, MNRAS, 175, 613

Tchernin, C., Aguilar, J. A., Neronov, A., & Montaruli, T. 2013, A&A, 555, A70

Urry, C. M., & Padovani, P. 1995, PASP, 107, 803