Draft Version November 15, 2016: Preprint Typeset Using L TEX Style Emulateapj v. 01/23/15

Draft version November 15, 2016
Preprint typeset using LATEX style emulateapj v. 01/23/15
THE CONTRIBUTION OF FERMI-2LAC BLAZARS TO THE DIFFUSE TEV-PEV NEUTRINO FLUX

IceCube Collaboration: M. G. Aartsen1 , K. Abraham2 , M. Ackermann3 , J. Adams4 , J. A. Aguilar5 , M. Ahlers6 ,
M. Ahrens7 , D. Altmann8 , K. Andeen9 , T. Anderson10 , I. Ansseau5 , G. Anton8 , M. Archinger11 , C. Arguelles12 ,
T. C. Arlen10 , J. Auffenberg13 , S. Axani12 , X. Bai14 , S. W. Barwick15 , V. Baum11 , R. Bay16 , J. J. Beatty17,18 ,
J. Becker Tjus19 , K.-H. Becker20 , S. BenZvi21 , P. Berghaus22 , D. Berley23 , E. Bernardini3 , A. Bernhard2 ,
D. Z. Besson24 , G. Binder25,16 , D. Bindig20 , M. Bissok13 , E. Blaufuss23 , S. Blot3 , D. J. Boersma26 , C. Bohm7 ,
M. Börner27 , F. Bos19 , D. Bose28 , S. Böser11 , O. Botner26 , J. Braun6 , L. Brayeur29 , H.-P. Bretz3 , A. Burgman26 ,
J. Casey30 , M. Casier29 , E. Cheung23 , D. Chirkin6 , A. Christov31 , K. Clark32 , L. Classen33 , S. Coenders2 ,
G. H. Collin12 , J. M. Conrad12 , D. F. Cowen10,34 , A. H. Cruz Silva3 , J. Daughhetee30 , J. C. Davis17 , M. Day6 ,
arXiv:1611.03874v1 [astro-ph.HE] 11 Nov 2016
J. P. A. M. de André35 , C. De Clercq29 , E. del Pino Rosendo11 , H. Dembinski36 , S. De Ridder37 , P. Desiati6 ,

K. D. de Vries29 , G. de Wasseige29 , M. de With38 , T. DeYoung35 , J. C. Dı́az-Vélez6 , V. di Lorenzo11 ,
H. Dujmovic28 , J. P. Dumm7 , M. Dunkman10 , B. Eberhardt11 , T. Ehrhardt11 , B. Eichmann19 , S. Euler26 ,
P. A. Evenson36 , S. Fahey6 , A. R. Fazely39 , J. Feintzeig6 , J. Felde23 , K. Filimonov16 , C. Finley7 , S. Flis7 ,
C.-C. Fösig11 , A. Franckowiak3 , T. Fuchs27 , T. K. Gaisser36 , R. Gaior40 , J. Gallagher41 , L. Gerhardt25,16 ,
K. Ghorbani6 , W. Giang42 , L. Gladstone6 , M. Glagla13 , T. Glüsenkamp3,* , A. Goldschmidt25 , G. Golup29 ,
J. G. Gonzalez36 , D. Góra3 , D. Grant42 , Z. Griffith6 , C. Haack13 , A. Haj Ismail37 , A. Hallgren26 , F. Halzen6 ,
E. Hansen43 , B. Hansmann13 , T. Hansmann13 , K. Hanson6 , D. Hebecker38 , D. Heereman5 , K. Helbing20 ,
R. Hellauer23 , S. Hickford20 , J. Hignight35 , G. C. Hill1 , K. D. Hoffman23 , R. Hoffmann20 , K. Holzapfel2 ,
A. Homeier44 , K. Hoshina6,54 , F. Huang10 , M. Huber2 , W. Huelsnitz23 , K. Hultqvist7 , S. In28 , A. Ishihara40 ,
E. Jacobi3 , G. S. Japaridze45 , M. Jeong28 , K. Jero6 , B. J. P. Jones12 , M. Jurkovic2 , A. Kappes33 , T. Karg3 ,
A. Karle6 , U. Katz8 , M. Kauer6,46 , A. Keivani10 , J. L. Kelley6 , J. Kemp13 , A. Kheirandish6 , M. Kim28 ,
T. Kintscher3 , J. Kiryluk47 , T. Kittler8 , S. R. Klein25,16 , G. Kohnen48 , R. Koirala36 , H. Kolanoski38 ,
R. Konietz13 , L. Köpke11 , C. Kopper42 , S. Kopper20 , D. J. Koskinen43 , M. Kowalski38,3 , K. Krings2 , M. Kroll19 ,
G. Krückl11 , C. Krüger6 , J. Kunnen29 , S. Kunwar3 , N. Kurahashi49 , T. Kuwabara40 , M. Labare37 ,
J. L. Lanfranchi10 , M. J. Larson43 , D. Lennarz35 , M. Lesiak-Bzdak47 , M. Leuermann13 , J. Leuner13 , L. Lu40 ,
J. Lünemann29 , J. Madsen50 , G. Maggi29 , K. B. M. Mahn35 , S. Mancina6 , M. Mandelartz19 , R. Maruyama46 ,
K. Mase40 , R. Maunu23 , F. McNally6 , K. Meagher5 , M. Medici43 , M. Meier27 , A. Meli37 , T. Menne27 , G. Merino6 ,
T. Meures5 , S. Miarecki25,16 , E. Middell3 , L. Mohrmann3 , T. Montaruli31 , M. Moulai12 , R. Nahnhauer3 ,
U. Naumann20 , G. Neer35 , H. Niederhausen47 , S. C. Nowicki42 , D. R. Nygren25 , A. Obertacke Pollmann20 ,
A. Olivas23 , A. Omairat20 , A. O’Murchadha5 , T. Palczewski51 , H. Pandya36 , D. V. Pankova10 , Ö. Penek13 ,
J. A. Pepper51 , C. Pérez de los Heros26 , C. Pfendner17 , D. Pieloth27 , E. Pinat5 , J. Posselt20 , P. B. Price16 ,
G. T. Przybylski25 , M. Quinnan10 , C. Raab5 , L. Rädel13 , M. Rameez31 , K. Rawlins52 , R. Reimann13 , M. Relich40 ,
E. Resconi2 , W. Rhode27 , M. Richman49 , B. Riedel42 , S. Robertson1 , M. Rongen13 , C. Rott28 , T. Ruhe27 ,
D. Ryckbosch37 , D. Rysewyk35 , L. Sabbatini6 , S. E. Sanchez Herrera42 , A. Sandrock27 , J. Sandroos11 ,
S. Sarkar43,53 , K. Satalecka3 , M. Schimp13 , P. Schlunder27 , T. Schmidt23 , S. Schoenen13 , S. Schöneberg19 ,
A. Schönwald3 , L. Schumacher13 , D. Seckel36 , S. Seunarine50 , D. Soldin20 , M. Song23 , G. M. Spiczak50 ,
C. Spiering3 , M. Stahlberg13 , M. Stamatikos17,55 , T. Stanev36 , A. Stasik3 , A. Steuer11 , T. Stezelberger25 ,
R. G. Stokstad25 , A. Stößl3 , R. Ström26 , N. L. Strotjohann3 , G. W. Sullivan23 , M. Sutherland17 , H. Taavola26 ,
I. Taboada30 , J. Tatar25,16 , S. Ter-Antonyan39 , A. Terliuk3 , G. Tešić10 , S. Tilav36 , P. A. Toale51 , M. N. Tobin6 ,
S. Toscano29 , D. Tosi6 , M. Tselengidou8 , A. Turcati2 , E. Unger26 , M. Usner3 , S. Vallecorsa31 ,
J. Vandenbroucke6 , N. van Eijndhoven29 , S. Vanheule37 , M. van Rossem6 , J. van Santen3 , J. Veenkamp2 ,
M. Vehring13 , M. Voge44 , M. Vraeghe37 , C. Walck7 , A. Wallace1 , M. Wallraff13 , N. Wandkowsky6 ,
Ch. Weaver42 , C. Wendt6 , S. Westerhoff6 , B. J. Whelan1 , S. Wickmann13 , K. Wiebe11 , C. H. Wiebusch13 ,
L. Wille6 , D. R. Williams51 , L. Wills49 , H. Wissing23 , M. Wolf7 , T. R. Wood42 , E. Woolsey42 , K. Woschnagg16 ,
D. L. Xu6 , X. W. Xu39 , Y. Xu47 , J. P. Yanez3 , G. Yodh15 , S. Yoshida40 , and M. Zoll7
Draft version November 15, 2016
ABSTRACT
The recent discovery of a diffuse cosmic neutrino flux extending up to PeV energies raises the ques-
tion of which astrophysical sources generate this signal. One class of extragalactic sources which may
produce such high-energy neutrinos are blazars. We present a likelihood analysis searching for cumu-
lative neutrino emission from blazars in the 2nd Fermi-LAT AGN catalogue (2LAC) using an IceCube
neutrino dataset 2009-12 which was optimised for the detection of individual sources. In contrast to
previous searches with IceCube, the populations investigated contain up to hundreds of sources, the
largest one being the entire blazar sample in the 2LAC catalogue. No significant excess is observed
and upper limits for the cumulative flux from these populations are obtained. These constrain the
maximum contribution of the 2LAC blazars to the observed astrophysical neutrino flux to be 27%
or less between around 10 TeV and 2 PeV, assuming equipartition of flavours at Earth and a single
power-law spectrum with a spectral index of −2.5. We can still exclude that the 2LAC blazars (and
sub-populations) emit more than 50% of the observed neutrinos up to a spectral index as hard as −2.2
in the same energy range. Our result takes into account that the neutrino source count distribution
is unknown, and it does not assume strict proportionality of the neutrino flux to the measured 2LAC
γ-ray signal for each source. Additionally, we constrain recent models for neutrino emission by blazars.
2 M. G. Aartsen et al.
1. INTRODUCTION with increasing accuracy (Aartsen et al. 2015a). The

The initial discovery of the astrophysical neutrino flux most recent results indicate a soft spectrum with a spec-
around PeV energies (Aartsen et al. 2013a) a few years tral index of −2.5±0.1 between around 10 TeV and 2 PeV
ago marked the beginning of high-energy neutrino as- with no significant deviation from an equal flavor com-
tronomy. Since then, its properties have been measured position at Earth (Aartsen et al. 2015b). The neutrino
signal has been found to be compatible with an isotropic
1 Department of Physics, University of Adelaide, Adelaide, distribution on the sky. This apparent isotropy suggests
5005, Australia that a significant fraction of the observed neutrinos is of
2 Physik-department, Technische Universität München, D-
extragalactic origin, a result which is also supported by
85748 Garching, Germany Ahlers et al. (2016). However, there are also indications
3 DESY, D-15735 Zeuthen, Germany
4 Dept. of Physics and Astronomy, University of Canterbury, for a 3-σ anisotropy (Neronov & Semikoz 2016) if low-
Private Bag 4800, Christchurch, New Zealand energy events (< 100 TeV) are omitted. Further data
5 Université Libre de Bruxelles, Science Faculty CP230, B-
are required to settle this issue.
1050 Brussels, Belgium Potential extragalactic source candidates are Active
6 Dept. of Physics and Wisconsin IceCube Particle Astro-
physics Center, University of Wisconsin, Madison, WI 53706, Galactic nuclei (AGN), where both radio-quiet (Stecker
USA et al. 1991) and radio-loud (Mannheim 1995) objects
7 Oskar Klein Centre and Dept. of Physics, Stockholm Univer-
have been considered for neutrino production since many
sity, SE-10691 Stockholm, Sweden
8 Erlangen Centre for Astroparticle Physics, Friedrich- years. Blazars, a subset of radio-loud active galactic nu-
Alexander-Universität Erlangen-Nürnberg, D-91058 Erlangen, clei with relativistic jets pointing towards Earth (Urry &
Germany Padovani 1995), are investigated in this paper. They are
9 Department of Physics, Marquette University, Milwaukee,
commonly classified based on the properties of the spec-
WI, 53201, USA
10 Dept. of Physics, Pennsylvania State University, University tral energy distribution (SED) of their electromagnetic
Park, PA 16802, USA emission. The blazar SED features two distinctive peaks:
11 Institute of Physics, University of Mainz, Staudinger Weg
a low-energy peak between infrared and X-ray energies,
7, D-55099 Mainz, Germany attributed to synchrotron emission of energetic electrons,
12 Dept. of Physics, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA and a high-energy peak at γ-ray energies, which can
13 III. Physikalisches Institut, RWTH Aachen University, D- be explained by several and possibly competing inter-
52056 Aachen, Germany
14 Physics Department, South Dakota School of Mines and
Technology, Rapid City, SD 57701, USA University, University Park, PA 16802, USA
15 Dept. of Physics and Astronomy, University of California, 35 Dept. of Physics and Astronomy, Michigan State Univer-
Irvine, CA 92697, USA sity, East Lansing, MI 48824, USA
16 Dept. of Physics, University of California, Berkeley, CA 36 Bartol Research Institute and Dept. of Physics and Astron-
94720, USA omy, University of Delaware, Newark, DE 19716, USA
17 Dept. of Physics and Center for Cosmology and Astro- 37 Dept. of Physics and Astronomy, University of Gent, B-
Particle Physics, Ohio State University, Columbus, OH 43210, 9000 Gent, Belgium
USA 38 Institut für Physik, Humboldt-Universität zu Berlin, D-
18 Dept. of Astronomy, Ohio State University, Columbus, OH
12489 Berlin, Germany
43210, USA 39 Dept. of Physics, Southern University, Baton Rouge, LA
19 Fakultät für Physik & Astronomie, Ruhr-Universität
70813, USA
Bochum, D-44780 Bochum, Germany 40 Dept. of Physics, Chiba University, Chiba 263-8522, Japan
20 Dept. of Physics, University of Wuppertal, D-42119 Wup- 41 Dept. of Astronomy, University of Wisconsin, Madison, WI
pertal, Germany 53706, USA
21 Dept. of Physics and Astronomy, University of Rochester, 42 Dept. of Physics, University of Alberta, Edmonton, Alberta,
Rochester, NY 14627, USA Canada T6G 2E1
22 National Research Nuclear University MEPhI (Moscow En- 43 Niels Bohr Institute, University of Copenhagen, DK-2100
gineering Physics Institute), Moscow, Russia Copenhagen, Denmark
23 Dept. of Physics, University of Maryland, College Park, MD 44 Physikalisches Institut, Universität Bonn, Nussallee 12, D-
20742, USA 53115 Bonn, Germany

24 Dept. of Physics and Astronomy, University of Kansas, 45 CTSPS, Clark-Atlanta University, Atlanta, GA 30314, USA
46 Dept. of Physics, Yale University, New Haven, CT 06520,
Lawrence, KS 66045, USA
25 Lawrence Berkeley National Laboratory, Berkeley, CA USA
47 Dept. of Physics and Astronomy, Stony Brook University,
94720, USA
26 Dept. of Physics and Astronomy, Uppsala University, Box Stony Brook, NY 11794-3800, USA
48 Université de Mons, 7000 Mons, Belgium
516, S-75120 Uppsala, Sweden
27 Dept. of Physics, TU Dortmund University, D-44221 Dort- 49 Dept. of Physics, Drexel University, 3141 Chestnut Street,
mund, Germany Philadelphia, PA 19104, USA

28 Dept. of Physics, Sungkyunkwan University, Suwon 440- 50 Dept. of Physics, University of Wisconsin, River Falls, WI
746, Korea 54022, USA

29 Vrije Universiteit Brussel, Dienst ELEM, B-1050 Brussels, 51 Dept. of Physics and Astronomy, University of Alabama,
Belgium Tuscaloosa, AL 35487, USA

30 School of Physics and Center for Relativistic Astrophysics, 52 Dept. of Physics and Astronomy, University of Alaska An-
Georgia Institute of Technology, Atlanta, GA 30332, USA chorage, 3211 Providence Dr., Anchorage, AK 99508, USA
31 Département de physique nucléaire et corpusculaire, Uni- 53 Dept. of Physics, University of Oxford, 1 Keble Road, Ox-
versité de Genève, CH-1211 Genève, Switzerland ford OX1 3NP, UK

32 Dept. of Physics, University of Toronto, Toronto, Ontario, 54 Earthquake Research Institute, University of Tokyo,
Canada, M5S 1A7 Bunkyo, Tokyo 113-0032, Japan

33 Institut für Kernphysik, Westfälische Wilhelms-Universität 55 NASA Goddard Space Flight Center, Greenbelt, MD 20771,
Münster, D-48149 Münster, Germany USA

34 Dept. of Astronomy and Astrophysics, Pennsylvania State * corresponding author: thorsten.gluesenkamp@fau.de
2LAC-blazar contribution to TeV-PeV neutrinos 3
action and radiation processes of high-energy electrons ermann et al. 2011) to define search positions for our
and high-energy nuclei (Böttcher et al. 2013). Several analysis (see section 2). The blazars in the 2LAC cata-
works suggest that blazar SEDs follow a sequence (Fos- logue comprise the majority (≈ 70%) of the total γ-ray
sati et al. 1998; Böttcher & Dermer 2002; Cavaliere & flux emitted from all GeV-blazars in the observable uni-
D’Elia 2002; Meyer et al. 2011), in which the peak energy verse between 100 MeV and 100 GeV (see appendix C).
of the synchrotron emission spectrum decreases with in- Compared to other Fermi catalogues starting at higher
creasing blazar luminosity. Accordingly, blazars can be energies, 1FHL (Ackermann et al. 2013) or 2FHL (Ack-
classified into low synchrotron peak (LSP), intermedi- ermann et al. 2016), the 2LAC contains more than twice
ate synchrotron peak (ISP) and high synchrotron peak the number of blazars. The goal is to look for a cumu-
(HSP) objects57 , a classification scheme introduced in lative neutrino flux excess from all 862 2LAC blazars or
Abdo et al. (2010a) which we use throughout this work. from specifically selected sub-populations using muon-
A second classifier is based on the prominence of emission track data with an angular resolution of about a degree
lines in the SED over the non-thermal continuum emis- in an unbinned maximum-likelihood stacking approach.
sion of the jet. Flat Spectrum Radio Quasars (FSRQs) We use two different ”weighting schemes” (see section
show Doppler-broadened optical emission lines (Stickel 4.2) to define the probability density functions (PDFs)
et al. 1991), while in so called BL Lac objects emission for the neutrino signal, expressing different assumptions
lines are hidden under a strong continuum emission. about the relative neutrino flux for each source. Each
Many calculations of high-energy neutrino emission weighting scheme represents its own test of the data.
from the jets of blazars can be found in the litera- The analysis differs from previous point source searches
ture. Neutrinos could be produced via charged pion most drastically in two points.
decay in interactions of high-energy protons with gas
(pp-interactions) in the jets (Schuster et al. 2002) or in 1. The blazar populations comprise nearly 2 orders of
interactions of protons with internal (Mannheim 1995) magnitude more sources.
or external (Atoyan & Dermer 2001) photon fields (pγ- 2. For the first time, we use a model-independent
interactions). Early models for the neutrino emission weighting scheme. In this test of the data, we make
from blazars made no explicit distinction based on the nearly no assumption about the exact ν/γ correla-
blazar class. Some of these have already been explicitly tion, except that the neutrino flux originates from
excluded at 90% C.L. by past diffuse neutrino flux mea- the defined blazar positions.
surements (Abbasi et al. 2011; Aartsen et al. 2014a), for
example the combined pp+pγ predictions in Mannheim Section 2 defines the five blazar populations consid-
(1995). More recent publications, on the other hand, ered in this analysis. Section 3 describes the muon track
differentiate between specific classes of blazars and are dataset used for this search. Section 4 summarizes the
largely not constrained by experiment, yet. The neu- analysis, including the technique of the unbinned stack-
trino production of BL Lac objects is modeled e.g. in ing search, a description of different weighting schemes,
Mücke et al. (2003); Tavecchio et al. (2014); Padovani the confidence interval construction, and a discussion on
et al. (2015) while neutrino production of FSRQs is cal- potential biases from non-hadronic contributions to the
culated e.g. in Becker et al. (2005); Murase et al. (2014). γ-ray flux. Section 5 presents the analysis results and
The models by Tavecchio et al. (2014) and Padovani et al. section 6 discusses their implications.
(2015) were in particular constructed to explain parts or
all of the astrophysical neutrino flux. With the analy- 2. THE BLAZAR POPULATIONS
sis presented here, we are able to test large parts of the The Fermi-LAT 2LAC catalogue (Ackermann et al.
parameter space of many of these models for the first 2011) contains 862 GeV-emitting blazars at high galac-
time. We do not consider theoretical calculations from tic latitudes |b| > 10◦ that are not affected by potential
the literature for individual objects, since these are not source confusion59 . The data for this catalog was taken
directly comparable to our results. between August 2008 and August 2010. We use the spec-
The neutrinos predicted by most models are produced troscopic classification into FSRQ and BL Lac objects
in charged pion decays which come with an associated (Stickel et al. 1991) and the independent classification
flux of γ-rays from neutral pion decays. Even if the into LSP, ISP and HSP objects (Abdo et al. 2010a) to
hadronic fraction is sub-dominant, one could on average define sub-populations of these 862 objects. We do not
expect a higher neutrino luminosity for a higher observed impose any other cuts (e.g. on the γ-ray flux) because
γ-luminosity (Murase et al. 2014). On a source-by-source the exact neutrino flux expectations are unknown as out-
basis, however, variations in the exact ν/γ correlation lined in section 1. The motivations for the particular
are likely. One strategy to cope with this uncertainty, sub-samples are described in the following.
which we follow in this paper, is to analyze large sam-
ples of objects and thereby to investigate average prop- All 2LAC Blazars (862 objects): The evolutionary
erties. We use the Fermi-LAT 2LAC catalogue58 (Ack- blazar sequence (Cavaliere & D’Elia 2002; Böttcher
& Dermer 2002) suggests that blazars form a con-
57 This scheme is a generalization of the XBL/RBL classification
tinuous spectrum of objects that are connected via
of BL Lac objects introduced by Padovani & Giommi (1995). cosmological evolution. A recent study by Ajello
58 The successor catalogue 3LAC (Ackermann et al. 2015) was
not yet published when this analysis was carried out. For the γ-
et al. (2014) supports this hypothesis. Since the
weighting scheme (see section 4.2.1), the results are expected to be corresponding evolution of the neutrino emission
nearly identical. The 2LAC sample already resolves the majority of
the GeV-blazar flux and the brightest blazars are also expected to 59 No source confusion means that the CLEAN flag from the cat-
be bright in the 3LAC catalogue in the quasi-steady approximation. alogue for a particular source is set.
is not known, the most unbiased assumption is to dataset all sky northern sky southern sky
group all blazars together. This is especially jus- IC-59 107011 42781 64230
tified for the analysis using the equal weighting IC-79 93720 48782 44938
scheme discussed in section 4.2. IC-86 136245 61325 74920
FSRQs (310 objects): The class of FSRQs show Table 1

Total number of data events in the respective datasets of IC-59,
strong, broad emission lines that potentially act as IC-79 and IC-86 for each celestial hemisphere. “Northern sky”
intense radiation targets for photomeson produc- means the zenith angle θ for the incoming particle directions is
tion of neutrinos (Atoyan & Dermer 2001; Murase equal to or larger than 90◦ . “Southern sky” means θ < 90◦ .
et al. 2014).
with the data aqcuisition system at the surface. The
LSPs (308 objects): The majority of FSRQs are LSP photo sensors detect Cherenkov light emitted by charged
objects. Giommi et al. (2012) argue that LSP- particles that are produced in the interactions of neutri-
BL Lacs are actually physically similar to FSRQs, nos with nuclei and electrons in the glacial ice. The ge-
but whose emission lines are overwhelmed by the ometry and sensitivity of the photo sensors leads to an ef-
strong jet continuum. This sample therefore groups fective energy threshold for neutrinos of about 100 GeV.
all LSP objects together. A more detailed description of the detector and the data
acquisition can be found in Abbasi et al. (2009).
ISP+HSPs (301 objects): HSP objects differ from Two main signatures can be distinguished for the
LSP objects in terms of their luminosities and recorded events: “track-like” and “shower-like”. Only
mainly consist of BL Lacs (Ajello et al. 2014). The track-like events are of interest for the analysis here.
peak-frequency boundary between LSP and HSP is They are the characteristic signature of muons produced
only defined artificially, with ISP objects filling the in the charged-current interactions of muon neutrinos60 .
gap. In order to have a larger sample of objects, IceCube was constructed between 2006 and 2011 with
the HSP objects are grouped together with the ISP a final configuration of 86 strings. We use data from
objects in one combined sample for this analysis. the 59-string (IC-59), 79-string (IC-79) and 86-string
LSP-BL Lacs (68 objects): Objects that share the (IC-86) configurations of the IceCube detector recorded
LSP and BL Lac classification have been specif- between May 2009 and April 2012. In contrast to pre-
ically considered for neutrino emission in Mücke vious publications, we do not include data from the 40-
et al. (2003). Therefore we test them as a separate string configuration here since the ice model description
sample. They form the smallest sub population in in the IC-40 MonteCarlo datasets is substantially dif-
this analysis. ferent and the sensitivity gain would be marginal. The
track event selection for the three years of data is sim-
The distribution of the sources on the sky for the ilar to the ones described in Aartsen et al. (2013b) and
largest sample (all 2LAC blazars) and smallest sample Aartsen et al. (2014b). The angular resolution of the ma-
(LSP-BL Lacs) are shown in figure 1. A modest LAT- jority of events in the track sample is better than 1◦ for
exposure deficit and lower sky-coverage by optical sur- events with reconstructed muon energies above 10 TeV
veys in the southern sky lead to a slight deficit of objects (Aartsen et al. 2014b). The angular reconstruction un-
in southern hemisphere (Ackermann et al. 2011). The certainty is calculated following the prescription given in
effect is most prominent for the BL Lac-dominated sam- Neunhoffer (2006). We apply one additional minor se-
ples. However, blazars without optical association are lection criterion for the estimated angular uncertainty of
also included in the 2LAC catalog and partly make up for the reconstructed tracks (σest. ≤ 5◦ ) for computational
this asymmetry in the total sample. For simplicity, we as- reasons. The removed events do not have any measurable
sume a quasi-isotropic source distribution for all popula- effect on the sensitivity. Event numbers for the individ-
tions (excluding the source-free region around the galac- ual datasets are summarized in table 1.
tic plane) for the calculation of quasi-diffuse fluxes. This The dataset is dominated by bundles of atmospheric
assumption also seems reasonable looking at the weight muons produced in cosmic-ray air shower interactions for
distribution of sources (equal weighting) in figures 9 (a)– tracks coming from the southern hemisphere (θ < 90◦ ).
(e), appendix D. Figure 2 shows the overlap between the Tracks from the northern hemisphere (θ ≥ 90◦ ) origi-
samples. The LSP-BL Lac, FSRQ and ISP+HSP sam- nate mostly from atmospheric neutrino interactions that
ple are nearly independent, with a small overlap of 3 produce muons. In order to reduce the overwhelming
sources between the FSRQs and ISP+HSP samples. The background of direct atmospheric muons to an accept-
largest overlap exists between the FSRQ and LSP sam- able level, it is necessary to impose a high-energy cut
ples, which share around 60% of their sources. The all- for events from the southern hemisphere. The cut raises
blazar sample contains 167 sources that are not shared the effective neutrino energy threshold to approximately
with any sub sample. These are sources that are either 100 TeV (Aartsen et al. 2014b), reducing the sensitivity
unclassified or only classified as BL Lac objects with no 60 We neglect track-like signals from ν + N → τ + X → µ + X,
τ
corresponding synchrotron peak classification. i.e. muons as end products of a ντ charged-current interaction
chain. The τ → µ decay happens with a branching fraction of
3. DATA SELECTION only 17% (Olive 2014), and the additional decay step lowers the
IceCube is a neutrino telescope located at the geo- outgoing muon energy, leading to an even further suppression of the
ντ contribution in a sample of track-like events. For hard fluxes
graphic South pole. It consists of about one km3 of (spectral index 1-2) above PeV energies, where the ντ -influence
Antarctic ice that is instrumented with 5160 optical becomes measurable due to ντ -regeneration (Bugaev et al. 2004),
photo sensors which are connected via cables (“strings”) this treatment is conservative.
All 2LAC Blazars LSP-BL Lacs
Figure 1. Distribution of sources in the sky for the largest and smallest sample of blazars (in equatorial Mollweide projection) — (left)
The largest sample, all 2LAC blazars (862 sources) — (right) The smallest sample, LSP-BL Lacs (68 sources). The excluded region of the
catalogue (|b| ≤ 10◦ ) is highlighted in red.
et al. (2014b).
All 2LAC Blazars 4. ANALYSIS
LSP-BL Lac
4.1. The likelihood function for unbinned ML stacking
167
The analysis is performed via an extended unbinned
maximum likelihood fit (Barlow 1990). The likelihood
19 function consists of two PDFs, one PDF B(x) for a back-
86
68 221 3 298 ground hypothesis and one PDF S(x) for a signal hy-
pothesis. Requiring the total number of observed events
to be the sum of the signal and background events, the
log-likelihood function can be written as
N n
s
X
ln(L){ns , ΓSI } = ln · S(δi , RAi , σi , εi ; ΓSI )
LSP ISP and HSP i=1
N (1)
FSRQ ns
+ 1− · B(sin(δi ), εi ) ,
Figure 2. Visualization of the source overlap between the differ- N
ent blazar populations.
where i indexes individual neutrino events. The likeli-
hood function depends on two free parameters: the nor-
to neutrino sources in this region by at least 1 order of
malization factor ns and spectral index ΓSI of the total
magnitude for spectra softer than E−2 . Only for harder
blazar signal. For computational reasons we assume that
spectra, the southern sky has a significant contribution to
each source of a given population shares the same spec-
the overall sensitivity. The northern sky does not require
tral index. The background evaluation for each event
such an energy cut, as upgoing tracks can only originate
depends on the reconstructed declination δi and the re-
from neutrino interactions, which have a much lower inci-
constructed muon energy εi . The signal part addition-
dence rate. However, at very high energies (again around
ally depends on the reconstructed right ascension RAi ,
100TeV), the Earth absorbs a substantial fraction of neu-
the angular error estimator σi and the power-law spectral
trinos, reducing also the expected astrophysical signal.
index ΓSI .
Charged-current νµ -interactions can happen far outside
The background PDF is constructed from binning the
the instrumented volume and still be detected, as high-
recorded data in reconstructed declination and energy.
energy muons may travel several kilometers through the
It is evaluated as
glacial ice before entering the detector. This effect in-
creases the effective detection area for certain arrival di- 1
B(sin(δi ), εi ) = · f (sin(δi ), εi ), (2)
rections, mostly around the horizon. 2π
The most sensitive region is therefore around the ce- 1
where 2π arises from integration over the right ascension
lestial equator, which does not require a high energy cut,
and f is the normalized joint probability distribution of
provides ample target material surrounding the detec-
the events in declination sin(δ) and energy ε.
tor, i.e. a large effective area, and does not suffer from
The signal PDF that describes a given blazar popula-
absorption of neutrinos above 100 TeV. However, these
tion is a superposition of the individual PDFs for each
zenith-dependent sensitivity changes are mostly impor-
source,
tant for the interpretation of the results (see e.g. section
5.3). The likelihood approach takes these differences into S(δi , RAi , σi , εi ; ΓSI )
account with the ”acceptance” term in eq. (6), section PNsrc
4.1, and a separation into several zenith-dependent anal- j=1 wj · Sj (δi , RAi , σi , εi ; ΓSI ) (3)
= ,
yses is not necessary. For more details on the properties PNsrc
j=1 wj
of the datasets and the zenith-dependent sensitivity be-
haviour, we refer to Aartsen et al. (2013b) and Aartsen where wj is a weight determining the relative normaliza-
tion of the PDF Sj for source j. This weight therefore and E > 100 GeV.
accounts for the relative contribution of source j to the 100GeV
dφγ,j
Z
combined signal. In general, different choices of wj are wj,model = Eγ dEγ (7)
possible. The two choices used in this work are discussed 100MeV dEγ
in section 4.2. Each term Sj in equation 3 is evaluated
as This is motivated by the fact that a similar amount of
energy is channeled into the neutrino and γ-ray emission
Sj (δi , RAi , σi , εi ; ΓSI ) if pion decay from pp or pγ interactions dominates the
1 1

Ψij [δi , RAi ]
2 !
(4) high-energy interaction. While the source environment
= ·exp · · gj (εi ; ΓSI ), is transparent to high-energy neutrinos, it might not be
2πσi 2 2 σi for γ-rays. Reprocessing of γ-rays due to γγ interac-
tions might then shift the energies of the photons to GeV
where the spatial term is expressed as a 2D symmetric and sub-GeV energies before they can leave the sources,
normal distribution and gj is the normalized PDF for the which would make them detectable by the Fermi-LAT.
reconstructed muon energy for source j. The term Ψij This might even be expected in pγ scenarios (Murase
is the angular separation between event i and source j. et al. 2015). Since a large fraction of blazars are located
at high redshifts z ≥ 1 61 , this reprocessing will also take
4.2. Weighting Schemes
place during propagation of the photons in the extra-
The term wj in equation 3 parametrizes the relative galactic background light (EBL), shifting γ-ray energies
contribution of source j to the combined signal. It cor- below a few hundred GeV for such sources (Domı́nguez
responds to the expected number of events for source j, et al. 2013). This places them potentially again in the
which can be expressed as energy range of the Fermi-LAT 2LAC catalogue. Even
Z Eν,max in the case that synchrotron contributions (e.g. muon
wj = Φ0,j · hj (Eν ) · Aeff (θj , Eν ) dEν , (5) or pion synchrotron radiation) dominate over pion de-
Eν,min cay in the MeV-GeV range, which has been considered
in particular for BL Lac objects (Mücke et al. 2003), one
where Aeff (θj , Eν ) is the effective area for incoming muon would expect the overall γ-ray emission to be propor-
neutrinos from a given source direction at a given energy, tional to the neutrino emission. This is also the case in
hj (Eν ) denotes the normalized neutrino energy spectrum models where inverse Compton processes dominate the
for source j, and Φ0,j its overall flux normalization. The high-energy γ-ray emission (Murase et al. 2014).
integration bounds Eν,min and Eν,max are set to 102 GeV The preceding arguments in favour of a γ-weighting
and 109 GeV respectively, except for the differential anal- scheme assume that all sources show equal proportion-
ysis (see section 4.3), in which they are defined for the ality. On a source-by-source basis, however, the propor-
given energy band. tionality factor can vary, as already mentioned in section
Under the assumption that all sources share the same 1.
spectral power-law shape, wj further simplifies via One contributing factor is the fact that Fermi probes
"Z # different sections of the blazar γ-ray peak for each source
Eν,max relative to the peak position. For simplicity, we do not
wj = [Φ0,j ] · h(Eν ; ΓSI ) · Aeff (θj , Eν ) dEν perform a spectral source-by-source fit in this paper,
Eν,min
leaving this aspect for potential future work. This is
= [C · wj,model ] · [wj,acc. (θj , ΓSI )] , (6) also mostly an issue for the ”All 2LAC-Blazar” sample,
since the other sub-classifications described in section 2
and splits into a “model” term, where wj,model is pro- depend on the peak position and this effect is largely mit-
portional to the expected relative neutrino flux of source igated. There are additional reasons for source-by-source
j, and into an “acceptance” term, which is fixed by the fluctuations in the γ/ν correlation due to EBL reprocess-
position of the source and the global energy spectrum. ing. First, the EBL absorption might not be sufficient
The term wj,model is not known, and its choice defines for close-by sources, such that emerging high-energy γ-
the “weighting scheme” for the stacking analysis. The rays are not reprocessed into the energy range of the
following two separate weighting schemes are used for 2LAC catalogue which ends at 100 GeV. Second, EBL
the signal PDF in the likelihood analysis, leading to two reprocessing differs between sources depending on the
different sets of tests. line-of-sight magnetic fields which deflect charged par-
ticle pairs produced in EBL cascades (Aharonian et al.
4.2.1. γ-weighting 1994) differently. Third, strong in-source γγ reprocessing
For this weighting scheme we first have to assume that could lead to γ-rays at even lower energies than 100 MeV
the γ-ray flux can be modeled as being quasi-steady be- (Murase et al. 2015) which would be below the 2LAC en-
tween 2008 and 2010, the time period which forms the ergy range.
basis for the 2LAC catalog. This makes it possible to All results presented in section 5 making use of the
extrapolate the flux expectation of each source to other γ-weighting scheme assume that the potential source-
time periods, e.g. into the non-overlapping part of the to-source fluctuations in the γ − ν correlation described
data-taking period of the IceCube data for this analy- here average out for large source populations and can
sis (2009-2012). Each model weight, i.e. the relative be neglected. More information on the distribution of
neutrino flux expected to arrive from a given source, is
then given by the source’s γ-ray energy flux observed by 61 With the exception of HSP objects, see Ackermann et al.
Fermi-LAT in the energy range between E > 100 MeV (2011).

the weights in dependence of declination can be found in 4.4. Simulations

figures 9 (a)–(e), appendix D. We estimate the sensitivity of our searches in both
weighting schemes using an ensemble of simulated
4.2.2. Equal weighting skymaps containing both background and signal events.
The γ-weighting scheme is optimal under the assump- We simulate the background by drawing events from
tion that the neutrino flux follows the measured γ-energy the experimental data sample, then randomizing their
flux exactly. Given the the uncertainties discussed in sec- right ascensions to remove any correlation with blazar
tion 4.2.1, we also use another weighting scheme, positions. This is the same method used in previous Ice-
Cube point source searches (Aartsen et al. 2013b, 2014b)
wj,model = 1, (8) and mitigates systematic uncertainties in the background
description due to the data-driven event injection.
which we expect to be more sensitive eventually if the The injection for signal differs depending on the
actual γ − ν correlation varies strongly from source weighting schemes. For the γ-weighting scheme, we in-
to source. It provides a complementary and model- ject signal events with the relative flux contribution of
independent test in which we are maximally agnostic to each source determined by the weight factors wj,model
the degree of correlation between γ-ray and neutrino lu- that are used in the PDF. In the equal weighting scheme,
minosities. following the same approach would lead to a simulated
We do not assume a specific neutrino emission in a signal of n equally bright sources, which is not realistic
given source when calculating the flux upper limits for for a population distributed widely in redshift and lumi-
the equal weighting scheme, in particular no equal emis- nosity. Therefore we inject events using a relative neu-
sion. We only assume, to some approximation, that the trino flux contribution that follows a realistic SCD. Since
differential source count distributions (SCDs) of γ-rays the neutrino dN/dS distribution of blazars is unknown,
and neutrinos have comparable shapes. The differen- we have chosen to use the blazar γ-ray SCD published
tial source count distribution, dN/dS, describes how the in Abdo et al. (2010c) as a template62 . Here we assume
energy-integrated flux S is distributed over all sources that for the population under investigation, the relative
and is a crucial property of any cosmological source pop- contributions to the total neutrino flux are distributed in
ulation. Section 4.4 provides more information on the a similar fashion as the relative contributions to the total
technical aspects of neutrino flux injection in the equal γ-ray flux. However, there are no assumptions about the
weighting test. Appendix A then discusses why the correlation of the neutrino and γ-ray flux for individual
methodology is robust against variations in the actual sources.
shape of the dN/dS distribution for the neutrino flux in There are two reasons to choose the γ-ray SCD as the
the IceCube energy range, and why the final result is primary template for the shape of the neutrino SCD. The
valid even if the neutrino SCD is different from the γ-ray first is that we select the populations based on their γ-ray
SCD. emission to start with. The second is that the form of the
high-energy γ-ray SCD is quite general, and has also been
4.3. Statistical tests observed with AGN detected in the radio (Hopkins et al.
2003) and x-ray (Georgakakis et al. 2008) bands. It starts
We perform statistical tests for each population of
blazars. The log-likelihood difference λ defines our test with quasi-Euclidean behavior (S 5/2 · dN/dS ≈ const.)
statistic (TS), given by at high fluxes, and then changes to a harder power law
index towards smaller flux values which ensures that the
λ = − 2 · log(L){ns = 0} total flux from the population remains finite.
(9) The skymap simulations are performed for many pos-
+ 2 · log(L){ns = ns,max , ΓSI = ΓSI,max }, sible SCD realizations by sampling from the dN/dS dis-
where ns,max and ΓSI,max are the number of signal events tribution. This is necessary since the number of sig-
and the signal spectral index that maximize the TS. nal events expected in IceCube for a given neutrino flux
We simulate an ensemble of background-only skymaps, varies greatly over the two hemispheres (see section 3).
where the TS distribution is compared with the TS value Thus, it matters how the neutrino flux is distributed over
obtained from the data. The p-value is then defined as the individual sources for the value of the resulting con-
the fraction of skymaps in the background ensemble that fidence interval. The shape of the SCD and the flux
has a larger TS value than the one observed. Ensembles sampling range have an additional impact. See appendix
of skymaps with different injected signal strengths are A for further details in the context of confidence interval
then used to calculate the resulting confidence interval. construction.
See section 4.4 for details on the skymap simulations.
5. RESULTS
In total we perform two distinct types of tests for which
p-values are calculated. The first (“integral”) assumes a 5.1. Observed p-values
power law spectrum for the blazar emission over the full Table 2 summarizes p-values for the “integral” test
energy range observable with IceCube (unless stated oth- (see section 4.3). Nine out of the ten tests show over-
erwise). The second (“differential”) assumes a neutrino fluctuations, but no significant excess. We find the
signal that is confined to a small energy range (half a
decade in energy) and has a power law spectrum with a 62 This blazar SCD strictly stems from the 1FGL catalogue
spectral index of −2 within this range. We perform the (Abdo et al. 2010b), but any SCD based on a newer catalog is
differential test for 14 energy ranges between 100 GeV not expected to change significantly since a large fraction of the
and 1 EeV. total γ-ray flux is already resolved in the 1FGL.
p-value of the background which is used for the construction of
Population
γ-weighting equal weighting differential flux upper limits.
All 2LAC blazars 36% (+0.4σ) 6% (+1.6σ) We give all further results in intensity units and calcu-
FSRQs 34% (+0.4σ) 34% (+0.4σ) late the quasi-diffuse flux63 for each population.
LSPs 36% (+0.4σ) 28% (+0.6σ) The flux upper limits in the equal weighting scheme
ISP/HSPs > 50% 11% (+1.2σ)
are calculated using multiple samplings from an assumed
LSP-BL Lacs 13% (+1.1σ) 7% (+1.5σ)
neutrino SCD for the blazars, as already outlined in sec-
Table 2 tion 4.4. Please refer to appendix A for further details
P-values and the corresponding significance in units of standard about the dependence of the flux upper limit on the
normal deviations in the power-law test. The table shows the choice of the SCD and a discussion of the robustness of
results for both weighting schemes.The values do not include a the equal-weighting results. In general, the equal weight-
trial factor correction.
ing upper limit results do not correspond to a single flux
value, but span a range of flux values.
For each upper limit64 we determine the valid energy
100 range according to the procedure in appendix B. This
0σ energy range specifies where IceCube has exclusion power
local p-value
for a particular model, and is also used for visualization
significance
10 −1 1σ
purposes in all figures.
2σ Systematic effects influencing the upper limits are
10−2 dominated by uncertainties on the absorption and scat-
tering properties of the Antarctic ice and the detection
10−3 3σ efficiency of the optical modules. Following Aartsen et al.
(2014b), the total systematic uncertainty on the upper
limits is estimated to be 21%. Since we are dealing with
10−4 2 3 4 5 6 7 8 9 upper limits only, we conservatively include the uncer-
10 10 10 10 10 10 10 10
tainty additively in all figures and tables.
Neutrino Energy [GeV]
5.3. Generic upper limits
Figure 3. Local p-values for the sample containing all 2LAC
blazars using the equal-weighting scheme (black) and γ-weighting Table 3 shows flux upper limits assuming a generic
scheme (green) in the differential test. power-law spectrum for the tested blazar populations,
strongest over-fluctuation, a 6% p-value, using the equal- calculated for the three different spectral indices −1.5,
weighting scheme for all 2LAC blazars. We omit a trial- −2.0, and −2.7.
factor correction because the populations have a large The distribution of the γ-ray energy flux among the
overlap and the result is not significant. sources in each population governs the flux upper limit
Figure 3 shows the p-values from the corresponding in the γ-weighting scheme. It is mostly driven by the dec-
“differential” test. The largest excess is visible in the lination of the strongest sources in the population, due to
5−10 TeV energy band with a pre-trial p-value of 4·10−3 . the strong declination dependence of IceCube’s effective
This outcome is totally compatible with a fluctuation of area. For FSRQs, the two sources with the largest γ-
the background, since the effect of multiple trials has to weights (3C 454.3 at DEC2000 = 16◦ and PKS1510-08 at
be taken into account which reduces the significance of DEC2000 = −9◦ ) carry around 15% of the total γ-weight
the observation substantially. An accurate calculation of of all FSRQs. Their positions close to the equator place
the trial-corrected p-value is again difficult, as neither the them in the most sensitive region for the IceCube detec-
five blazar samples, nor the 14 tested energy ranges per tor, and the γ-weighting upper limits for FSRQs are more
sample, are independent. We again omit it for simplicity. than a factor of 2 lower than the corresponding equal-
Comparing the differential p-value plot of all 2LAC weighting limits. For the LSP-BL Lacs, the two strongest
blazars with the other populations (see figures 10 (a)– sources (PKS 0426-380 at DEC2000 = −38◦ and PKS
(e) in appendix D), one finds that the overfluctuation is 0537-441 at DEC2000 = −44◦ ) carry nearly 30% of the
caused by the LSP-BL Lac-, FSRQ- and ISP/HSP pop- total γ-weight but are located in the southern sky, where
ulation, which are nearly independent and each show a IceCube is not very sensitive. The γ-weighting upper
small excess in 5 TeV - 20 TeV region. In the γ-weighting limit is therefore comparable to the equal-weighting up-
scheme, the ISP/HSP p-value distribution is nearly flat, per limit. The reader is referred to appendix D for more
which leads to the weaker overfluctuation in the ”all information on the weight distribution.
2LAC blazar” sample compared to the equal weighting Figure 4 shows the differential upper limit in compari-
scenario. son to the median sensitivity for all 2LAC blazars using
the equal-weighting scheme. This population showed the
5.2. Flux upper limits largest overfluctuation. We plot here the upper limit de-
rived from the median SCD sampling outcome, since in
Since no statistically significant neutrino emission from general the equal weighting upper limit depends on the
the analyzed source populations was found, we calculate neutrino flux realization of the SCD (see appendix A).
flux upper limits using various assumptions about their As expected, the differential limit is slightly higher, by
energy spectrum. We use the CLs upper limit construc-
tion (Read 2000). It is more conservative than a standard 63 The flux divided by the solid angle of the sky above 10 degrees
Neyman construction, e.g. used in Aartsen et al. (2014b), galactic latitude, i.e. 0.83 × 4π. See section 2 for a justification.
64 With the exception of the differential upper limit.
but allows for a proper evaluation of under-fluctuations
Spectrum: Φ0 · (E/GeV)−1.5
2LAC Blazar Upper Limit equal weighting
Φ0 90% [GeV−1 cm−2 s−1 sr−1 ] 10−6
[GeVs−1 cm−2 sr−1 ]

Blazar Class ΓSI = −2.5, Eν > 10 TeV γ-weighting
γ-weighting equal weighting ΓSI = −2.2, Eν > 10 TeV
All 2LAC Blazars 1.6 × 10−12 4.6 (3.8 − 5.3) × 10−12
10−7
FSRQs 0.8 × 10−12 2.1 (1.0 − 3.1) × 10−12
LSPs 1.0 × 10−12 1.9 (1.2 − 2.6) × 10−12
ISPs/HSPs 1.8 × 10−12 2.6 (2.0 − 3.2) × 10−12 10−8
LSP-BL Lacs 1.1 × 10−12 1.4 (0.5 − 2.3) × 10−12 19% − 27%
νµ
10−9
dΦ
E 2 dE
Spectrum: Φ0 · (E/GeV)−2.0
Φ0 90% [GeV−1 cm−2 s−1 sr−1 ] Astrophysical Diffuse Flux
7%
Blazar Class 10−10
γ-weighting equal weighting
102 10 3
104
10 5
10 106 7
108 109
All 2LAC Blazars 1.5 × 10−9 4.7 (3.9 − 5.4) × 10−9
FSRQs 0.9 × 10−9 1.7 (0.8 − 2.6) × 10−9
LSPs 0.9 × 10−9 2.2 (1.4 − 3.0) × 10−9
Figure 5. 90% C.L. flux upper limits for all 2LAC blazars in
ISPs/HSPs 1.3 × 10−9 2.5 (1.9 − 3.1) × 10−9 comparison to the observed astrophysical diffuse neutrino flux. The
LSP-BL Lacs 1.2 × 10−9 1.5 (0.5 − 2.4) × 10−9 latest combined diffuse neutrino flux results from Aartsen et al.
(2015b) are plotted as the best-fit power-law with spectral index
Spectrum: Φ0 · (E/GeV)−2.7 −2.5 , and as a differential flux unfolding using 68% central and
Φ0 90% [GeV−1 cm−2 s−1 sr−1 ] 90% U.L. confidence intervals. The flux upper limit is shown using
Blazar Class
γ-weighting equal weighting both weighting schemes for a power-law with spectral index −2.5
(blue). Percentages denote the fraction of the upper limit compared
All 2LAC Blazars 2.5 × 10−6 8.3 (7.0 − 9.7) × 10−6 to the astrophysical best fit value. The equal-weighting upper limit
FSRQs 1.7 × 10−6 3.3 (1.6 − 5.1) × 10−6 for a flux with a harder spectral index of −2.2 is shown in green.
LSPs 1.6 × 10−6 3.8 (2.4 − 5.2) × 10−6
teresting to observe this excess with future IceCube data.
ISPs/HSPs 1.6 × 10−6 4.6 (3.5 − 5.6) × 10−6
For information on the differential upper limits from the
LSP-BL Lacs 2.2 × 10−6 2.8 (1.0 − 4.6) × 10−6
other samples the reader is referred to appendix D.
Table 3
90% C.L. upper limits on the diffuse (νµ + ν µ )-flux from the 5.4. The maximal contribution to the diffuse
different blazar populations tested. The table contains results for astrophysical flux
power-law spectra with spectral indices −1.5, −2.0, and −2.7.
The equal-weighting column shows the median flux upper limit The astrophysical neutrino flux is observed between
and the 90% central interval of different sample realizations of the 10 TeV and 2 PeV (Aartsen et al. 2015b). Its spectrum
Fermi-LAT source count contribution (in parentheses). All values has been found to be compatible with a single power-law
include systematic uncertainties.
and a spectral index of −2.5 over most of this energy
range. Accordingly, we use a power-law with the same
spectral index and a minimum neutrino energy of 10 TeV
10−5
for the signal injected into the simulated skymaps when
Sensitivity
calculating the upper limit for a direct comparison. Fig-
10−6 ±1σ
ure 5 shows the flux upper limit for an E −2.5 power-law
±2σ
spectrum starting at 10 TeV for both weighting schemes
10−7 in comparison to the most recent global fit of the astro-
physical diffuse neutrino flux, assuming an equal compo-
10−8 sition of flavors arriving at Earth.
The equal-weighting upper limit results in a maximally
νµ
19%-27% contribution of the total 2LAC blazar sample

dΦ
E 2 dE
10−9 All 2LAC blazars (equal weighting)

to the observed best fit value of the astrophysical neu-
90% C.L. Upper Limit (median SCD outcome)
10−10 2 trino flux, including systematic uncertainties. This limit
10 103 104 105 106 107 108 109 is independent of the detailed correlation between the
Neutrino Energy [GeV] γ-ray and neutrino flux from these sources. The only as-
sumption is that the respective neutrino and γ-ray SCDs
Figure 4. Differential 90% C.L. upper limit on the (νµ + ν µ )-flux have similar shapes (see section 5.2 for details on signal
using equal weighting for all 2LAC blazars. The ±1σ and ±2σ injection). We use the Fermi-LAT blazar SCD as pub-
null expectation is shown in green and yellow, respectively. The
upper limit and expected regions correspond to the median SCD lished in Abdo et al. (2010c) as a template for sampling.
sampling outcome. However, we find that even if the shape of the SCD dif-
fers from this template, the upper limit still holds and
a factor of about 2, than the median outcome in the en- is robust. In appendix A we discuss the effect of differ-
ergy range between 5 TeV and 10 TeV where the largest ent SCD shapes and discuss how the combination with
excess is observed. This is the average behavior for a soft existing point source constraints (Aartsen et al. 2015c)
flux with spectral index of about −3.0 65 , if one assumes leads to a nearly SCD-independent result, since a point
a simple power-law fit to explain the data. While such a source analysis and a stacking search with equal weights
physical interpretation can not be made yet, it will be in- effectively trace opposite parts of the available parameter
space for the dN/dS distribution.
65 This can be read off in figure 8. The ratio function indicates in
In case we assume a proportionality between the γ-ray
which energy range a given flux function appears first, on average. and neutrino luminosities of the sources, the γ-weighting
Population
weighting scheme (see appendix B, figure 8) by the hard component above
equal γ γ (extrapol.) a PeV is negligible. In such a scenario we expect the con-
all 2LAC blazars 19% − 27% 7% 10% straint to be rather similar to our result from the simple
FSRQs 5% − 17% 5% 7% power-law test with spectral index −2.5.
LSPs 6% − 15% 5% 7%
ISP/HSPs 9% − 15% 5% 7% 5.5. Upper limits on models for diffuse neutrino
LSP-BL Lacs 3% − 13% 6% 9% emission
Table 4
For the experimental constraints on existing theoreti-
Maximal contributions to the best-fit diffuse flux from Aartsen cal calculations, we only considered models for the dif-
et al. (2015b) assuming equipartition of neutrino flavours. The fuse emission from blazar populations, not predictions
equal-weighting case shows this maximal contribution for the 90% for specific objects. These include the calculations by
central outcomes of potential dN/dS realizations. The last column
shows the maximal contribution of the integrated emission from Mannheim (1995), Halzen & Zas (1997) and Protheroe
the total parent population in the observable universe exploiting (1997) for generic blazars, the calculations by Becker
the γ-ray completeness of the 2LAC blazars (see appendix C). et al. (2005) and Murase et al. (2014) for FSRQs and
calculations by Mücke et al. (2003),Tavecchio et al.
(2014),Tavecchio & Ghisellini (2015) and Padovani et al.
limit constrains the maximal flux contribution of all (2015) for BL Lacs.
2LAC blazars to 7% of the observed neutrino flux in the The upper limits in this section are calculated using
full 10 TeV to 2 PeV range. Since the blazars resolved in the γ-weighting scheme and therefore assume a correla-
the 2LAC account for 70 % of the total γ-ray emission tion between the neutrino flux and the measured γ-ray
from all GeV blazars (Ajello et al. 2015)this further im- energy flux. This allows us to account for the fraction of
plies that at most 10% of the astrophysical neutrino flux the neutrino emission that arises from the blazars not de-
stems from all GeV blazars extrapolated to the whole tected in γ-rays. The fraction of γ-ray emission from re-
universe, again in the full 10 TeV to 2 PeV range and solved 2LAC blazars in general (including BL Lacs), and
assuming the γ-weighting is an appropriate weighting of FSRQs in particular, is about 70% (Ajello et al. 2015,
assumption. Table 4 summarizes the maximal contribu- 2012). Therefore, the flux upper limits for the entire
tions for all populations, including the γ-weighting result population are a factor 1/0.7 ≈ 1.43 weaker than those
scaled to the total respective total population of sources derived for the quasi-diffuse flux of the 2LAC blazars.
in the observable universe. See appendix C for more details on this factor.
It is interesting to compare these numbers directly to Table 5 summarizes model rejection factors (Hill &
the γ-ray sector. Ajello et al. (2015) show that GeV Rawlins 2003)66 for all considered models. Many of these
blazars (100 MeV − 100 GeV) contribute approximately models can be constrained by this analysis. Figure 6 a)-
50% to the extragalactic gamma-ray background (EGB). d) visualizes the flux upper limits in comparison to the
The resolved 1FGL (Abdo et al. 2010b) blazar compo- neutrino flux predictions.
nent in particular contributes around 35%. This estimate In the early models (before the year 2000) the neu-
should be rather similar for the 2LAC blazars studied trino flux per source is calculated as being directly pro-
here, which are defined based on the more recent 2FGL portional to the γ-ray flux in the energy range Eγ >
catalogue (Nolan et al. 2012) (see appendix C for a dis- 100 MeV(Mannheim 1995) (A), Eγ > 1 MeV (Mannheim
cussion). The 2LAC blazar contribution to the astro- 1995) (B), 20 MeV < Eγ < 30 GeV (Halzen & Zas 1997)
physical neutrino flux is therefore by at least a factor and Eγ > 100 MeV (Protheroe 1997). The γ-weighting
0.75 smaller than the corresponding extragalactic con- scheme is therefore almost implicit in all these calcu-
tribution in the γ-regime. The difference of this contri- lations, although the energy ranges vary slightly from
bution between the two sectors becomes substantial (7% the 100 MeV − 100 GeV energy range used for the γ-
maximally allowed contribution for neutrinos versus 35% weighting.
for γ-rays) if one assumes a γ/ν-correlation. From the newer models, only Padovani et al. (2015)
Figure 5 also shows the equal-weighting constraint for uses a direct proportionality between neutrino and γ-ray
a harder neutrino spectrum with a spectral index of −2.2. flux (for Eγ > 10 GeV), where the proportionality fac-
This harder spectral index is about 3 standard deviations tor encodes a possible leptonic contribution. In all other
away from the best fit value derived in Aartsen et al. publications a direct correlation to γ-rays is not used
(2015b), and can be used as an extremal case given the for the neutrino flux calculation. Since all these models
current observations. The comparison of this upper limit assume that p/γ-interactions dominate the neutrino pro-
with the hard end of the ”butterfly” shows that even in duction, the resulting neutrino fluxes are calculated via
this case less than half of the bulk emission can originate the luminosity in the target photon fields. In Becker et al.
in the 2LAC blazars with minimal assumptions about (2005) the neutrino flux is proportional to the target ra-
the relative neutrino emission strengths. Due to the low- dio flux which in turn is connected to the disk luminosity
count status of the data, we omit multi power-law spectra via the model from Falcke & Biermann (1995). In Mücke
tests at this point. However, one can estimate the con- et al. (2003) it is directly proportional to the radiation
straints for more complicated models using figure 8 in of the synchrotron peak. In Murase et al. (2014) the
appendix B, which shows the energy range for a given neutrino flux is connected to the x-ray luminosity, which
spectrum that contributes the dominant fraction to the in turn is proportional to the luminosity in various tar-
sensitivity. The sensitivity for a possible two-component get photon fields. In Tavecchio et al. (2014) the neu-
model, for example, having a soft component at TeV en-
ergies and a hard component in the PeV range, would 66 The flux upper limit divided by the flux predicted in the
be dominated by the soft regime, as the ”ratio function” model.

Type Model MRF

(A) 1.30
(Mannheim 1995)
Generic blazars < 0.1 (B)
(Halzen & Zas 1997) < 0.1
(Protheroe 1997) < 0.1
(Becker et al. 2005) 2.28
ΓSI = −2.0 (BLR) ξCR < 12
FSRQs ΓSI = −2.0 (blazar) ξCR < 21
(Murase et al. 2014) ΓSI = −2.3 (BLR) ξCR < 153
ΓSI = −2.3 (blazar) ξCR < 241
HSP (optimistic) 76.29
(Mücke et al. 2003)
LSP (optimistic) 5.78
BL Lacs HSP-dominated (1) 1.06
(Tavecchio et
a al. 2014) HSP-dominated (2) 0.35
(Tavecchio & Ghisellini 2015) LSP-dominated 0.21
(Padovani et al. 2015) HSP (baseline) 0.75
a Predictions from Tavecchio et al. (2014); Tavecchio & Ghisellini (2015) enhanced
by a factor 3 in correspondence with the authors.
Table 5
Summary of constraints and model rejection factors for the diffuse neutrino flux predictions from blazar populations. The values include a
correction factor for unresolved sources (see appendix C) and systematic uncertainties. For models involving a range of flux predictions
we calculate the MRF with respect to the lower flux of the optimistic templates (Mücke et al. 2003) or constraints on baryon to photon
luminosity ratios ξCR (Murase et al. 2014).
10−5 10−4 Models 90% C.L.

Astrophysical Upper Limit
−6 Mücke et al. LSP (optimistic)

10 Diffuse Flux 10−5
Mücke et al. HSP (optimistic)
10−7 10−6 Tavecchio et al. LSP
Padovani et al. HSP
10−7
10−8
Models 10−8
10−9 Mannheim (A,B)
(B) (A) 10−9
10−10 Protheroe
νµ
νµ
10−10
dΦ
E 2 dE
Halzen
dΦ
E 2 dE
10−11 10−11
90% C.L. Differential Upper Limit
10−12 2 10−12 2
10 103 104 105 106 107 108 109 1010 1011 1012 10 103 104 105 106 107 108 109 1010 1011 1012
Neutrino Energy [GeV] Neutrino Energy [GeV]
(a) generic blazars (b) BL Lacs

10−5 10−5
Models Models
10−6 Murase BLR-Zone, Γp = 2.0 10−6 Murase BLR-Zone, Γp = 2.3

Murase Blazar-Zone, Γp = 2.0 Murase Blazar-Zone, Γp = 2.3
10−7 10−7 Becker ‘05
−8 −8
10 10
−9
10 10−9
10−10 10−10
νµ
νµ
dΦ
dΦ
E 2 dE
E 2 dE
10−11 10−11
90% C.L. Upper Limit 90% C.L. Upper Limit
10−12 2 10−12 2
10 103 104 105 106 107 108 109 1010 1011 1012 10 103 104 105 106 107 108 109 1010 1011 1012
(c) FSRQs - 1 (d) FSRQs - 2

Figure 6. 90% C.L. upper limits on the (νµ + ν µ )-flux for models of the neutrino emission from (a) generic blazars (Mannheim 1995;
Halzen & Zas 1997; Protheroe 1997), (b) BL Lacs (Mücke et al. 2003; Tavecchio & Ghisellini 2015; Padovani et al. 2015) and (c)+(d) FSRQs
(Becker et al. 2005; Murase et al. 2014). The upper limits include a correction factor that takes into account the flux from unresolved
sources (see appendix C) and systematic uncertainties. The astrophysical diffuse neutrino flux measurement (Aartsen et al. 2015b) is shown
in green for comparison.
trino luminosity is calculated using target photon fields ality between the γ-ray luminosity in the 2LAC
from the inner jet “spine-layer”. However, a correlation energy range and the neutrino luminosity, we can
to the γ-ray flux in these latter models may still exist, extend the constraint to the parent population of
even in the case that leptonic γ-ray contributions dom- all GeV blazars in the observable universe. The
inate. This is mentioned in Murase et al. (2014), which corresponding maximal contribution from all GeV
explicitly predicts the strongest γ-ray emitters to be also blazars is then around 10%, or 5 − 10% from the
the strongest neutrino emitters, even though the model other blazar sub-populations. In each case we use
contains leptonically produced γ-ray emission. It should the same power-law assumption as before in or-
be noted that an independent IceCube analysis studying der to compare to the observed flux. For FSRQs
the all-flavor diffuse neutrino flux at PeV energies and our analysis allows for a 7% contribution to the
beyond (Aartsen et al. 2016a) recently also put strong diffuse flux, which is in general agreement with a
constraints on some of the flux predictions discussed in result found by Wang & Li (2015) who indepen-
this section. dently estimated that FSRQs do not contribute
more than 10% to the diffuse flux using our earlier
6. SUMMARY AND OUTLOOK small-sample stacking result for 33 FSRQs (Aart-
In this paper, we have analyzed all 862 Fermi-LAT sen et al. 2014b).
2LAC blazars and 4 spectrally selected sub populations
via an unbinned likelihood stacking approach for a cu- 2. We calculated upper limits using the γ-weighting
mulative neutrino excess from the given blazar direc- scheme for 15 models of the diffuse neutrino emis-
tions. The study uses 3 years of IceCube data (2009- sion from blazar populations found in the litera-
2012) amounting to a total of around 340000 muon-track ture. For most of these models, the upper limit
events. constrains the model prediction, for some of them
Each of the 5 populations were analyzed with two by more than an order of magnitude. The implicit
weighting schemes which encode the assumptions about assumption in all these upper limits is a proportion-
the relative neutrino flux from each source in a given ality between the source-by-source γ-ray luminosity
population. The first weighting scheme uses the energy in the 2LAC energy range and its corresponding
flux observed in γ-rays as weights, the second scheme neutrino luminosity. All models published before
gives each source the same weight. This resulted in a the year 2000, and the model by Padovani et al.
total of 10 statistical tests which were in turn analyzed (2015) implicitly contain this assumption, although
in two different ways. The first is an “integral” test, in some of their energy ranges differ from the exact
which a power-law flux with a variable spectral index is energy range in the 2LAC catalogue. Even for the
fitted over the full energy range that IceCube is sensi- other models the proportionality assumption may
tive to. The second is a differential analysis, in which 14 still hold, as indicated by Murase et al. (2014).
energy segments between 102 GeV and 109 GeV, each
spanning half a decade in energy, are fit independently Kadler et al. (2016) recently claimed a 5% chance prob-
with a constant spectral index of −2. ability for a PeV IceCube event to have originated close
Nine from ten integral tests show over-fluctuations, but to blazar PKS B1424-418 during a high-fluence state.
none of them are significant. The largest overfluctua- While 5% is not yet statistical evidence, our results do
tion, a 6% p-value, is observed for all 862 2LAC blazars not contradict such single PeV-event associations, es-
combined using the model-independent equal-weighting pecially since a dominant fraction of the sensitivity of
scheme. The differential test for all 2LAC blazars using our analysis comes from the sub-PeV energy range. The
equal source weighting reveals that the excess appears in same authors also show that the measured all-sky PeV
the 5-10 TeV region with a local p-value of 2.6σ. No cor- neutrino flux can not be compatible with an origin in
rection for testing multiple hypotheses is applied, since a pure FSRQ population that has a peaked spectrum
even without a trial correction this excess cannot be con- around PeV energies, as it would over-predict the num-
sidered significant. ber of observed events. Instead, one has to invoke ad-
Given the null results we then calculated flux upper ditional assumptions, for example a certain contribution
limits. The two most important results of this paper are: from BL Lacs, leptonic contributions to the SED, or a
spectral broadening of the arriving neutrino flux down to
1. We calculated a flux upper limit for a power-law TeV energies due to Doppler shifts from the jets and the
spectrum starting at 10 TeV with a spectral index intrinsic redshift distribution of the blazars. Our results
of −2.5 for all 2LAC blazars. We compared this suggest that the last assumption, a spectrum broaden-
upper limit to the diffuse astrophysical neutrino ing down to TeV energies, only works if the resulting
flux observed by IceCube (Aartsen et al. 2015b). power-law spectral index is harder than around −2.2,
We found that the maximal contribution from all as the flux is otherwise in tension with our γ-weighting
2LAC blazars in the energy range between 10 TeV upper limit. A hard PeV spectrum is interestingly also
and 2 PeV is at most 27%, including systematic seen by a recent IceCube analysis (Aartsen et al. 2016b)
effects and with minimal assumptions about the that probes the PeV range with muon neutrinos. Re-
neutrino/γ-ray correlation in each source. Chang- gardless of these speculations, the existing sub-PeV data
ing the spectral index of the tested flux to −2.2, a requires an explanation beyond the 2LAC sample from a
value allowed at about 3 standard deviations given yet unidentified galactic or extra-galactic source class.
the current global fit result (Aartsen et al. 2015b), Our results do not provide a solution to explain the
weakens this constraint by about a factor of two. bulk emission of the astrophysical diffuse neutrinos, but
If we assume for each source a similar proportion- they provide robust constraints that might help to con-
struct the global picture. Recently, Murase et al. (2015) All examples are shown for a population of 862 objects,
argued that current observations favour sources that are the size of the “All 2LAC blazars” sample. In the left
opaque to γ-rays. This would for example be expected panel (“extended”), the SCD template is the measured
in the cores of AGN. Our findings on the 2LAC blazars blazar γ-ray SCD from Abdo et al. (2010c) and is extrap-
mostly probe the emission from relativistically beamed olated five orders of magnitude to lower flux values. The
AGN jets and are in line with these expectations. We minimum flux value is arbitrarily chosen, but it is small
also do not constrain neutrinos from blazar classes that enough such that the distribution is a scale-free power-
are not part of the 2LAC catalogue, for example extreme law and an extension towards even smaller flux values
HSP objects. These sources might emit up to 30% of the does not make a difference in terms of average sample
diffuse flux (Padovani et al. 2016), and studies in this outcomes. In the second panel (“standard”),the SCD
direction with other catalogues are in progress. is exactly the Fermi-LAT blazar SCD from Abdo et al.
While the slight excess in the 5-10 TeV region is not (2010c) and spans three orders of magnitude in flux. In
yet significant, further observations by IceCube may clar- the third panel the SCD is of Euclidean form and ex-
ify if we see an emerging soft signal or just a statistical tended over a flux range such that the cumulative SCD
fluctuation. equals to the number count of the “standard” distribu-
tion in the second panel. In the fourth panel the SCD is
a delta distribution, which gives an equal weight to each
We acknowledge the support from the following agen- source, i.e. the assumption that is used in the weighting
cies: U.S. National Science Foundation-Office of Polar of the PDF for the statistical test. The lower row dis-
Programs, U.S. National Science Foundation-Physics Di- plays the respective 90% central interval for upper limit
vision, University of Wisconsin Alumni Research Foun- outcomes from this analysis and constraints from 6-year
dation, the Grid Laboratory Of Wisconsin (GLOW) grid single point source search discovery potentials.
infrastructure at the University of Wisconsin - Madi- As we sample the relative source contributions from a
son, the Open Science Grid (OSG) grid infrastructure; growing flux range, corresponding to the column order
U.S. Department of Energy, and National Energy Re- 4 → 3 → 2 → 1, flux upper limit variations increase
search Scientific Computing Center, the Louisiana Opti- and the stacking analysis constraint weakens. At the
cal Network Initiative (LONI) grid computing resources; same time, the constraints from the single point source
Natural Sciences and Engineering Research Council search become stronger with a single source more and
of Canada, WestGrid and Compute/Calcul Canada; more dominating the total population.
Swedish Research Council, Swedish Polar Research Sec- The first and fourth columns correspond to limiting
retariat, Swedish National Infrastructure for Comput- cases, neither of which are appropriate to use for this
ing (SNIC), and Knut and Alice Wallenberg Foun- analysis, but are just shown to illustrate the general be-
dation, Sweden; German Ministry for Education and havior of the procedure. The delta-peak SCD (4th col-
Research (BMBF), Deutsche Forschungsgemeinschaft umn) is unphysical, since it corresponds to an equal flux
(DFG), Helmholtz Alliance for Astroparticle Physics per source. The extrapolated SCD (1st column) yields
(HAP), Research Department of Plasmas with Complex an extreme spread of signal contributions, which roughly
Interactions (Bochum), Germany; Fund for Scientific Re- corresponds to a random draw from the entire population
search (FNRS-FWO), FWO Odysseus programme, Flan- in the universe, assuming that the faint end corresponds
ders Institute to encourage scientific and technological re- to the weakest blazars that can in principle be detected.
search in industry (IWT), Belgian Federal Science Policy Since the random draw mostly exists of sources from
Office (Belspo); University of Oxford, United Kingdom; the faint-flux end, it corresponds to a situation where
Marsden Fund, New Zealand; Australian Research Coun- the neutrino flux is anti-proportional to the γ-ray flux,
cil; Japan Society for Promotion of Science (JSPS); the which is unphysical. All results in this paper make use
Swiss National Science Foundation (SNSF), Switzerland; of the 2nd column SCD. The cross-over of point source
National Research Foundation of Korea (NRF); Villum constraints and constraints from this stacking analysis
Fonden, Danish National Research Foundation (DNRF), is reached for realizations drawn from a dN/dS distri-
Denmark bution that lies in between the SCD from the 1st and
2nd column. For illustrative purposes we also include
APPENDIX a particular flux scenario with equal intrinsic luminos-
ity, where the neutrino flux per source is proportional to
DEPENDENCE OF FLUX UPPER LIMITS ON THE
SOURCE COUNT DISTRIBUTION SAMPLING
1/dL 2 67 . Missing redshifts for BL Lacs are drawn ran-
domly from the BL Lac distribution where the redshifts
The equal-weighting limits use source count distribu- are known, taking into account the synchroton peak in-
tions to model the neutrino injection. The SCD serves formation. The results for two such realizations of miss-
as a PDF template from which relative neutrino injec- ing redshift sources lie in the range of SCD draws from
tion weights are drawn. Depending on the shape of the the 2nd column. Since the two realizations do not dif-
SCD, and the range of flux values in which the SCD is fer significantly, we conclude that the resulting estimate
being sampled, the resulting central neutrino upper limit is robust, even though some of the BL Lac redshifts are
value shifts and the range of calculated flux upper limits unknown.
broadens. This is illustrated in figure 7, which shows the
upper limits derived from ensemble simulations drawn
from four different SCDs. It also contains standard 6-
year point source constraints (Aartsen et al. 2015c) for 67 Using standard ΛCDM cosmological parameters as deter-
similar flux realizations. mined by the Planck mission (Planck Collaboration 2015).
7000
107
105
euclidean 6000 ”delta-peak”
106
dN/dlog10 (S) 103 5000
105 104
4000
104
103 3000
103 102
2000
10 2 extended standard 102 1000
0
−14 −10 −6 −8 −7 −6 −8 −7 −6 −8 −7 −6
log10 (S[a.u.])
90 % central interval from dN/dS sample realizations

u.l. this stacking analysis 6y PS disc. pot. constraints
u.l. this analysis - particular realizations

equal lum., wj ∝ 1/dL 2 (1) equal lum., wj ∝ 1/dL 2 (2)
10−4
10−5
10−6
νµ
dΦ
E 2.5 dE
10−7
Figure 7. Comparison of equal-weighting upper limits for different SCDs which are used to sample relative source injection weights,
shown for the population of all 2LAC blazars. The upper row shows the SCDs and the lower row the respective constraints for an E −2.5
flux starting at 10 TeV. The source flux S on the x-axis is shown in arbitrary units, since only the relative neutrino flux counts, but orients
itself by the integrated γ-ray flux from Abdo et al. (2010c). The light blue band marks the 90% central interval of upper limit outcomes
for random samplings of the given SCD, the light red band marks the constraints from the 6-year PS search for similar random samplings.
Two specific realizations modeling equal intrinsic luminosity, taking into account the luminosity distance and randomly drawn redshifts for
missing-redshift BL Lacs, are shown in green and magenta.
DETERMINATION OF ENERGY RANGE FOR UPPER UPPER LIMIT CORRECTION FACTORS FOR
LIMITS UNRESOLVED SOURCES
We determine the energy range for which the upper The γ-weighting scheme implicitly assumes a propor-
limit is supported by IceCube data based on the differ- tionality between neutrino and γ-ray luminosities. In this
ential sensitivity. The procedure is illustrated in figure 8 case, the fraction of the total neutrino flux of the popu-
for three different generic power-law spectra using the γ- lation that originates from the source resolved in γ-rays
weighting scheme. The ratio of the differential sensitivity is equal to the fraction of the total γ-ray flux that orig-
curve with a convex energy spectrum generally forms a inates from the resolved sources. It has been estimated
function with a single maximum that falls off towards the that 70% of the total diffuse γ-ray flux from blazars
sides. The central interval enclosing 90% of the area un- between 100 MeV and 100GeV has been resolved in
der the ratio function is used to define the energy range 1FGL blazars at high galactic latitudes |b| > 15◦ (Ajello
that contributes 90% to the sensitivity for a given en- et al. 2015), and similarly for 1FGL FSRQs (Ajello et al.
ergy spectrum. This area-sensitivity relation has been 2012) in particular. The 2LAC catalogue used in this
checked empirically. The methodology is an extension to work contains blazars at galactic latitudes |b| > 10◦ .
a previous method, e.g. used in Aartsen et al. (2014b), At the extra galactic latitudes (10◦ < |b| < 15◦ ) the
which only uses the 90 % central interval of signal events detection efficiency might be worse. However, even if
and thereby neglects the background rate. it unrealistically sharply drops to zero, one can esti-
mate that the total resolved fraction does only shrink
by an amount that is proportional to the ratio of the
10◦ < |b| < 15◦ sky fraction (≈ 8% of the total sky)
with respect to the rest (≈ 74% of the total sky ), i.e.
Γ = −2.7 Γ = −2.0 Γ = −1.5

−5
10
[GeVs−1 cm−2 sr]

sensitivity curve (black)
10−6 and upper limits (colored)
10−7
10−8
νµ
dΦ
E 2 dE 10−9
10−10
0.8
0.7 ratio limit to sensitivity
0.6
0.5 90 % central interval
a. u.
0.4
0.3
0.2
0.1
0.0
[GeVs−1 cm−2 sr]
10−7
boundary-corrected upper limits
10−8
10−9
νµ
dΦ
E 2 dE
10−10
2 3 4 5 6 7 8 9
log10 (Eν /GeV)
Figure 8. Determination of the energy range that contributes 90% to the total sensitivity of IceCube for a neutrino flux of a given
spectrum. The construction is shown for the total 2LAC-blazar population using the γ-energy flux weighting scheme for three power-law
spectra with spectral indices −1.5, −2.0 and −2.7.
to 0.08·0+0.74·70%
0.82 ≈ 63%. Since this estimate is conser- SUPPLEMENTARY FIGURES
vative, is still within the error on the quoted 70% value,
and since the 2LAC sample is based on the more sen-
sitive 2FGL catalogue, we conclude that the 70% com-
pleteness is a reasonable estimate to choose for the 2LAC
sources. Accordingly, in a first step one can use a scaling
factor for the neutrino flux upper limits of 1.4 ≈ 1/0.7
to account for the contributions of the blazars that are
not in our sample. In the scenario that high-energy γ-
rays from blazars are “dissipated”, i.e. isotropized in
EBL-induced γγ-cascades due to intergalactic magnetic
fields (Aharonian et al. 1994), the fraction of neutrinos
emitted from the sources not resolved in γ-rays could
be higher. A simple estimation based on numbers from
Ajello et al. (2015) shows, however, that even then the
scaling factor must be less than ≈ 2.8. The total EGB
has an intensity of 11.3 × 10−6 photons/cm2 s sr, of
which 4.1 × 10−6 photons/cm2 s sr originates from re-
solved blazars and the other 7.2 × 10−6 photons/cm2 s sr
from the IGRB (isotropic γ-ray background). If we as-
sume that the entire contribution to the IGRB stems
from EBL induced γγ-cascades from γ-rays emitted by
blazars, i.e. unresolvable isotropized γ-rays, the result-
ing ratio between the total emission from blazars and the
emission from resolved blazars would be 11.3/4.1 ≈ 2.8.
Since the intensity of the IGRB can be well explained by
contributions from unresolved - but in principle resolv-
able - blazars, from star-forming galaxies, and from radio
galaxies (Ajello et al. 2015), we deem this maximum cas-
cade emission scenario unlikely and use the factor of 1.43
throughout this work.
rel. weight / bin [%]

6
5 10
4 8
3 6
2 4
1 2
0 0
−1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0
sin(dec) sin(dec)
(a) All 2LAC Blazars (b) FSRQs


6 10
5 8
4 6
3
2 4
1 2
0 0
−1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0
sin(dec) sin(dec)
(c) ISP/HSPs (d) LSPs

15
10
0
−1.0 −0.5 0.0 0.5 1.0
sin(dec)
(e) LSP-BL Lacs

Figure 9. Relative contribution to the total sum of all source weights for a given declination bin (in percent). The γ-weighting scheme is
shown in green and the equal weighting scheme is shown in black. The binning is chosen such that at most 5 sources fall into a bin for the
largest sample.
REFERENCES Ajello, M., Gasparrini, D., Sánchez-Conde, M., et al. 2015, The
Astrophysical Journal, 800, L27
Ajello, M., Romani, R., Gasparrini, D., et al. 2014, The
Aartsen, M. G., et al. 2013a, Science, 342, 1242856 Astrophysical Journal, 780, 73
—. 2013b, The Astrophysical Journal, 779, 132 Ajello, M., Shaw, M., Romani, R., et al. 2012, The Astrophysical
—. 2014a, Physical Review, D89, 062007 Journal, 751, 108
—. 2014b, The Astrophysical Journal, 796, 109 Atoyan, A., & Dermer, C. 2001, Physical Review Letters, 87,
—. 2015a, Physical Review, D91, 022001 221102
—. 2015b, The Astrophysical Journal, 809, 98 Barlow, R. 1990, Nuclear Instruments and Methods in Physics
—. 2015c, arXiv:1510.05222 Research, A297, 496
—. 2016a, accepted by Physical Review Letters, arXiv:1607.05886 Becker, J., Biermann, P., & Rhode, W. 2005, Astroparticle
—. 2016b, accepted by Astrophysical Journal, arXiv:1607.08006 Physics, 23, 355
Abbasi, R., et al. 2009, Nuclear Instruments and Methods in Böttcher, M., & Dermer, C. 2002, The Astrophysical Journal,
Physics Research, A601, 294 564, 86
—. 2011, Physical Review, D84, 082001 Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, The
Abdo, A., Ackermann, M., Ajello, M., et al. 2010a, The Astrophysical Journal, 768, 54
Astrophysical Journal, 716, 30 Bugaev, E., Montaruli, T., Shlepin, Y., & Sokalski, I. 2004,
—. 2010b, The Astrophysical Journal Supplement Series, 188, 405 Astroparticle Physics, 21, 491
—. 2010c, The Astrophysical Journal, 720, 435 Cavaliere, A., & D’Elia, V. 2002, The Astrophysical Journal, 571,
Ackermann, M., Ajello, M., Allafort, A., et al. 2011, The 226
Astrophysical Journal, 743, 171 Domı́nguez, A., Finke, J., Prada, F., et al. 2013, The
—. 2013, The Astrophysical Journal Supplement Series, 209, 34 Astrophysical Journal, 770, 77
Ackermann, M., Ajello, M., Atwood, W. B., et al. 2015, The Falcke, H., & Biermann, P. 1995, Astronomy and Astrophysics,
Astrophysical Journal, 810, 14 293, 665
—. 2016, The Astrophysical Journal Supplement Series, 222, 5 Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini,
Aharonian, F. A., Coppi, P. S., & Volk, H. J. 1994, The G. 1998, Monthly Notices of the Royal Astronomical Society,
Astrophysical Journal, 423, L5 299, 433
Ahlers, M., Bai, Y., Barger, V., & Lu, R. 2016, Physical Review,
93, 013009
10−5 10−5
[GeVs−1 cm−2 sr−1 ] All 2LAC blazars FSRQs

10−6 10−6
10−7 10−7
10−8 10−8
νµ
νµ
10−9 10−9
dΦ
dΦ
E 2 dE
E 2 dE
equal weighting (median SCD outcome) equal weighting (median SCD outcome)
γ-weighting γ-weighting
10−10 2 10−10 2
10 103 104 105 106 107 108 109 10 103 104 105 106 107 108 109
100 100
local p-value
local p-value
0σ 0σ
significance
significance
10−1 1σ 10−1 1σ
2σ 2σ
10−2 10−2
10−3 3σ 10−3 3σ
(a) All 2LAC Blazars (b) FSRQs

−5 −5
10 10
ISP/HSPs LSPs

10−6 10−6
10−7 10−7
10−8 10−8
νµ
νµ
10−9 10−9
dΦ
dΦ
E 2 dE
E 2 dE
equal weighting (median SCD outcome) equal weighting (median SCD outcome)
γ-weighting γ-weighting
10−10 2 10−10 2
10 103 104 105 106 107 108 109 10 103 104 105 106 107 108 109
100 100
local p-value
local p-value
0σ 0σ
significance
significance
10−1 1σ 10−1 1σ
2σ 2σ
10−2 10−2
10−3 3σ 10−3 3σ
(c) ISP/HSPs (d) LSPs

−5
10
LSP-BL Lacs
90% C.L. Upper Limit

10−6
10−7
10−8
νµ
10−9
dΦ
E 2 dE
equal weighting (median SCD outcome)

γ-weighting
10−10 2
10 103 104 105 106 107 108 109
100
local p-value
0σ
significance
10−1 1σ
2σ
10−2
10−3 3σ
(e) LSP-BL Lacs

Figure 10. 90% C.L. differential upper limits on the diffuse (νµ + ν µ )-flux and corresponding local p-values for the different blazar
populations. The equal-weighting upper limits represent the median SCD sampling outcome.
Georgakakis, A., Nandra, K., Laird, E. S., Aird, J., & Trichas, M. Padovani, P., Petropoulou, M., Giommi, P., & Resconi, E. 2015,
2008, Monthly Notices of the Royal Astronomical Society, 388, Monthly Notices of the Royal Astronomical Society, 452, 1877
1205 Padovani, P., Resconi, E., Giommi, P., Arsioli, B., & Chang,
Giommi, P., Padovani, P., Polenta, G., et al. 2012, Monthly Y. L. 2016, Monthly Notices of the Royal Astronomical Society,
Notices of the Royal Astronomical Society, 420, 2899 457, 3582
Halzen, F., & Zas, E. 1997, The Astrophysical Journal, 488, 669 Planck Collaboration. 2015, accepted by Astronomy &
Hill, G., & Rawlins, K. 2003, Astroparticle Physics, 19, 393 Astrophysics, arXiv:1502.01589
Hopkins, A. M., Afonso, J., Chan, B., et al. 2003, The Protheroe, R. 1997, ASP Conference Series, 121, 585
Astronomical Journal, 125, 465 Read, A. 2000
Kadler, M., Krauß, F., Mannheim, K., et al. 2016, Nature Physics Schuster, C., Pohl, M., & Schlickeiser, R. 2002, Astronomy and
Mannheim, K. 1995, Astroparticle Physics, 3, 295 Astrophysics, 382, 829
Meyer, E. T., Fossati, G., Georganopoulos, M., & Lister, M. L. Stecker, F. W., Done, C., Salamon, M. H., & Sommers, P. 1991,
2011, The Astrophysical Journal, 740, 98 Physical Review Letters, 66, 2697, [Erratum: Phys. Rev.
Mücke, A., Protheroe, R., Engel, R., Rachen, J., & Stanev, T. Lett.69,2738(1992)]
2003, Astroparticle Physics, 18, 593 Stickel, M., Fried, J. W., Kuehr, H., Padovani, P., & Urry, C. M.
Murase, K., Guetta, D., & Ahlers, M. 2015, Physical Review 1991, The Astrophysical Journal, 374, 431
Letters, 116, 071101 Tavecchio, F., & Ghisellini, G. 2015, Monthly Notices of the
Murase, K., Inoue, Y., & Dermer, C. 2014, Physical Review, D90, Royal Astronomical Society, 451, 1502
023007 Tavecchio, F., Ghisellini, G., & Guetta, D. 2014, The
Neronov, A., & Semikoz, D. 2016, Astroparticle Physics, 75, 6063 Astrophysical Journal, 793, L18
Neunhoffer, T. 2006, Astroparticle Physics, 25, 220 Urry, C., & Padovani, P. 1995, Publications of the Astronomical
Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012, The Society of the Pacific, 107, 803
Astrophysical Journal Supplement Series, 199, 31 Wang, B., & Li, Z. 2015, Science China Physics, Mechanics &
Olive, K. 2014, Chinese Physics C, 38, 090001 Astronomy, 59, 1
Padovani, P., & Giommi, P. 1995, The Astrophysical Journal,
444, 567

Draft Version November 15, 2016: Preprint Typeset Using L TEX Style Emulateapj v. 01/23/15

Uploaded by

Copyright:

Available Formats

Draft Version November 15, 2016: Preprint Typeset Using L TEX Style Emulateapj v. 01/23/15

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Draft Version November 15, 2016: Preprint Typeset Using L TEX Style Emulateapj v. 01/23/15

Uploaded by

Copyright:

Available Formats

Draft version November 15, 2016

Preprint typeset using LATEX style emulateapj v. 01/23/15

THE CONTRIBUTION OF FERMI-2LAC BLAZARS TO THE DIFFUSE TEV-PEV NEUTRINO FLUX

J. P. A. M. de André35 , C. De Clercq29 , E. del Pino Rosendo11 , H. Dembinski36 , S. De Ridder37 , P. Desiati6 ,

1. INTRODUCTION with increasing accuracy (Aartsen et al. 2015a). The

20742, USA 53115 Bonn, Germany

mund, Germany Philadelphia, PA 19104, USA

746, Korea 54022, USA

Belgium Tuscaloosa, AL 35487, USA

versité de Genève, CH-1211 Genève, Switzerland ford OX1 3NP, UK

Canada, M5S 1A7 Bunkyo, Tokyo 113-0032, Japan

Münster, D-48149 Münster, Germany USA

FSRQs (310 objects): The class of FSRQs show Table 1

All 2LAC Blazars LSP-BL Lacs

Fermi-LAT in the energy range between E > 100 MeV (2011).

the weights in dependence of declination can be found in 4.4. Simulations

for a particular model, and is also used for visualization

[GeVs−1 cm−2 sr−1 ]

19%-27% contribution of the total 2LAC blazar sample

10−9 All 2LAC blazars (equal weighting)

be dominated by the soft regime, as the ”ratio function” model.

Type Model MRF

10−5 10−4 Models 90% C.L.

−6 Mücke et al. LSP (optimistic)

(a) generic blazars (b) BL Lacs

[GeVs−1 cm−2 sr−1 ]

10−6 Murase BLR-Zone, Γp = 2.0 10−6 Murase BLR-Zone, Γp = 2.3

(c) FSRQs - 1 (d) FSRQs - 2

90 % central interval from dN/dS sample realizations

u.l. this analysis - particular realizations

Γ = −2.7 Γ = −2.0 Γ = −1.5

[GeVs−1 cm−2 sr]

rel. weight / bin [%]

(a) All 2LAC Blazars (b) FSRQs

rel. weight / bin [%]

(c) ISP/HSPs (d) LSPs

(e) LSP-BL Lacs

[GeVs−1 cm−2 sr−1 ]

(a) All 2LAC Blazars (b) FSRQs

[GeVs−1 cm−2 sr−1 ]

(c) ISP/HSPs (d) LSPs

90% C.L. Upper Limit

equal weighting (median SCD outcome)

(e) LSP-BL Lacs

You might also like