Analysis of the cumulative neutrino flux from Fermi-LAT blazar populations using 3 years of IceCube data

The recent discovery of a diffuse neutrino flux up to PeV energies raises the question of which populations of astrophysical sources contribute to this diffuse signal. One extragalactic candidate source population to produce high-energy neutrinos are Blazars. We present results from a likelihood analysis searching for cumulative neutrino emission from Blazar populations selected with the 2nd Fermi-LAT AGN catalog (2LAC) using an IceCube data set that has been optimized for the detection of individual sources. In contrast to previous searches with IceCube, the investigated populations contain up to hundreds of sources, the biggest one being the entire Blazar sample measured by the Fermi-LAT. No significant neutrino signal was found from any of these populations. Some implications of this non-observation for the origin of the observed PeV diffuse signal will be discussed.


Introduction
In 2013, the IceCube collaboration reported on the detection of an astrophysical diffuse neutrino flux in the 40 TeV -2 PeV energy range [1]. Since then it has been measured down to about 10 TeV in energy and the spectrum has been shown to be consistent with a power law of index −2.5 [2]. The origin of the observed neutrino signal is so far unidentified. A prominent candidate is the population of Blazars, active galactic nuclei (AGN) with their relativistic jets pointing along the line-of-sight. Neutrino emission from Blazars has been discussed extensively in hadronic jet models originating from p-p or p-γ interactions (e.g. [3]) which also produce gamma rays from neutral pion decays. The extragalactic non-thermal radiation background at GeV energies and above (EGB) is dominated by the emission from Blazars [4]. An experimental approach is presented here to answer the question if the Blazars that dominate the high energy EGB emission also produce a large fraction of the astrophysical neutrino flux. The idea is to look for directional clustering in a large sample of muon tracks collected by IceCube, around the directions of gamma-ray sources associated with Blazars in the Fermi-2LAC catalog [5].
This analysis is performed as a maximum likelihood stacking analysis looking for a cumulative signal from multiple sources. In contrast to previous stacking searches which looked at 10-30 sources [6] [7], the populations in this work usually comprise more than 100 sources, the biggest population being 862 sources when all Blazars from the 2LAC catalog are combined. Any search for a cumulative signal from multiple sources needs to make assumptions about the relative contributions ("weights") a e-mail: thorsten.gluesenkamp@desy.de arXiv:1502.03104v1 [astro-ph.HE] 10 Feb 2015 The Journal's name of the individual sources within the population. Here, this weighting scheme is set up in a way to be minimally biased with respect to γ−ν correlation assumptions. Section 2 covers the motivation for this type of search and the technical setup. In Section 3 the sample of muon tracks used in this analysis is described. In Section 4 the results of this search are summarized. Conclusions are presented in Section 5.

Motivation and Method
In previous stacking searches for a cumulative neutrino signal from Blazars in IceCube data [6] [7], the measured gamma-ray flux was taken as an estimate for the expected neutrino flux. However, the γ and ν spectra may not be strongly correlated: Other contributions to the gamma flux almost surely exist, e.g. inverse Compton photons [8] from high-energy electrons or synchrotron photons from muons and protons [9]. Moreover, radiation fields at the sources and in the interstellar medium absorb and reprocess the gamma rays originally produced together with the neutrinos. In order to overcome this bias as much as possible, two modifications with respect to previous IceCube stacking searches have been implemented: • There is no sub-selection of sources within the catalog based on their gamma-ray flux The populations tested for a neutrino signal are only defined via their spectral classification, the largest one being just all Blazars of the 2LAC catalog. Due to the competing production processes and different environmental conditions described above, a source with a relatively low gamma-ray flux could still contribute significantly to the neutrino signal.
• Two complementary weighting schemes are used, including a uniform relative weight for all sources Note that F γ is the energy flux. Gamma-rays produced in neutral pion decays at multi-TeV energies (i.e. by the same process that also produces a neutrino signal) are absorbed quickly by IR photons present in the source environment and reprocessed to lower energies. In this process they would approximately maintain their bolometric luminosity.
Equal weighting is used as the most unbiased weighting scheme, completely ignoring potential correlations. Although it is in a strict sense unrealistic, because no population consist of sources which all emit equally strong, the resulting limits are quite independent of the unknown correlation between neutrino and gamma-ray flux. They even hold for the case that no correlation exists.
The analysis is performed with an unbinned maximum likelihood Ansatz. The log-likelihood function is constructed as where N is the total number of neutrino events, S i and B i are the signal probability distribution function (PDF) and background PDF evaluations for event i. The background PDF is taken from the data distribution in declination and reconstructed energy and is uniform in right ascension. The signal PDF is defined as RICAP-14 The Roma International Conference On Astroparticle Physics where j denotes the individual sources within the population. The term S j (x i ; σ i ) denotes the point spread function (PSF) of source j evaluated at celestial coordinate x i , depending on the event-specific angular error estimate σ i , while the term ε j (E i ) denotes the energy PDF of source j evaluated for the reconstructed energy E i . The PSF is locally modeled as an analytic 2-d radially symmetric gaussian while the energy PDF ε is estimated from Monte Carlo simulations. Each individual source within the signal PDF is weighted with a declination-dependent detector response weight w dec. and the relative source weight w source described above. Equation 1 is minimized with respect to n s and the likelihood ratio with respect to n s = 0 is calculated for each population and each weighting scheme. It is tested if the observed likelihood ratio is compatible with statistical fluctuations of the background events by repeating the process for many randomized sky maps. A randomized sky map is generated by randomizing the right ascension of the observed events (a more detailed description of this procedure is found in [7]). In case no significant deviations from background fluctuations are observed, upper limits on the neutrino flux of the tested population are calculated via the CL s method [10]. This is performed twice, once for each weighting scheme.

Muon Track Data
A muon-track selection for the years 2009-2011 (IC59-IC86/1) is used, which is similar to the one described in [7] 1 . The median angular resolution of muon-neutrino tracks with this selection is around 1 • at 1 TeV and 0.5 • at 100 TeV. One additional cut is applied which keeps only events that have an estimated angular resolution better than 5 degrees. In the southern hemisphere (downgoing events in IceCube) this sample is dominated by muons from air showers, while in the upgoing region the majority of events come from atmospheric muon-neutrinos.

Populations tested
The populations to be tested are based on the 2nd Fermi AGN catalog (2LAC) [5] and are classified purely by spectral properties. Only sources that do not suffer from potential confusion are considered (CLEAN flag set to True). The populations are listed in table 1. It should be noted that some of the populations are overlapping and the p-values correlated (the largest overlap exists between the LSP and FSRQ populations, which is about 60 %).

Results
For each population two tests were performed, one for each of the two different source weighting schemes. Table 2 summarizes the results 2 . All test outcomes are compatible with statistical fluctuations of the background. The smallest p-value is 6 % for the "All 2LAC Blazar" sample when tested with the "equal weighting" scheme. Figure 1 shows the resulting neutrino flux upper limits for an E −2.5 spectrum for the "All 2LAC-Blazars" population and for the FSRQ population. Notice that the flux upper limit using the equal-weighting scheme is represented as a band, instead of a line. It serves as an indicator of how robust this flux limit is towards changes of the actual source count distribution over the sky, which is unknown. Random realizations are created assuming a source count distribution for the neutrino emission that has the same shape as the source count distribution in gamma rays  1 2 62 Motivated by work in [9] 1 FSRQ/BL-LAC based on optical line equivalent width 2 Low/Intermediate/High Synchroton Peaked Object (LSP/ISP/HSP), based on position of synchroton peak   Figure 1: Neutrino flux upper limits for an E −2.5 spectrum (blue) compared to diffuse bestfit (black solid) from [2]. a) shows this comparison for the "All 2LAC Blazar" sample, b) for the FSRQs. Percentages with arrows denote the fraction with respect to the diffuse flux.
RICAP-14 The Roma International Conference On Astroparticle Physics (parametrization from [14]). For each source in the population a relative source strength is drawn from the source count distribution and assigned to the source. The band encompasses the central 90 % of the flux upper limits found from different such realizations. In comparison to the limits the best fit diffuse flux above 10 TeV is shown from [2] with a spectral index of −2.46. Assuming an isotropic flux and 1/1/1 composition of flavors, the maximal contribution from the Blazar populations with respect to the diffuse flux can be read off. This contribution to the diffuse neutrino flux is at most around 20 %, even if there is a very weak correlation with the observed gamma flux. For FSRQ's, these numbers are lower as the population size is smaller.