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

The initial discovery of astrophysical neutrino flux around PeV energies (Aartsen et al. 2013a) a few years ago marked the beginning of high-energy neutrino astronomy. Since then, the properties of the flux have been measured with increasing accuracy (Aartsen et al. 2015a). The most recent results indicate a soft spectrum with a spectral index of −2.5 ± 0.1 between around 10 TeV and 2 PeV with no significant deviation from an equal flavor composition on Earth (Aartsen et al. 2015b). The neutrino signal has been found to be compatible with an isotropic distribution in the sky. This apparent isotropy suggests that a significant fraction of the observed neutrinos are of extragalactic origin, a result which is also supported by Ahlers et al. (2016). However, there are also indications for a 3-σ anisotropy (Neronov & Semikoz 2016) if low-energy events ($\lt 100$ TeV) are omitted. Further data are required to settle this issue.

Potential extragalactic sources are active galactic nuclei (AGNs): both radio-quiet (Stecker et al. 1991) and radio-loud (Mannheim 1995) objects have been considered as the origin of neutrino production for many years. Blazars, a subset of radio-loud AGNs with relativistic jets pointing towards Earth (Urry & Padovani 1995), are investigated in this paper. They are commonly classified according to the properties of the spectral energy distribution (SED) of their electromagnetic emission. The blazar SED features two distinctive peaks: a low-energy peak between infrared and X-ray energies, attributed to the synchrotron emission of energetic electrons, and a high-energy peak at γ-ray energies, which can be explained by several and possibly competing interaction and radiation processes of high-energy electrons and high-energy nuclei (Böttcher et al. 2013). Several works suggest that blazar SEDs follow a sequence (Fossati et al. 1998; Böttcher & Dermer 2002; Cavaliere & D'Elia 2002; Meyer et al. 2011) in which the peak energy of the synchrotron emission spectrum decreases with increasing blazar luminosity. Accordingly, blazars can be classified into low synchrotron peak (LSP), intermediate synchrotron peak (ISP), and high synchrotron peak (HSP) objects⁵⁶ , a classification scheme introduced in Abdo et al. (2010a) which we use throughout this work. A second classifier is based on the prominence of emission lines in the SED over the non-thermal continuum emission of the jet. Flat-spectrum radio quasars (FSRQs) show Doppler-broadened optical emission lines (Stickel et al. 1991), while in so-called BL Lac objects, emission lines are hidden under a strong continuum emission.

Many calculations of high-energy neutrino emission from the jets of blazars can be found in the literature. Neutrinos could be produced via charged pion decay in interactions of high-energy protons with gas (pp-interactions) in the jets (Schuster et al. 2002) or in interactions of protons with internal (Mannheim 1995) or external (Atoyan & Dermer 2001) photon fields (pγ-interactions). Early models for the neutrino emission from blazars made no explicit distinction based on the blazar class. Some of these classes have already been explicitly excluded at 90% C.L. by past diffuse neutrino flux measurements (Abbasi et al. 2011; Aartsen et al. 2014a)—for example, the combined pp+pγ predictions in Mannheim (1995). More recent publications, on the other hand, differentiate between specific classes of blazars and are largely not yet constrained by experiment. The neutrino production of BL Lac objects is modeled, e.g., in Mücke et al. (2003), Tavecchio et al. (2014), and Padovani et al. (2015), while the neutrino production of FSRQs is calculated, e.g., in Becker et al. (2005) and Murase et al. (2014). The models by Tavecchio et al. (2014) and Padovani et al. (2015) were in particular constructed to explain parts or all of the astrophysical neutrino flux. With the analysis presented here, we are able to test large parts of the parameter space of many of these models for the first time. We do not consider theoretical calculations from the literature for individual objects, since these are not directly comparable to our results.

The neutrinos predicted by most models are produced in charged pion decays, which come with an associated flux of γ-rays from neutral pion decays. Even if the hadronic fraction is subdominant, one could on average expect a higher neutrino luminosity for a higher observed γ-luminosity (Murase et al. 2014). On a source-by-source basis, however, variations in the exact $\nu /\gamma$ correlation are likely. One strategy to cope with this uncertainty, which we follow in this paper, is to analyze large samples of objects and thereby to investigate average properties. We use the Fermi-LAT 2LAC catalog⁵⁷ (Ackermann et al. 2011) to define search positions for our analysis (see Section 2). The blazars in the 2LAC catalog constitute the majority ($\approx 70 \%$) of the total γ-ray flux emitted from all GeV blazars in the observable universe between 100 MeV and 100 GeV (see Appendix C). The 2LAC contains more than twice the number of blazars compared to other Fermi catalogs starting at higher energies, such as 1FHL (Ackermann et al. 2013) or 2FHL (Ackermann et al. 2016). The goal is to look for a cumulative neutrino flux excess from all 862 2LAC blazars or from specifically selected subpopulations using muon-track data with an angular resolution of about a degree in an unbinned maximum-likelihood stacking approach. We use two different "weighting schemes" (see Section 4.2) to define the probability density functions (PDFs) for the neutrino signal, expressing different assumptions about the relative neutrino flux for each source. Each weighting scheme represents its own test of the data.

The analysis most drastically differs from previous point source searches in two points:

• 1.  
  The blazar populations are comprised of nearly 2 orders of magnitude more sources.
• 2.  
  For the first time, we use a model-independent weighting scheme. In this test of the data, we make no assumption about the exact $\nu /\gamma$ correlation except that the neutrino flux originates from the defined blazar positions.

Section 2 defines the five blazar populations considered in this analysis. Section 3 describes the muon-track data set used for this search. Section 4 summarizes the analysis, including the technique of the unbinned stacking search, a description of different weighting schemes, the confidence interval construction, and a discussion on potential biases from non-hadronic contributions to the γ-ray flux. Section 5 presents the analysis results, and Section 6 discusses their implications.

The Fermi-LAT 2LAC catalog (Ackermann et al. 2011) contains 862 GeV-emitting blazars at high galactic latitudes $| b| \gt 10^\circ$ that are not affected by potential source confusion.⁵⁸ The data for this catalog were taken between August 2008 and August 2010. We use the spectroscopic classification into FSRQ and BL Lac objects (Stickel et al. 1991) and the independent classification into LSP, ISP, and HSP objects (Abdo et al. 2010a) to define subpopulations of these 862 objects. We do not impose any other cuts (e.g., on the γ-ray flux) because the exact neutrino flux expectations are unknown, as outlined in Section 1. The motivations for the particular subsamples are described in the following.

• All 2LAC Blazars (862 objects). The evolutionary blazar sequence (Cavaliere & D'Elia 2002; Böttcher & Dermer 2002) suggests that blazars form a continuous spectrum of objects that are connected via cosmological evolution. A recent study by Ajello et al. (2014) supports this hypothesis. Since the corresponding evolution of the neutrino emission is not known, the most unbiased assumption is to group all blazars together. This is especially justified for the analysis using the equal-weighting scheme discussed in Section 4.2.
• FSRQs (310 objects). The class of FSRQs shows strong, broad emission lines that potentially act as intense radiation targets for the photomeson production of neutrinos (Atoyan & Dermer 2001; Murase et al. 2014).
• LSPs (308 objects). The majority of FSRQs are LSP objects. Giommi et al. (2012) argue that LSP-BL Lacs are actually physically similar to FSRQs but have emission lines that are overwhelmed by the strong jet continuum. This sample therefore groups all LSP objects together.
• ISP+HSPs (301 objects). HSP objects differ from LSP objects in terms of luminosity and mainly consist of BL Lacs (Ajello et al. 2014). The peak-frequency boundary between LSP and HSP objects is only defined artificially, with ISP objects filling the gap. In order to have a larger sample of objects, the HSP objects are grouped together with the ISP objects in one combined sample for this analysis.
• LSP-BL Lacs (68 objects). Objects that share the LSP and BL Lac classification have been specifically considered for neutrino emission in Mücke et al. (2003). Therefore, we test them as a separate sample. They form the smallest subpopulation in this analysis.

The distribution of the sources in the sky for the largest sample (all 2LAC blazars) and smallest sample (LSP-BL Lacs) is shown in Figure 1. A modest LAT-exposure deficit and lower sky coverage by optical surveys in the southern sky lead to a slight deficit of objects in the southern hemisphere (Ackermann et al. 2011). The effect is most prominent for the BL Lac-dominated samples. However, blazars without optical association are also included in the 2LAC catalog and partly make up for this asymmetry in the total sample. For simplicity, we assume a quasi-isotropic source distribution for all populations (excluding the source-free region around the galactic plane) for the calculation of quasi-diffuse fluxes. This assumption also seems reasonable given the weight distribution of sources (equal weighting) in Figures 9(a)–(e) and Appendix D. Figure 2 shows the overlap between the samples. The LSP-BL Lac, FSRQ, and ISP+HSP samples are nearly independent of one another, with a small overlap of 3 sources between the FSRQ and ISP+HSP samples. The largest overlap exists between the FSRQ and LSP samples, which share around 60% of their sources. The all-blazar sample contains 167 sources that are not shared with any subsample. These are sources that are either unclassified or only classified as BL Lac objects with no corresponding synchrotron peak classification.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

IceCube is a neutrino telescope located at the geographic South Pole. It consists of about one ${\mathrm{km}}^{3}$ of Antarctic ice that is instrumented with 5160 optical photosensors which are connected via cables ("strings") with the data aqcuisition system at the surface. The photosensors detect Cherenkov light emitted by charged particles that are produced in the interactions of neutrinos with nuclei and electrons in the glacial ice. The geometry and sensitivity of the photosensors lead to an effective energy threshold of about 100 GeV for neutrinos. A more detailed description of the detector and the data acquisition can be found in Abbasi et al. (2009).

Two main signatures can be distinguished for the recorded events: "track-like" and "shower-like." Only track-like events are of interest for the analysis here. They are the characteristic signature of muons produced in the charged-current interactions of muon neutrinos.⁵⁹

IceCube was constructed between 2006 and 2011 with a final configuration of 86 strings. We use data from the 59-string (IC-59), 79-string (IC-79), and 86-string (IC-86) configurations of the IceCube detector recorded between May 2009 and April 2012. In contrast to previous publications, we do not include data from the 40-string configuration here since the ice model description in the IC-40 Monte Carlo data sets is substantially different and the sensitivity gain would be marginal. The track event selection for the three years of data is similar to the one described in Aartsen et al. (2013b, 2014b). The angular resolution of the majority of events in the track sample is better than 1° for events with reconstructed muon energies above 10 TeV (Aartsen et al. 2014b). The angular reconstruction uncertainty is calculated following the prescription given in Neunhoffer (2006). We apply one additional minor selection criterion for the estimated angular uncertainty of the reconstructed tracks (${\sigma }_{\mathrm{est}.}\leqslant 5^\circ$) for computational reasons. The removed events do not have any measurable effect on the sensitivity. Event numbers for the individual data sets are summarized in Table 1.

The data set is dominated by bundles of atmospheric muons produced in cosmic-ray air shower interactions for tracks coming from the southern hemisphere ($\theta \lt 90^\circ$). Tracks from the northern hemisphere ($\theta \geqslant 90^\circ$) originate mostly from atmospheric neutrino interactions that produce muons. In order to reduce the overwhelming background of direct atmospheric muons to an acceptable level, it is necessary to impose a high-energy cut for events from the southern hemisphere. The cut raises the effective neutrino energy threshold to approximately 100 TeV (Aartsen et al. 2014b), reducing the sensitivity to neutrino sources in this region by at least 1 order of magnitude for spectra softer than ${{\rm{E}}}^{-2}$. Only for harder spectra does the southern sky have a significant contribution to the overall sensitivity. The northern sky does not require such an energy cut, as upgoing tracks can only originate from neutrino interactions, which have a much lower incidence rate. However, at very high energies (again around $100\,\mathrm{TeV}$), the Earth absorbs a substantial fraction of neutrinos, reducing the expected astrophysical signal as well. Charged-current ${\nu }_{\mu }$ interactions can happen far outside the instrumented volume and still be detected, as high-energy muons may travel several kilometers through the glacial ice before entering the detector. This effect increases the effective detection area for certain arrival directions, mostly around the horizon.

The most sensitive region is therefore around the celestial equator, which does not require a high-energy cut, provides ample target material surrounding the detector, i.e., a large effective area, and does not suffer from absorption of neutrinos above 100 TeV. However, these zenith-dependent sensitivity changes are mostly important for the interpretation of the results (see, e.g., Section 5.3). The likelihood approach takes these differences into account with the "acceptance" term in Equation (6), Section 4.1, and a separation into several zenith-dependent analyses is not necessary. For more details on the properties of the data sets and the zenith-dependent sensitivity behavior, we refer the reader to Aartsen et al. (2013b) and Aartsen et al. (2014b).

4.1. The Likelihood Function for Unbinned ML Stacking

Analysis is performed via an extended unbinned maximum-likelihood fit (Barlow 1990). The likelihood function consists of two PDFs, one PDF $B(\overline{x})$ for a background hypothesis and one PDF $S(\overline{x})$ for a signal hypothesis. 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

$$\begin{eqnarray}\mathrm{ln}(L)\{{n}_{s},{{\rm{\Gamma }}}_{\mathrm{SI}}\} & = & \displaystyle \sum _{i=1}^{N}\mathrm{ln}\left(\displaystyle \frac{{n}_{s}}{N}\cdot S({\delta }_{i},R.A{.}_{i},{\sigma }_{i},{\varepsilon }_{i};{{\rm{\Gamma }}}_{\mathrm{SI}})\right.\\ & & \left.+\left(1-\displaystyle \frac{{n}_{s}}{N}\right)\cdot B(\sin ({\delta }_{i}),{\varepsilon }_{i})\right),\end{eqnarray} \tag{ 1 }$$

where i indexes individual neutrino events. The likelihood function depends on two free parameters: the normalization factor n_s and spectral index ${{\rm{\Gamma }}}_{\mathrm{SI}}$ of the total blazar signal. For computational reasons we assume that each source of a given population shares the same spectral index. The background evaluation for each event depends on the reconstructed declination ${\delta }_{i}$ and the reconstructed muon energy ${\varepsilon }_{i}$. The signal part additionally depends on the reconstructed right ascension ${\rm{R}}.{\rm{A}}{.}_{i}$, the angular error estimator ${\sigma }_{i}$, and the power-law spectral index ${{\rm{\Gamma }}}_{\mathrm{SI}}$.

The background PDF is constructed from binning the recorded data in the reconstructed declination and energy. It is evaluated as

$$\begin{eqnarray}&&B(\sin ({\delta }_{i}),{\varepsilon }_{i})=\displaystyle \frac{1}{2\pi }\cdot f(\sin ({\delta }_{i}),{\varepsilon }_{i}),\end{eqnarray} \tag{ 2 }$$

where $\tfrac{1}{2\pi }$ arises from integration over the right ascension and f is the normalized joint probability distribution of the events in declination $\sin (\delta )$ and energy ε.

The signal PDF that describes a given blazar population is a superposition of the individual PDFs for each source,

$$\begin{eqnarray}&&S({\delta }_{i},{\rm{R}}.{\rm{A}}{.}_{i},{\sigma }_{i},{\varepsilon }_{i};{{\rm{\Gamma }}}_{\mathrm{SI}})\\ &&\quad =\displaystyle \frac{{\displaystyle \sum }_{j=1}^{{N}_{\mathrm{src}}}{w}_{j}\cdot {S}_{j}({\delta }_{i},{\rm{R}}.{\rm{A}}{.}_{i},{\sigma }_{i},{\varepsilon }_{i};{{\rm{\Gamma }}}_{\mathrm{SI}})}{{\displaystyle \sum }_{j=1}^{{N}_{\mathrm{src}}}{w}_{j}},\end{eqnarray} \tag{ 3 }$$

where w_j is a weight determining the relative normalization of the PDF S_j for source j. This weight therefore accounts for the relative contribution of source j to the combined signal. In general, different choices of w_j are possible. The two choices used in this work are discussed in Section 4.2. Each term S_j in Equation (3) is evaluated as

$$\begin{eqnarray}&&{S}_{j}({\delta }_{i},{\rm{R}}.{\rm{A}}{.}_{i},{\sigma }_{i},{\varepsilon }_{i};{{\rm{\Gamma }}}_{\mathrm{SI}})\\ &&\quad =\displaystyle \frac{1}{2\pi {{\sigma }_{i}}^{2}}\cdot \exp \left(\displaystyle \frac{1}{2}\cdot {\left(\displaystyle \frac{{{\rm{\Psi }}}_{{ij}}[{\delta }_{i},{\rm{R}}.{\rm{A}}{.}_{i}]}{{\sigma }_{i}}\right)}^{2}\right)\cdot {g}_{j}({\varepsilon }_{i};{{\rm{\Gamma }}}_{\mathrm{SI}}),\end{eqnarray} \tag{ 4 }$$

where the spatial term is expressed as a 2D symmetric normal distribution and g_j is the normalized PDF for the reconstructed muon energy for source j. The term ${{\rm{\Psi }}}_{{ij}}$ is the angular separation between event i and source j.

4.2. Weighting Schemes

The term w_j in Equation (3) parametrizes the relative contribution of source j to the combined signal. It corresponds to the expected number of events for source j and can be expressed as

$$\begin{eqnarray}&&{w}_{j}={\int }_{{E}_{\nu ,\min }}^{{E}_{\nu ,\max }}{{\rm{\Phi }}}_{0,j}\cdot {h}_{j}({E}_{\nu })\cdot {A}_{\mathrm{eff}}({\theta }_{j},{E}_{\nu })\ {{dE}}_{\nu },\end{eqnarray} \tag{ 5 }$$

where ${A}_{\mathrm{eff}}({\theta }_{j},{E}_{\nu })$ is the effective area for incoming muon neutrinos from a given source direction at a given energy, ${h}_{j}({E}_{\nu })$ denotes the normalized neutrino energy spectrum for source j, and ${{\rm{\Phi }}}_{0,j}$ is its overall flux normalization. The integration bounds ${E}_{\nu ,\min }$ and ${E}_{\nu ,\max }$ are set to ${10}^{2}$ and ${10}^{9}\ \mathrm{GeV}$, respectively, except for the differential analysis (see Section 4.3), in which they are defined for the given energy band.

Under the assumption that all sources share the same spectral power-law shape, w_j is further simplified via

$$\begin{eqnarray}{w}_{j} & = & [{{\rm{\Phi }}}_{0,j}]\cdot \left[{\displaystyle \int }_{{E}_{\nu ,\min }}^{{E}_{\nu ,\max }}h({E}_{\nu };{{\rm{\Gamma }}}_{\mathrm{SI}})\cdot {A}_{\mathrm{eff}}({\theta }_{j},{E}_{\nu })\ {{dE}}_{\nu }\right]\\ & = & [C\cdot {w}_{j,\mathrm{model}}]\cdot [{w}_{j,\mathrm{acc}.}({\theta }_{j},{{\rm{\Gamma }}}_{\mathrm{SI}})]\end{eqnarray} \tag{ 6 }$$

and splits into a "model" term ${w}_{j,\mathrm{model}}$—which is proportional to the expected relative neutrino flux of source j—and an "acceptance" term, which is fixed by the position of the source and the global energy spectrum. The term ${w}_{j,\mathrm{model}}$ is not known, and its choice defines the "weighting scheme" for the stacking analysis. The following two separate weighting schemes are used for the signal PDF in the likelihood analysis, leading to two different sets of tests.

4.2.1. γ-weighting

For this weighting scheme, we first have to assume that the γ-ray flux can be modeled to be quasi-steady between 2008 and 2010, the time period which forms the basis for the 2LAC catalog. This makes it possible to extrapolate the flux expectation of each source to other time periods, e.g., into the non-overlapping part of the data-taking period of the IceCube data for this analysis (2009–2012). Each model weight, i.e., the relative neutrino flux expected to arrive from a given source, is then given by the source's γ-ray energy flux observed by Fermi-LAT in the energy range between $E\gt 100\ \mathrm{MeV}$ and $E\gt 100\ \mathrm{GeV}$:

$$\begin{eqnarray}&&{w}_{j,\mathrm{model}}={\int }_{100\mathrm{MeV}}^{100\mathrm{GeV}}{E}_{\gamma }\displaystyle \frac{d{\phi }_{\gamma ,j}}{{{dE}}_{\gamma }}\ {{dE}}_{\gamma }.\end{eqnarray} \tag{ 7 }$$

This is motivated by the fact that a similar amount of energy is channeled into the neutrino and γ-ray emission if pion decay from pp or $p\gamma$ interactions dominates the high-energy interaction. While the source environment is transparent to high-energy neutrinos, it might not be for γ-rays. The reprocessing of γ-rays due to $\gamma \gamma$ interactions might then shift the energies of the photons to GeV and sub-GeV energies before they can leave the sources, which would make them detectable by the Fermi-LAT. This might even be expected in $p\gamma$ scenarios (Murase et al. 2015). Since a large fraction of blazars are located at high redshifts $z\geqslant 1$,⁶⁰ this reprocessing will also take place during the propagation of photons in extragalactic background light (EBL), shifting γ-ray energies below a few hundred GeV for such sources (Domínguez et al. 2013). This again potentially places them in the energy range of the Fermi-LAT 2LAC catalog. Even in the case where synchrotron contributions (e.g., muon or pion synchrotron radiation) dominate pion decay in the $\mathrm{MeV}$–$\mathrm{GeV}$ range, which has been considered for BL Lac objects in particular (Mücke et al. 2003), one would expect the overall γ-ray emission to be proportional to the neutrino emission. This is also the case in models where inverse Compton processes dominate the high-energy γ-ray emission (Murase et al. 2014).

The preceding arguments in favor of a γ-weighting scheme assume that all sources show equal proportionality. On a source-by-source basis, however, the proportionality factor can vary, as already mentioned in Section 1.

One contributing factor is the fact that Fermi probes different sections of the blazar γ-ray peak for each source relative to the peak position. For simplicity, we do not perform a spectral source-by-source fit in this paper, leaving this aspect for potential future work. This is also mostly an issue for the "all 2LAC blazar" sample, since the other subclassifications described in Section 2 depend on the peak position and this effect is largely mitigated. There are additional reasons for source-by-source fluctuations in the $\gamma /\nu$ correlation due to EBL reprocessing. First, EBL absorption might not be sufficient for close-by sources, such that emerging high-energy γ-rays are not reprocessed into the energy range of the 2LAC catalog, which ends at 100 GeV. Second, EBL reprocessing differs between sources depending on the line-of-sight magnetic fields, which deflect charged particle pairs produced in EBL cascades differently (Aharonian et al. 1994). Third, strong in-source $\gamma \gamma$ reprocessing could lead to γ-rays at even lower energies than 100 MeV (Murase et al. 2015), which would be below the 2LAC energy range.

All results presented in Section 5 based on the γ-weighting scheme assume that the potential source-to-source fluctuations in the $\gamma \mbox{--}\nu$ correlation described here average out for large source populations and can be neglected. More information on the distribution of weights according to declination can be found in Figures 9(a)–(e) and Appendix D.

4.2.2. Equal Weighting

The γ-weighting scheme is optimal under the assumption that the neutrino flux follows the measured γ-energy flux exactly. Given the uncertainties discussed in Section 4.2.1, we also use another weighting scheme,

$$\begin{eqnarray}&&{w}_{j,\mathrm{model}}=1,\end{eqnarray} \tag{ 8 }$$

which we expect to be more sensitive eventually if the actual $\gamma \mbox{--}\nu$ correlation varies strongly from source to source. It provides a complementary and model-independent test which is maximally agnostic to the degree of correlation between γ-ray and neutrino luminosities.

We do not assume a specific neutrino emission—equal emission in particular—in a given source when calculating the flux upper limits for the equal-weighting scheme. We only assume, to some approximation, that the differential source count distributions (SCDs) of γ-rays and neutrinos have comparable shapes. The differential SCD, dN/dS, describes how the energy-integrated flux S is distributed over all sources, and is a crucial property of any cosmological source population. Section 4.4 provides more information on the technical aspects of neutrino flux injection in the equal-weighting test. Appendix A then discusses why the methodology is robust against variations in the actual shape of the dN/dS distribution for the neutrino flux in the IceCube energy range and why the final result is valid even if the neutrino SCD is different from the γ-ray SCD.

4.3. Statistical Tests

We perform statistical tests for each population of blazars. The log-likelihood difference λ defines our test statistic (TS), given by

$$\begin{eqnarray}\lambda & = & -\,2\cdot \mathrm{log}(L)\{{n}_{s}=0\}\\ & & +2\cdot \mathrm{log}(L)\{{n}_{s}={n}_{s,\max },{{\rm{\Gamma }}}_{\mathrm{SI}}={{\rm{\Gamma }}}_{\mathrm{SI},\max }\},\end{eqnarray} \tag{ 9 }$$

where ${n}_{s,\max }$ and ${{\rm{\Gamma }}}_{\mathrm{SI},\max }$ are the number of signal events and the signal spectral index that maximize the TS. We simulate an ensemble of background-only skymaps where the TS distribution is compared with the TS value obtained from the data. The p-value is then defined as the fraction of skymaps in the background ensemble that has a larger TS value than the one observed. Ensembles of skymaps with different injected signal strengths are then used to calculate the resulting confidence interval. See Section 4.4 for details on the skymap simulations.

In total we perform two distinct types of tests for which p-values are calculated. The first ("integral") assumes a power-law spectrum for the blazar emission over the full energy range observable with IceCube (unless stated otherwise). The second ("differential") assumes a neutrino signal that is confined to a small energy range (half a decade in energy), and has a power-law spectrum with a spectral index of −2 within this range. We perform the differential test for 14 energy ranges between 100 GeV and 1 EeV.

4.4. Simulations

5.1. Observed p-Values

Table 2 summarizes the p-values for the "integral" test (see Section 4.3). Nine out of the ten tests show overfluctuations but no significant excess. We find the strongest overfluctuation, a 6% p-value, using the equal-weighting scheme for all 2LAC blazars. We omit a trial-factor correction because the populations have a large overlap and the result is not significant.

Table 2.  p-Values and the Corresponding Significance in Units of Standard Normal Deviations in the Power-Law Test

Population	p-value
	γ-weighting	Equal Weighting
All 2LAC blazars	$36 \%$ $(+0.4\sigma )$	$6 \%$ $(+1.6\sigma )$
FSRQs	$34 \%$ $(+0.4\sigma )$	$34 \%$ $(+0.4\sigma )$
LSPs	$36 \%$ $(+0.4\sigma )$	$28 \%$ $(+0.6\sigma )$
ISP/HSPs	$\gt 50 \%$	$11 \%$ $(+1.2\sigma )$
LSP-BL Lacs	$13 \%$ $(+1.1\sigma )$	$7 \%$ $(+1.5\sigma )$

Note. The table shows the results for both weighting schemes. The values do not include a trial-factor correction.

Download table as:  ASCIITypeset image

Figure 3 shows the p-values from the corresponding "differential" test. The largest excess is visible in the 5–10 TeV energy band with a pre-trial p-value of $4\cdot {10}^{-3}$. This outcome is totally compatible with a fluctuation of the background, since the effect of multiple trials has to be taken into account, reducing the significance of the observation substantially. Accurate calculation of the trial-corrected p-value is again difficult, as neither the five blazar samples nor the 14 tested energy ranges per sample are independent. We again omit it for simplicity.

Zoom In Zoom Out Reset image size

Comparing the differential p-value plot of all 2LAC blazars with that of the other populations (see Figures 10(a)–(e) in Appendix D), one finds that the overfluctuation is caused by the LSP-BL Lac, FSRQ, and ISP/HSP populations, which are nearly independent of one another and show a small excess in the 5 TeV–20 TeV region. In the γ-weighting scheme, the ISP/HSP p-value distribution is nearly flat, which leads to an overfluctuation in the all 2LAC blazar sample that is weaker than that in the equal-weighting scenario.

5.2. Flux Upper Limits

5.3. Generic Upper Limits

Table 3 shows flux upper limits assuming a generic power-law spectrum for the tested blazar populations, calculated for the three different spectral indices: −1.5, −2.0, and −2.7.

Table 3.  90% C.L. Upper Limits on the Diffuse (${\nu }_{\mu }+{\overline{\nu }}_{\mu }$) Flux from the Different Blazar Populations Tested

Spectrum: ${{\rm{\Phi }}}_{0}\cdot {(E/\mathrm{GeV})}^{-1.5}$
Blazar Class	${{{\rm{\Phi }}}_{0}}^{90 \% }[{\mathrm{GeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}]$
	γ-weighting	Equal Weighting
All 2LAC Blazars	$1.6\times {10}^{-12}$	$4.6\ (3.8\mbox{--}5.3)\times {10}^{-12}$
FSRQs	$0.8\times {10}^{-12}$	$2.1\ (1.0\mbox{--}3.1)\times {10}^{-12}$
LSPs	$1.0\times {10}^{-12}$	$1.9\ (1.2\mbox{--}2.6)\times {10}^{-12}$
ISPs/HSPs	$1.8\times {10}^{-12}$	$2.6\ (2.0\mbox{--}3.2)\times {10}^{-12}$
LSP-BL Lacs	$1.1\times {10}^{-12}$	$1.4\ (0.5\mbox{--}2.3)\times {10}^{-12}$
Spectrum: ${{\rm{\Phi }}}_{0}\cdot {(E/\mathrm{GeV})}^{-2.0}$
Blazar Class	${{{\rm{\Phi }}}_{0}}^{90 \% }[{\mathrm{GeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}]$
	γ-weighting	Equal Weighting
All 2LAC Blazars	$1.5\times {10}^{-9}$	$4.7\ (3.9\mbox{--}5.4)\times {10}^{-9}$
FSRQs	$0.9\times {10}^{-9}$	$1.7\ (0.8\mbox{--}2.6)\times {10}^{-9}$
LSPs	$0.9\times {10}^{-9}$	$2.2\ (1.4\mbox{--}3.0)\times {10}^{-9}$
ISPs/HSPs	$1.3\times {10}^{-9}$	$2.5\ (1.9\mbox{--}3.1)\times {10}^{-9}$
LSP-BL Lacs	$1.2\times {10}^{-9}$	$1.5\ (0.5\mbox{--}2.4)\times {10}^{-9}$
Spectrum: ${{\rm{\Phi }}}_{0}\cdot {(E/\mathrm{GeV})}^{-2.7}$
Blazar Class	${{{\rm{\Phi }}}_{0}}^{90 \% }[{\mathrm{GeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}]$
	γ-weighting	Equal Weighting
All 2LAC Blazars	$2.5\times {10}^{-6}$	$8.3\ (7.0\mbox{--}9.7)\times {10}^{-6}$
FSRQs	$1.7\times {10}^{-6}$	$3.3\ (1.6\mbox{--}5.1)\times {10}^{-6}$
LSPs	$1.6\times {10}^{-6}$	$3.8\ (2.4\mbox{--}5.2)\times {10}^{-6}$
ISPs/HSPs	$1.6\times {10}^{-6}$	$4.6\ (3.5\mbox{--}5.6)\times {10}^{-6}$
LSP-BL Lacs	$2.2\times {10}^{-6}$	$2.8\ (1.0\mbox{--}4.6)\times {10}^{-6}$

Note. The table contains results for power-law spectra with spectral indices of −1.5, −2.0, and −2.7. The equal-weighting column shows the median flux upper limit and the 90% central interval of different sample realizations of the Fermi-LAT source count contribution (in parentheses). All values include systematic uncertainties.

Download table as:  ASCIITypeset image

The distribution of the γ-ray energy flux among the sources in each population governs the flux upper limit in the γ-weighting scheme. It is mostly driven by the declination of the strongest sources in the population, due to the strong declination dependence of IceCube's effective area. For FSRQs, the two sources with the largest γ-weights (3C 454.3 at $\mathrm{decl}{.}_{2000}=16^\circ$ and PKS1510-08 at $\mathrm{decl}{.}_{2000}=-9^\circ$) carry around 15% of the total γ-weight of all FSRQs. Their positions close to the equator place them in the most sensitive region for the IceCube detector, and the γ-weighting upper limits for FSRQs are more than a factor of 2 lower than the corresponding equal-weighting limits. For the LSP-BL Lacs, the two strongest sources (PKS 0426-380 at $\mathrm{decl}{.}_{2000}=-38^\circ$ and PKS 0537-441 at $\mathrm{decl}{.}_{2000}=-44^\circ$) carry nearly 30% of the total γ-weight but are located in the southern sky, where IceCube is not very sensitive. The γ-weighting upper limit is therefore comparable to the equal-weighting upper limit. The reader is referred to Appendix D for more information on the weight distribution.

Figure 4 shows the differential upper limit in comparison to the median sensitivity for all 2LAC blazars using the equal-weighting scheme. This population showed the largest overfluctuation. We plot here the upper limit derived from the median SCD sampling outcome, since in general the equal-weighting upper limit depends on the neutrino flux realization of the SCD (see Appendix A). As expected, the differential limit is slightly higher, by a factor of about 2, than the median outcome in the energy range between 5 and 10 TeV, where the largest excess is observed. This is the average behavior for a soft flux with a spectral index of about −3.0⁶⁴ if one assumes a simple power-law fit to explain the data. While such a physical interpretation cannot be made yet, it will be interesting to observe this excess with future IceCube data. For information on the differential upper limits from the other samples, the reader is referred to Appendix D.

Zoom In Zoom Out Reset image size

5.4. The Maximal Contribution to the Diffuse Astrophysical Flux

Astrophysical neutrino flux is observed between 10 TeV and 2 PeV (Aartsen et al. 2015b). Its spectrum has been found to be compatible with a single power law 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 for the signal injected into the simulated skymaps when calculating the upper limit for a direct comparison. Figure 5 shows the flux upper limit for an ${E}^{-2.5}$ power-law spectrum starting at 10 TeV for both weighting schemes in comparison to the most recent global fit of the astrophysical diffuse neutrino flux, assuming an equal composition of flavors arriving on Earth.

Zoom In Zoom Out Reset image size

The equal-weighting upper limit results in a maximal 19%–27% contribution of the total 2LAC blazar sample to the observed best-fit value of the astrophysical neutrino flux, including systematic uncertainties. This limit is independent of the detailed correlation between the γ-ray and neutrino flux from these sources. The only assumption is that the respective neutrino and γ-ray SCDs have similar shapes (see Section 5.2 for details on the signal injection). We use the Fermi-LAT blazar SCD published in Abdo et al. (2010c) as a template for sampling. However, we find that even if the shape of the SCD differs from the shape of this template, the upper limit still holds and is robust. In Appendix A we discuss the effect of different SCD shapes and how combination with existing point source constraints (Aartsen et al. 2015c) leads to a nearly SCD-independent result, since a point source analysis and a stacking search with equal weights effectively trace opposite parts of the available parameter space for the dN/dS distribution.

If we assume proportionality between the γ-ray and neutrino luminosities of the sources, the γ-weighting limit constrains the maximal flux contribution of all 2LAC blazars to 7% of the observed neutrino flux in the full 10 TeV to 2 PeV range. Since the blazars resolved in the 2LAC account for 70% of the total γ-ray emission from all GeV blazars (Ajello et al. 2015), this further implies that at most 10% of the astrophysical neutrino flux stems from all GeV blazars extrapolated to the whole universe, again in the full 10 TeV to 2 PeV range and assuming the γ-weighting is an appropriate weighting assumption. Table 4 summarizes the maximal contributions for all populations, including the γ-weighting result scaled to the respective total population of sources in the observable universe.

Table 4.  Maximal Contributions to the Best-fit Diffuse Flux from Aartsen et al. (2015b) Assuming the Equipartition of Neutrino Flavors

Population	Weighting Scheme
	Equal	γ	γ (Extrapol.)
all 2LAC blazars	19%–27%	$7 \%$	$10 \%$
FSRQs	$5 \% \mbox{--}17 \%$	$5 \%$	$7 \%$
LSPs	$6 \% \mbox{--}15 \%$	$5 \%$	$7 \%$
ISP/HSPs	$9 \% \mbox{--}15 \%$	$5 \%$	$7 \%$
LSP-BL Lacs	$3 \% \mbox{--}13 \%$	$6 \%$	$9 \%$

Note. The equal-weighting case shows this maximal contribution for the 90% central outcomes of potential dn/ds realizations. The last column shows the maximal contribution of the integrated emission from the total parent population in the observable universe exploiting the γ-ray completeness of the 2LAC blazars (see Appendix C).

Download table as:  ASCIITypeset image

It is interesting to compare these numbers directly to the γ-ray sector. Ajello et al. (2015) show that GeV blazars (100 MeV–100 GeV) contribute approximately 50% to the extragalactic gamma-ray background. The resolved 1FGL (Abdo et al. 2010b) blazar component in particular contributes around 35%. This estimate should be rather similar for the 2LAC blazars studied here, which are defined based on the more recent 2FGL catalog (Nolan et al. 2012) (see Appendix C for a discussion). The 2LAC blazar contribution to the astrophysical neutrino flux is therefore at least a factor of 0.75 smaller than the corresponding extragalactic contribution in the γ-regime. The difference of this contribution between the two sectors becomes substantial (7% maximally allowed contribution for neutrinos versus 35% for γ-rays) if one assumes a γ/ν-correlation.

Figure 5 also shows the equal-weighting constraint for a harder neutrino spectrum with a spectral index of −2.2. This harder spectral index is about 3 standard deviations away from the best-fit value derived in Aartsen et al. (2015b) and can be used as an extremal case given the current observations. Comparison of this upper limit with the hard end of the "butterfly" shows that even in this case less than half of the bulk emission can originate in the 2LAC blazars with minimal assumptions about the relative neutrino emission strengths. Due to the low-count status of the data, we omit multi power-law spectrum tests at this point. However, one can estimate the constraints for more complicated models using Figure 8 in Appendix B, which shows the energy range for a given spectrum that contributes the dominant fraction to the sensitivity. The sensitivity for a possible two-component model that has, for example, a soft component at TeV energies and a hard component in the PeV range would be dominated by the soft regime, as the "ratio function" (see Appendix B, Figure 8) by the hard component above a PeV is negligible. In such a scenario, we expect the constraint to be rather similar to our result from the simple power-law test with a spectral index of −2.5.

5.5. Upper Limits on Models for Diffuse Neutrino Emission

For experimental constraints on existing theoretical calculations, we only considered models for diffuse emission from blazar populations, not predictions for specific objects. These include the calculations by Mannheim (1995), Halzen & Zas (1997), and Protheroe (1997) for generic blazars; the calculations by Becker et al. (2005) and Murase et al. (2014) for FSRQs; and the calculations by Mücke et al. (2003), Tavecchio et al. (2014), Tavecchio & Ghisellini (2015), and Padovani et al. (2015) for BL Lacs.

The upper limits in this section are calculated using the γ-weighting scheme and therefore assume a correlation between the neutrino flux and the measured γ-ray energy flux. This allows us to account for the fraction of the neutrino emission that arises from blazars not detected in γ-rays. The fraction of γ-ray emission from resolved 2LAC blazars in general (including BL Lacs) and from FSRQs in particular is about 70% (Ajello et al. 2015, 2012). Therefore, the flux upper limits for the entire population are a factor of $1/0.7\approx 1.43$ weaker than those derived for the quasi-diffuse flux of the 2LAC blazars. See Appendix C for more details on this factor.

Table 5 summarizes the model rejection factors (Hill & Rawlins 2003)⁶⁵ for all considered models. Many of these models can be constrained by this analysis. Figures 6(a)–(d) visualize the flux upper limits in comparison to the neutrino flux predictions.

Zoom In Zoom Out Reset image size

Table 5.  Summary of Constraints and Model Rejection Factors for the Diffuse Neutrino Flux Predictions from Blazar Populations

Type	Model		MRF
Generic blazars	(Mannheim 1995)	(A)	1.30
		(B)	<0.1
	(Halzen & Zas 1997)		<0.1
	(Protheroe 1997)		<0.1
FSRQs	(Becker et al. 2005)		2.28
	(Murase et al. 2014)	${{\rm{\Gamma }}}_{\mathrm{SI}}=-2.0$ (BLR)	${\xi }_{\mathrm{CR}}\lt 12$
		${{\rm{\Gamma }}}_{\mathrm{SI}}=-2.0$ (blazar)	${\xi }_{\mathrm{CR}}\lt 21$
		${{\rm{\Gamma }}}_{\mathrm{SI}}=-2.3$ (BLR)	${\xi }_{\mathrm{CR}}\lt 153$
		${{\rm{\Gamma }}}_{\mathrm{SI}}=-2.3$ (blazar)	${\xi }_{\mathrm{CR}}\lt 241$
BL Lacs	(Mücke et al. 2003)	HSP (optimistic)	76.29
		LSP (optimistic)	5.78
	(Tavecchio et al. 2014)	HSP-dominated (1)	1.06
	a	HSP-dominated (2)	0.35
	(Tavecchio & Ghisellini 2015)	LSP-dominated	0.21
	(Padovani et al. 2015)	HSP (baseline)	0.75

Notes. 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 (Mücke et al. 2003) or to constraints on the baryon-to-photon luminosity ratios ${\xi }_{\mathrm{CR}}$ (Murase et al. 2014).

^aPredictions from Tavecchio et al. (2014) and Tavecchio & Ghisellini (2015) enhanced by a factor of 3 in correspondence with the authors.

Download table as:  ASCIITypeset image

In early models (before the year 2000) the neutrino flux per source is calculated to be directly proportional to the γ-ray flux in the energy range ${E}_{\gamma }\gt 100\ \mathrm{MeV}$ (Mannheim 1995) (A), ${E}_{\gamma }\gt 1\ \mathrm{MeV}$ (Mannheim 1995) (B), $20\ \mathrm{MeV}\lt {E}_{\gamma }\ \lt 30\ \mathrm{GeV}$ (Halzen & Zas 1997), and ${E}_{\gamma }\gt 100\ \mathrm{MeV}$ (Protheroe 1997). The γ-weighting scheme is therefore almost implicit in all these calculations, although the energy ranges vary slightly from the $100\ \mathrm{MeV}\mbox{--}100\ \mathrm{GeV}$ energy range used for the γ-weighting.

Among the newer models, only Padovani et al.'s (2015) uses a direct proportionality between the neutrino and γ-ray flux (for ${E}_{\gamma }\gt 10\ \mathrm{GeV}$), where the proportionality factor encodes a possible leptonic contribution. In all other publications a direct correlation to γ-rays is not used for the neutrino flux calculation. Since all these models assume that p/γ-interactions dominate the neutrino production, the resulting neutrino fluxes are calculated via the luminosity in the target photon fields. In Becker et al. (2005) the neutrino flux is proportional to the target radio flux, which in turn is connected to the disk luminosity via the model from Falcke & Biermann (1995). In Mücke et al. (2003) the neutrino flux is directly proportional to the radiation of the synchrotron peak. In Murase et al. (2014) the neutrino flux is connected to the X-ray luminosity, which in turn is proportional to the luminosity in various target photon fields. In Tavecchio et al. (2014) the neutrino luminosity is calculated using target photon fields from the inner jet "spine layer." However, a correlation to the γ-ray flux in these latter models may still exist, even in the case where leptonic γ-ray contributions dominate. This is mentioned in Murase et al. (2014), who explicitly predict the strongest γ-ray emitters to also be the strongest neutrino emitters, even though the model contains leptonically produced γ-ray emission. It should be noted that an independent IceCube analysis studying the all-flavor diffuse neutrino flux at PeV energies and beyond (Aartsen et al. 2016a) recently also put strong constraints on some of the flux predictions discussed in this section.

In this paper, we have analyzed all 862 Fermi-LAT 2LAC blazars and 4 spectrally selected subpopulations via an unbinned likelihood stacking approach for a cumulative neutrino excess from the given blazar directions. The study uses 3 years of IceCube data (2009–2012), which amount to a total of around 340000 muon-track events.

Each of the 5 populations was analyzed with two weighting schemes which encode assumptions about the relative neutrino flux from each source in a given population. The first weighting scheme uses the energy flux observed in γ-rays as weights, whereas the second scheme gives each source the same weight. This resulted in a total of 10 statistical tests, which were in turn analyzed in two different ways. The first is an "integral" test, in which a power-law flux with a variable spectral index is fitted over the full energy range that IceCube is sensitive to. The second is a differential analysis, in which 14 energy segments between ${10}^{2}$ and ${10}^{9}\ \mathrm{GeV}$, each spanning half a decade in energy, are fit independently with a constant spectral index of −2.

Nine of the ten integral tests show overfluctuations, but none of them are significant. The largest overfluctuation, a 6% p-value, is observed for all 862 2LAC blazars combined using the model-independent equal-weighting scheme. The differential test for all 2LAC blazars using equal source weighting reveals that the excess appears in the 5–10 TeV region with a local p-value of $2.6\sigma$. No correction for testing multiple hypotheses is applied, since even without a trial correction this excess cannot be considered significant.

Given the null results, we then calculated the flux upper limits. The two most important results of this paper are as follows:

• 1.  
  We calculated a flux upper limit for a power-law spectrum starting at $10\ \mathrm{TeV}$ with a spectral index of −2.5 for all 2LAC blazars. We compared this upper limit to the diffuse astrophysical neutrino flux observed by IceCube (Aartsen et al. 2015b). We found that the maximal contribution from all 2LAC blazars in the energy range between 10 TeV and 2 PeV is $27 \%$, including systematic effects and with minimal assumptions about the neutrino/γ-ray correlation in each source. Changing the spectral index of the tested flux to −2.2, a value allowed at about 3 standard deviations given the current global fit result (Aartsen et al. 2015b), weakens this constraint by about a factor of two. If we assume for each source a similar proportionality between the γ-ray luminosity in the 2LAC energy range and the neutrino luminosity, we can extend the constraint to the parent population of all GeV blazars in the observable universe. The corresponding maximal contribution is then around $10 \%$ from all GeV blazars, or 5%–10% from the other blazar subpopulations. In each case, we use the same power-law assumption as before in order to compare it to the observed flux. For FSRQs our analysis allows for a $7 \%$ contribution to the diffuse flux, which is in general agreement with a result obtained by Wang & Li (2015), who independently estimated that FSRQs do not contribute more than $10 \%$ to the diffuse flux using our earlier small-sample stacking result for 33 FSRQs (Aartsen et al. 2014b).
• 2.  
  We calculated upper limits using the γ-weighting scheme for 15 models of the diffuse neutrino emission from blazar populations found in the literature. For most of these models, the upper limit constrains the model prediction—by more than an order of magnitude for some of them. The implicit assumption in all these upper limits is a proportionality between the source-by-source γ-ray luminosity in the 2LAC energy range and its corresponding neutrino luminosity. All models published before the year 2000 and the model by Padovani et al. (2015) implicitly contain this assumption, although some of their energy ranges differ from the exact energy range in the 2LAC catalog. Even for the other models the proportionality assumption may still hold, as indicated by Murase et al. (2014).

Kadler et al. (2016) recently claimed a $5 \%$ probability for a PeV IceCube event to have originated close to blazar PKS B1424-418 during a high-fluence state. While $5 \%$ is not yet statistical evidence, our results do not contradict such single PeV-event associations, especially since a dominant fraction of the sensitivity of our analysis comes from the sub-PeV energy range. The same authors also show that the measured all-sky PeV neutrino flux cannot be compatible with an origin in a pure FSRQ population that has a peaked spectrum around PeV energies, as it would overpredict the number of observed events. Instead, one has to invoke additional assumptions—for example, a certain contribution from BL Lacs, leptonic contributions to the SED, or a spectral broadening of the arriving neutrino flux down to TeV energies due to Doppler shifts from the jets and the intrinsic redshift distribution of the blazars. Our results suggest that the last assumption, a spectrum broadening down to TeV energies, only works if the resulting power-law spectral index is harder than around −2.2, as the flux is otherwise in tension with our γ-weighting upper limit. A hard PeV spectrum is interestingly also seen by a recent IceCube analysis (Aartsen et al. 2016b) that probes the PeV range with muon neutrinos. Regardless of these speculations, the existing sub-PeV data require an explanation beyond the 2LAC sample from a yet unidentified galactic or extragalactic source class.

Our results do not provide a solution to explain the bulk emission of astrophysical diffuse neutrinos, but they provide robust constraints that might help to construct the global picture. Recently, Murase et al. (2015) argued that current observations favor sources that are opaque to γ-rays. This would, for example, be expected in the cores of AGNs. Our findings on the 2LAC blazars mostly provide a basis for probing the emission from relativistically beamed AGN jets and are in line with these expectations. We also do not constrain neutrinos from blazar classes that are not part of the 2LAC catalog—for example, extreme HSP objects. These sources may emit up to 30% of the diffuse flux (Padovani et al. 2016), and studies in this direction with other catalogs are in progress.

While the slight excess in the 5–10 TeV region is not yet significant, further observations by IceCube may clarify if what we see is an emerging soft signal or just a statistical fluctuation.

We acknowledge support from the following agencies: the U.S. National Science Foundation–Office of Polar Programs, the U.S. National Science Foundation–Physics Division, the University of Wisconsin Alumni Research Foundation, the Grid Laboratory of Wisconsin's grid infrastructure at the University of Wisconsin–Madison, and the Open Science Grid's grid infrastructure; the U.S. Department of Energy's National Energy Research Scientific Computing Center and the Louisiana Optical Network Initiative's grid computing resources; the Natural Sciences and Engineering Research Council of Canada, WestGrid, and Compute/Calcul Canada; the Swedish Research Council, the Swedish Polar Research Secretariat, the Swedish National Infrastructure for Computing, and the Knut and Alice Wallenberg Foundation, Sweden; the German Ministry for Education and Research, Deutsche Forschungsgemeinschaft, the Helmholtz Alliance for Astroparticle Physics, and the Research Department of Plasmas with Complex Interactions (Bochum), Germany; the Fund for Scientific Research, FWO Odysseus program, Flanders Institute (to encourage scientific and technological research in industry) and the Belgian Federal Science Policy Office; the University of Oxford, United Kingdom; the Marsden Fund, New Zealand; the Australian Research Council; Japan Society for Promotion of Science; the Swiss National Science Foundation, Switzerland; the National Research Foundation of Korea (NRF); and Villum Fonden, Danish National Research Foundation, Denmark.

The equal-weighting limits use source count distributions to model the neutrino injection. The SCD serves as a PDF template from which relative neutrino injection weights are drawn. Depending on the shape of the SCD and the range of flux values in which the SCD is being sampled, the resulting central neutrino upper limit value shifts, and the range of calculated flux upper limits broadens. This is illustrated in Figure 7, which shows the upper limits derived from ensemble simulations drawn from four different SCDs. It also includes standard 6-year point source constraints (Aartsen et al. 2015c) for similar flux realizations.

Zoom In Zoom Out Reset image size

All examples are shown for a population of 862 objects, the size of the "all 2LAC blazars" sample. In the left panel ("extended"), the SCD template is the measured blazar γ-ray SCD from Abdo et al. (2010c) and is extrapolated by five orders of magnitude to lower flux values. The minimum flux value is arbitrarily chosen, but it is small enough such that the distribution is a scale-free power law and an extension toward even smaller flux values does not make a difference in terms of average sample outcomes. In the second panel ("standard"), the SCD is exactly the Fermi-LAT blazar SCD from Abdo et al. (2010c) and spans three orders of magnitude in flux. In the third panel the SCD is of Euclidean form and is extended over a flux range such that the cumulative SCD equals the number count of the "standard" distribution in the second panel. In the fourth panel the SCD is a delta distribution, which gives an equal weight to each source, i.e., the assumption that is used in the weighting of the PDF for the statistical test. The lower row displays the respective 90% central interval for upper limit outcomes from this analysis and constraints from 6-year single-point-source search discovery potentials.

As we sample the relative source contributions from a growing flux range, corresponding to the column order $4\to 3\to 2\to 1$, flux upper limit variations increase, and the stacking analysis constraint weakens. At the same time, the constraints from the single-point-source search become stronger with a single source increasingly dominating the total population.

The first and fourth columns correspond to limiting cases, neither of which is appropriate to use for this analysis, but they are just shown to illustrate the general behavior of the procedure. The delta-peak SCD (4th column) is unphysical, since it corresponds to an equal flux per source. The extrapolated SCD (1st column) yields an extreme spread of signal contributions, which roughly corresponds to a random draw from the entire population in the universe, assuming that the faint end corresponds to the weakest blazars that can in principle be detected. Since the random draw mostly consists of sources from the faint-flux end, it corresponds to a situation where the neutrino flux is anti-proportional to the γ-ray flux, which is unphysical. All results in this paper make use of the 2nd-column SCD. A cross-over of point source constraints and constraints from this stacking analysis is reached for realizations drawn from a dN/dS distribution that lies between the SCDs from the 1st and 2nd columns. For illustrative purposes we also include a particular flux scenario with equal intrinsic luminosity, where the neutrino flux per source is proportional to $1/{{d}_{L}}^{2}$.⁶⁶ Missing redshifts for BL Lacs are drawn randomly from the BL Lac distribution where the redshifts are known, taking into account the synchroton peak information. The results for two such realizations of missing-redshift sources lie in the range of SCD draws from the 2nd column. Since the two realizations do not differ significantly, we conclude that the resulting estimate is robust, even though some of the BL Lac redshifts are unknown.

We determine the energy range for which the upper limit is supported by IceCube data based on the differential sensitivity. The procedure is illustrated in Figure 8 for three different generic power-law spectra using the γ-weighting scheme. The ratio of the differential sensitivity curve to a convex energy spectrum generally forms a function with a single maximum that falls off toward the sides. The central interval enclosing 90% of the area under the ratio function is used to define the energy range that contributes 90% to the sensitivity for a given energy spectrum. This area–sensitivity relation has been checked empirically. The methodology is an extension to a previous method used in Aartsen et al. (2014b), which only uses the 90% central interval of signal events and thereby neglects the background rate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The γ-weighting scheme implicitly assumes a proportionality between the neutrino and γ-ray luminosities. In this case, the fraction of the total neutrino flux of the population that originates from the source resolved in γ-rays is equal to the fraction of the total γ-ray flux that originates from the resolved sources. It has been estimated that 70% of the total diffuse γ-ray flux from blazars between 100 MeV and $100\mathrm{GeV}$ has been resolved in 1FGL blazars at high galactic latitudes $| b| \gt 15^\circ$ (Ajello et al. 2015), particularly in 1FGL FSRQs (Ajello et al. 2012). The 2LAC catalog used in this work contains blazars at galactic latitudes $| b| \gt 10^\circ$. At extragalactic latitudes ($10^\circ \lt | b| \lt 15^\circ$) the detection efficiency might be worse. However, even if it unrealistically sharply drops to zero, one can estimate that the total resolved fraction only shrinks by an amount that is proportional to the ratio of the $10^\circ \lt | b| \lt 15^\circ$ sky fraction ($\approx 8 \%$ of the total sky) with respect to the rest ($\approx 74 \%$ of the total sky ), i.e., to $\tfrac{0.08\cdot 0+0.74\cdot 70 \% }{0.82}\approx 63 \%$. Since this estimate is conservative and still within the error on the quoted 70% value and since the 2LAC sample is based on the more sensitive 2FGL catalog, we conclude that 70% completeness 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\approx 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 $\gamma \gamma$-cascades due to intergalactic magnetic fields (Aharonian et al. 1994), the fraction of neutrinos emitted from 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 $\approx 2.8$. The total extragalactic gamma-ray background has an intensity of $11.3\times {10}^{-6}\ \mathrm{photons}/{\mathrm{cm}}^{2}\ {\rm{s}}\ \mathrm{sr}$, of which $4.1\times {10}^{-6}\ \mathrm{photons}/{\mathrm{cm}}^{2}\ {\rm{s}}\ \mathrm{sr}$ originates from resolved blazars and the other $7.2\times {10}^{-6}\ \mathrm{photons}/{\mathrm{cm}}^{2}\ {\rm{s}}\ \mathrm{sr}$ from the IGRB (isotropic γ-ray background). If we assume that the entire contribution to the IGRB stems from EBL-induced $\gamma \gamma$-cascades from γ-rays emitted by blazars, i.e., unresolvable isotropized γ-rays, the resulting ratio between the total emission from blazars and the emission from resolved blazars would be $11.3/4.1\approx 2.8$. Since the intensity of the IGRB can be well explained by contributions from unresolved—but in principle resolvable—blazars, from star-forming galaxies, and from radio galaxies (Ajello et al. 2015), we deem this maximum cascade emission scenario unlikely and use a factor of 1.43 throughout this work.

Figures 9 and 10 show the individual weight distribution in declination, p-values and upper limits for all blazar populations. This information is supplementary material for Section 4.2.1 and Section 5.