Correlating High-energy IceCube Neutrinos with 5BZCAT Blazars and RFC Sources

The quest for the sources of the highest-energy cosmic radiation is closely connected with the search for the production sites of high-energy astrophysical neutrinos, generated in cosmic-ray interactions. In this pursuit, the IceCube Neutrino Observatory ⁶ (Aartsen et al. 2017c) has been playing a leading role for over a decade, from the detection of hundreds of astrophysical neutrinos with energies reaching above the PeV (10¹⁵ eV), which marked the birth of high-energy (TeV–PeV) neutrino astrophysics (IceCube Collaboration 2013), to the recently reported evidence for neutrino emission from the Seyfert II galaxy NGC 1068 (IceCube Collaboration et al. 2022). The first association with an extragalactic source of high-energy neutrinos happened in 2017 September, when IceCube detected an event with very high energy (∼290 TeV) and a high probability of having astrophysical origin. Follow-up multimessenger observations by the Fermi-LAT ⁷ and MAGIC ⁸ Collaborations confirmed that the neutrino event had a γ-ray counterpart in the blazar ⁹ TXS 0506+056, which was also exhibiting an enhanced γ-ray activity at the time of the neutrino detection (IceCube Collaboration et al. 2018a). The chance coincidence probability of correlated neutrino and γ-ray emissions was found to be at the 3σ level. An a posteriori analysis of IceCube archival data at the position of TXS 0506+056 detected a neutrino flare between 2014 September and 2015 March incompatible with the background at the 3.5σ level (IceCube Collaboration et al. 2018b). The TXS 0506+056 association immediately triggered great interest in the hypothesis of a correlation between blazars and astrophysical neutrinos, which had been already suggested in several studies antecedent to the 2017 September event (e.g., Padovani & Resconi 2014; Padovani et al. 2016; Huber & Krings 2017; Lucarelli et al. 2017, and references therein). Since then, many studies have addressed this topic, trying to find statistically robust associations, mostly with inconclusive (≲3σ) results (see, e.g., the review by Giommi & Padovani 2021; Abbasi et al. 2023a, and references therein). Buson et al. (2022a) have recently cross-matched the 5th edition of the Roma-BZCAT Multifrequency Catalogue of Blazars (5BZCAT: Massaro et al. 2015) with a neutrino p-value map covering 7 yr of observations made by IceCube between 2008 and 2015 (Aartsen et al. 2017b). The neutrino sky published by IceCube maps the local probability of the observation to be consistent with the background only, at each location in the sky. Due to its unique location at the geographic South Pole, neutrinos coming from the northern sky have to cross Earth to reach the IceCube detector. Above PeV energies, they start to be absorbed in Earth at declinations δ > 30°. Therefore, Buson et al. (2022a) focused only on the southern sky (i.e., −85° < δ < − 5°), arguing that this region is the most sensitive to an astrophysical signal in the PeV–EeV energy range. They came to the conclusion that blazars were associated with high-energy astrophysical neutrinos with a post-trial chance probability of 6 × 10⁻⁷, later scaled down to 2 × 10⁻⁶ (Buson et al. 2022b), equivalent to ∼4.5σ. The main purpose of this paper is to test the results of Buson et al. (2022a, 2022b) and expand their study first to the 7 yr northern sky and then to the whole sky using a recently released, ∼60% larger and improved, IceCube data set covering 10 yr of observations (IceCube Collaboration et al. 2021).

5BZCAT contains coordinates and multifrequency data for 3561 sources, either confirmed blazars or exhibiting blazar-like characteristics. When it was first released, this list of blazars was likely the most comprehensive one available, as it included the vast majority of sources well documented in the literature. By construction, however, 5BZCAT is a very inhomogeneous sample, as it collects sources detected in various frequency bands and different surveys and is not flux limited, i.e., it does not include all objects above a certain flux in a given band. Following Buson et al. (2022a), we exclude sources overlapping with the Galactic plane (Galactic latitude ∣b_II∣ < 10°), consider only ∣δ∣ < 85°, and do not include the 92 objects classified as "BL Lac candidates" in 5BZCAT.

To assess the impact of using a flux-limited and uniformly distributed sample, as compared to an inhomogeneous list of targets, we also use a subset of the very long baseline interferometry (VLBI) based Radio Fundamental Catalogue ¹⁰ (RFC), which is a radio-selected sample of 3411 active galactic nuclei (AGNs) with an average 8 GHz VLBI flux density ≥150 mJy used by Plavin et al. (2021). This was originally selected by Plavin et al. (2020), who noted that a special effort was made to ensure completeness above this limit. (As shown in Section 6.2, the RFC sources are distributed more uniformly in the sky than the ones in the 5BZCAT, especially in the northern sky.) For consistency, we use here the same version as in the latter paper (rfc_2020b), which tested for an association between the RFC and IceCube neutrinos (see Section 6.2). As done for 5BZCAT, we exclude RFC sources with ∣b_II∣ < 10° and ∣δ∣ > 85°. The RFC is dominated by Doppler-boosted AGNs, i.e., blazars. In fact, when cross-correlating it with a "master blazar sample" composed of the blazars in the 5BZCAT, 3HSP (Chang et al. 2019), and Fermi 4FGL-DR3 (Abdollahi et al. 2022) catalogs, which includes 6425 unique sources, we find 59% matches. Since not all blazars in the Universe have been identified, this provides a robust lower limit to the RFC blazar fraction. This percentage increases to ≳84% when we make a cut at 0.5 Jy instead of 0.15 Jy, which shows that brighter radio blazars are more likely to have been identified.

The IceCube Collaboration has recently released 10 yr of candidate neutrino events detected between 2008 April 6 and 2018 July 10 (IceCube Collaboration et al. 2021). ¹¹ The sample is optimized for the search of neutrino point sources, thus including track-like events, primarily due to muon neutrino candidates (muon tracks), from all directions in the sky. In comparison to the previous 7 yr sample (Aartsen et al. 2017b), utilized to generate the neutrino p-value sky map used by Buson et al. (2022a), the number of events is increased by ∼60%, and an updated processing procedure and a new angular reconstruction are applied to the events recorded between 2012 April and 2015 May (and to the events recorded afterward). Notably, the enhanced directional reconstruction improved the average angular resolution by more than 10% for events with reconstructed muon energies above 10 TeV (Carver 2019; Aartsen et al. 2020). In the northern sky, the range of reconstructed energies goes from a few tens of GeV to ∼1 PeV, with most events reconstructed at 1 TeV. The southern sky, instead, covers a wider spectrum of energies, going up to 10 PeV. Indeed, the absorption of the highest-energy neutrinos in Earth reduces the detection efficiency of IceCube for neutrino energies above 1 PeV at δ ≳ 30°. On the other hand, the southern neutrino sky is dominated by the background of muons and muon bundles created in atmospheric cosmic-ray showers—resembling high-energy tracks (Aartsen et al. 2016a)—which are completely absorbed above δ = 5°. These considerations leave the horizon (−5° ≲ δ ≲ 30°) as the region where IceCube is the most sensitive to an astrophysical flux of neutrinos.

The 1,134,450 events in the released 10 yr sample were previously used by the IceCube Collaboration for a time-integrated search for neutrino point sources in both hemispheres (Aartsen et al. 2020). A machine-readable version of the neutrino p-value map in Aartsen et al. (2020) is not publicly available. However, the data release includes the experimental data with their per-event measured physical properties (time, direction, and energy), the binned effective areas, and the binned instrument responses for five data-taking seasons. The effective areas are published as 2D histograms in neutrino energy and decl., while the detector response functions are provided as 5D histograms mapping a neutrino to the probability of the reconstructed observables in the detector. With this information it is possible to analyze the data to search for astrophysical neutrino emission. We produce our own neutrino p-value map with the public 10 yr sample using an unbinned maximum likelihood ratio method, following the prescription outlined in Braun et al. (2008) (see Appendix A for details). However, it is worth mentioning that only three bins for the true incoming neutrino decl. are provided in the detector response matrices, covering the sky from −90° to −10°, from −10° to +10°, and from +10° to +90°, respectively. The coarseness of this information causes a decrease of ∼20%–50% in the sensitivity of the analysis to an astrophysical neutrino flux with respect to what is published in Aartsen et al. (2020). This discrepancy may be attributed to the provided detector response being averaged over wide decl. ranges, where it undergoes considerable variations (see Appendix B). It should also be noted that the likelihood formalism adopted in this work does not exactly conform to the one used by the IceCube Collaboration (see Appendix A for details on the differences). Therefore, the neutrino p-value map based on the 10 yr data release produced within this work should not be considered an exact reproduction of the one published in Aartsen et al. (2020), but should instead be seen as the outcome of a likelihood-based search for neutrino sources that finds results mostly statistically compatible with what has been published by the IceCube Collaboration (see Appendix B).

From the 10 yr neutrino p-value map that we obtain from the likelihood scan of the sky, we can then select the most interesting neutrino spots for the cross-correlation analyses described in the following section.

We perform various positional correlation analyses similar to the one described in Buson et al. (2022a). Their test correlates the sources in 5BZCAT (see Section 2) with the southern neutrino sky by using 7 yr of muon tracks recorded by the IceCube detector (Aartsen et al. 2017b). These have been published in the form of pre-trial local p-values on a grid of equivalent pixels of ∼(011)² constructed using a HEALPix (Górski et al. 2005) projection with Nside parameter 512. ¹² We follow Plavin et al. (2021) and Buson et al. (2022a) and denote the negative logarithm of local p-values mapping the neutrino sky as $L=-{\mathrm{log}}_{10}p$ to avoid confusion with the p-values of the correlation analyses. The L values are based on a maximum likelihood ratio comparing a background and a signal hypothesis, where in the former case the neutrino emission is due to atmospheric background and diffuse astrophysical emission, while in the latter an additional component comes from the source and clusters around it. Hence, a larger value of L indicates a lower probability that the neutrino data are consistent with the background-only hypothesis (see Appendix A).

We extend the analysis to the 7 yr northern sky. Furthermore, we search for correlations between the same neutrino sample and the RFC sources (see Section 2). Finally, we increase the size of the neutrino sample by using an updated neutrino sky map based on the 10 yr data release (see Section 3) and search for correlations between the newer neutrino sample and the two aforementioned catalogs.

As we want to follow the same statistical method developed in Buson et al. (2022a), we also search for a spatial correlation between the blazars and the IceCube neutrinos using the hot spots in the neutrino sky map of L values. The hot spots are defined as all independent pixels with local significance L above some predefined threshold ${L}_{\min }$. For this type of search, the hot spots have to be selected so as not to double-count possible interesting neutrino signals. In fact, the pixels in a sky map are necessarily correlated because of the same events contributing to the likelihood at neighboring locations. To identify independent hot spots, we select pixels above the chosen significance threshold whose centers are at least 1° apart from each other (Aartsen et al. 2017b, 2019). The test statistic (TS) of the correlation analysis is the number of hot spots with significance above a certain local threshold ${L}_{\min }$ having at least one source closer than some association radius r_assoc. A set of ${L}_{\min }$ and r_assoc values is defined a priori, and a TS is computed for each possible combination. For each pair of the ${L}_{\min }$ and r_assoc parameters, the p-value of the correlations is defined as the chance probability to get the observed TS value. We derive it by running the correlation analysis on pseudo-experiments generated by randomizing the position of the sources within 10° from their original position. In this way, the original structure of the catalog is preserved in the simulated ones (see Buson et al. 2022a for more details). The combination of the analysis parameters yielding the smallest p-value is the best-fit result. The best p-value then needs to be corrected for the "look-elsewhere effect," i.e., for having tested all possible combinations of ${L}_{\min }$ and r_assoc. We choose different values of the analysis parameters depending on the neutrino sample that is being tested, as detailed in the following.

4.1. Correlation with the 10 yr Neutrino Sample

For the correlation test using the 10 yr neutrino L map, we redefine the analysis parameters ${L}_{\min }$ and r_assoc. Given the reduced sensitivity of our point-source analysis compared to the one published in Aartsen et al. (2020), we decided to start scanning the ${L}_{\min }$ parameter from a lower significance of 3.0, instead of 3.5, in bins of 0.5 up to 4.5. As for r_assoc, we chose its range based on the capability of our point-source likelihood analysis to localize a neutrino source in the sky. To estimate this localization performance, we simulate astrophysical neutrino sources emitting with an E⁻² energy spectrum at various locations in the sky. For each simulation, we search the associated hot spot in a circle of 1° radius centered at the source location and calculate the angular distance between the two. Figure 1 displays the distribution of these angular distances. We do it separately in the northern and southern sky, as we expect different outcomes owing to the different energy ranges and data selection procedures used for the two hemispheres (IceCube Collaboration et al. 2021). To simulate the sources, signal events are injected by sampling them from the detector response matrices, which map the true incoming neutrino energy and decl. to the distribution of reconstructed muons that the detector would see. This procedure intrinsically takes into account the smearing due to the detector resolutions for direction and energy reconstructions. To locate the neutrino hot spot, we divide the circle around the source into pixels of ∼(011)², maximize the likelihood ratio (see Appendix B), and calculate the local significance L at each of them. The pixel with the maximum local L value is taken as the hot spot for that simulation. By repeating this procedure ∼10,000 times in each hemisphere, we find that 50% (90%) of the hot spots with L > 3.0 are localized closer than 017 (080) from the injection location in the northern sky, ¹³ while in the southern sky 50% (90%) of the selected hot spots are closer than 013 (035) from the simulated source position (see Figure 1). A smaller association region corresponds to a reduced chance probability of associating a neutrino hot spot with a source. Hence, to always include at least 50% of the distributions and maximally reduce the background, we scan r_assoc from 020 and go up to 070 (the maximum r_assoc values used by Buson et al. 2022a) in steps of 005.

Zoom In Zoom Out Reset image size

5.1. 5BZCAT−7 yr Hot Spot Correlation

Figure 2 shows the 5BZCAT−hot spot spatial correlation local (pre-trial) p-value versus r_assoc for the 7 yr southern sky (left panel) using the same ${L}_{\min }$ values as in Buson et al. (2022a). This is almost exactly the same as their Figure 1. We also derive the same minimum post-trial significance (∼2 × 10⁻⁶, ∼4.5σ) as Buson et al. (2022b) for ${L}_{\min }=4.0$ and r_assoc = 055. The associated blazars are the same 10 as those listed in Table 2 of Buson et al. (2022a), i.e., 5 flat-spectrum radio quasars (FSRQs), 3 BL Lacs, and 2 blazars of uncertain type, with relatively large mean angular separations between the blazar and hot spot positions 〈ψ〉 ∼ 040 and mean redshift 〈z〉 ∼ 1.3. All of them are by definition radio detected, and three are γ-ray sources.

Zoom In Zoom Out Reset image size

In the 7 yr northern sky, 112, 43, 18, and 7 neutrino hot spots have $L\gt {L}_{\min }$, with ${L}_{\min }$ ranging from 3.0 to 4.5 in steps of 0.5 (see Section 4.1). The correlation analysis results are compatible with the background (p-value >0.1%, which is equivalent to a significance lower than 3σ), with a minimum local p-value ∼1% (Figure 2, right panel) and a post-trial one of ∼10%. Note that in this case, based on the results shown in Sections 4 and 4.1, r_assoc starts from 02 and the additional significance threshold ${L}_{\min }=3.0$ is used.

5.2. RFC−7 yr Hot Spot Correlation

Figure 3 shows the RFC−neutrino hot spot spatial correlation local (pre-trial) p-value versus r_assoc for the 7 yr southern sky (left panel; using the same ${L}_{\min }$ and r_assoc thresholds as in Buson et al. 2022a), which shows a minimum local p-value 1.4 × 10⁻⁴, which becomes 6.9 × 10⁻⁴ (3.2σ) post-trial, achieved for ${L}_{\min }=4.5$ and r_assoc = 055, the latter being the same value as for 5BZCAT. This is easily explained by the fact that the RFC and 5BZCAT sources responsible for this minimum are the same. Since for 5BZCAT ${L}_{\min }$ was 4.0, the associated blazars in this case are fewer, only five, i.e., four FSRQs and one blazar of uncertain type, with 〈r_assoc〉 ∼ 041 and 〈z〉 ∼ 1.8. Only one of them is a γ-ray source. However, similarly to the 5BZCAT case, in the 7 yr northern sky (Figure 3, right panel) the results are not significant, with a minimum local p-value ∼3% and a post-trial one of ∼17%.

Zoom In Zoom Out Reset image size

5.3. 5BZCAT−10 yr Hot Spot Correlation

In the 10 yr southern sky, 149, 58, 18, and 7 hot spots have $L\gt {L}_{\min }$ with ${L}_{\min }$ in [3.0, 3.5, 4.0, 4.5], while the northern sky has 95, 36, 7, and 2 hot spots with significance above the same set of thresholds. We note that while Aartsen et al. (2020) reported a local significance L = 3.7 at the location of the blazar TXS 0506+065, this object has a lower significance L = 2.7 in this analysis. This is interpreted as being due to the detector response matrix, which averages the detector smearing of the incoming neutrinos in a bin ranging from −10° to +10°, where IceCube's data selection and processing change (see Appendix B for more details). Having a significance below the minimum tested threshold of ${L}_{\min }=3.0$, TXS 0506+056 is not included in our correlation analysis. However, even the inclusion of a hot spot coincident with TXS 0506+056 does not make the correlation of either of the two tested catalogs incompatible with a chance coincidence (see Appendix C).

Figure 4 shows the 5BZCAT−hot spot spatial correlation local (pre-trial) p-value versus r_assoc for the 10 yr southern (left panel) and northern (right panel) skies. The results are not significant, with minimum local (post-trial) p-values ∼25% (∼75%) and ∼3% (∼19%), respectively.

Zoom In Zoom Out Reset image size

5.4. RFC−10 yr Hot Spot Correlation

Figure 5 shows the RFC−hot spot spatial correlation local (pre-trial) p-value versus r_assoc for the 10 yr southern (left panel) and northern (right panel) skies. Again the results are not significant, with minimum local (post-trial) p-values ∼32% (∼80%) and ∼7% (∼34%), respectively.

Zoom In Zoom Out Reset image size

The results from the different correlation analyses described in this section are summarized in Table 1.

Table 1. Summary of the Results

	Southern Sky		Northern Sky
	7 yr ν	10 yr ν	7 yr ν	10 yr ν
5BZCAT	2.5 × 10−6	0.75	0.10	0.19
RFC	6.9 × 10−4	0.80	0.17	0.34

Note. Post-trial p-values for the different positional correlation analyses performed in this work.

Download table as:  ASCIITypeset image

6.1. 5BZCAT

We have reproduced the very significant result of Buson et al. (2022a, 2022b) on the neutrino hot spot−5BZCAT blazar correlation using 7 yr of IceCube data in the southern sky. However, this apparently strong result is not robust since the correlation becomes insignificant when applying the same type of analysis to the 7 yr northern sky. Moreover, when using the larger and improved 10 yr IceCube sample covering the whole sky, no correlation is found in either hemisphere.

Given the nature of the 5BZCAT catalog (Section 2), its utilization for a statistical test in place of other available flux-limited samples, such as the RFC catalog, is an arguable choice. Moreover, being 9 yr old, 5BZCAT is also by now quite incomplete: for example, while the latest incremental Fermi-4LAT catalog (Abdollahi et al. 2022) contains more than 3700 blazars, 5BZCAT includes only ∼1000 γ-ray-emitting sources. Indeed, the left panel of Figure 6 shows that the decl. distribution of 5BZCAT is very different from uniform, as the Southern Hemisphere is severely undersampled, while there are two notable excesses at δ ∼ 0° and δ ∼ 40°. The former point is well-known and due to a radio−optical identification bias toward the northern sky (e.g., Massaro et al. 2015). We test the hypothesis that the decl. distribution of the sources is a realization of the uniform distribution (excluding the Galactic plane) using the Kolmogorov–Smirnov (K-S) and Wilcoxon–Mann–Whitney (WMW) tests. Both tests result in strong incompatibility, with a significance larger than 8σ.

Zoom In Zoom Out Reset image size

This situation is similar to that of Pierre Auger Collaboration et al. (2008), in which the authors had found an excess (chance probability 1.7 × 10⁻³) in the correlation of ultra−high-energy cosmic rays (UHECRs) above 57 EeV lying within 31 of sources closer than 75 Mpc in the AGN catalog of Véron-Cetty & Véron (2006; this catalog is also a very inhomogeneous list of sources). However, the fraction of arrival directions that were fulfilling the criteria above went from ${69}_{-13}^{+11} \%$ (to be compared with 21% expected for isotropic UHECRs), down to ${38}_{-6}^{+7} \%$ (Abreu et al. 2010), and then to 33% ± 5% and finally ${28.1}_{-3.6}^{+3.8} \%$ (Aab et al. 2015). To quote the latter paper, "the high level of correlation found initially was probably affected by a statistical fluctuation." We have shown here that this is also the case for the IceCube−5BZCAT correlation.

The distribution of the angular distance between the hot spots and the 10 associated sources in Buson et al. (2022a) deserves some scrutiny. The mean offset is 04, which is ∼3 times larger than the median localization capability of the point-source analysis used here to analyze the 10 yr neutrino sample in the southern sky (see Section 4.1 and Figure 1). Seven out of 10 hot spots are found at an angular distance ψ ≥ 04 from the source they are associated with. The binomial probability that this happens using our point-source analysis on the 10 yr sample is ∼2 × 10⁻⁶. With respect to the 7 yr data used by Buson et al. (2022a, 2022b), the average angular resolution on the single event in the 10 yr data set has improved by ∼10% for the past 6 yr only, while the first 4 yr are treated exactly as before. Hence, the localization capability would only slightly worsen when applied to the 7 yr sample (for which no neutrino events are publicly available). However, even considering an extremely bad localization performance, e.g., assuming that the angular resolution of the source localization with the 7 yr neutrino sample is a factor of 2 larger than with the 10 yr sample, the binomial probability defined above would still be extremely small, ∼10⁻³. This makes the observed correlation inconsistent with a hypothetical signal and provides additional and independent evidence that the correlation claimed by Buson et al. (2022a, 2022b) could be due to the blazar sample used and statistical fluctuations.

6.2. RFC

6.3. The Blazar Contribution to IceCube Neutrinos

We have reproduced and extended the work of Buson et al. (2022a, 2022b), who tested the association of blazars in the 5BZCAT catalog with an IceCube neutrino map covering 7 yr of observations in the southern sky, to the northern sky and to a ∼60% larger and improved, recently released all-sky IceCube data set covering 10 yr. While we managed to reproduce the ∼4.5σ equivalent p-value of Buson et al. (2022a, 2022b), only insignificant results were obtained in the 7 yr northern sky and in the 10 yr data for both Northern and Southern Hemispheres. The original finding appears, then, to be only a statistical fluctuation, which might be partly due to the nature of 5BZCAT, which is not a flux-limited, well-defined catalog but just a (by now 9 yr old) list of sources detected in various bands.

We have also applied the same methodology to a complete, blazar-dominated AGN sample with 8 GHz VLBI flux density ≥150 mJy. In this case, we have derived a ∼3.3σ equivalent result in the 7 yr southern sky but, contrary to Plavin et al. (2021), who used a different approach but the same catalog, got an insignificant result in the 7 yr northern sky. However, as for 5BZCAT, no correlation was found in the 10 yr whole sky.

All of the above does not exclude a blazar contribution to the IceCube signal, which is supported by the neutrino association with TXS 0506+056, but is consistent with a scenario where only a small fraction of the γ-ray- and radio-selected blazars are high-energy neutrino sources, in agreement with previous works.

We thank the referee for the insightful and very detailed comments. C.B. thanks S. Coenders and M. Wolf for their help in developing parts of the analysis used in this work and M. Karl and H. Niederhausen for proofreading this manuscript. This work is supported by the Deutsche Forschungsgemeinschaft through grant SFB 1258 "Neutrinos and Dark Matter in Astro- and Particle Physics."

Facility: IceCube South Pole Neutrino Observatory - .

Software: SkyLLH ¹⁴ (Bellenghi et al. 2023), healpy (Zonca et al. 2019).

We search the entire sky for continuous neutrino emission using 10 yr of muon tracks recorded and processed by the IceCube Collaboration. We use an adaptation of the unbinned maximum likelihood method described by Braun et al. (2008) to search for neutrino point-like sources. We search for clusters of events in the sky while also using the information on their energy to further discriminate from random accumulations produced by the background. Hence, the statistical model is a mixture of a signal and a background model so that the data can be described by one of the two following hypotheses:

• 1.  
  H₀: the recorded events only consist of background muons, produced in the atmosphere or by the astrophysical diffuse neutrino flux;
• 2.  
  H₁: the recorded events are a mixture of background muons and of an additional component originating from an astrophysical point-like neutrino source, which emits with some strength n _s according to a power-law energy spectrum ∝E^−γ .

The test statistic is the negative logarithm of the likelihood ratio:

$$\begin{eqnarray}&&\mathrm{TS}=-2\mathrm{log}\displaystyle \frac{{ \mathcal L }({H}_{0}| \mathrm{Data})}{{ \mathcal L }({H}_{1}| \mathrm{Data})},\end{eqnarray} \tag{ A1 }$$

which is maximized with respect to the two free parameters in the signal model, n _s ≥ 0 and γ ∈ [1, 5].

The null hypothesis H₀ is a particular case of the alternative hypothesis H₁, i.e., H₁ reduces to H₀ when the point-like emission strength is zero (n _s = 0). Therefore, the two hypotheses are nested, and the likelihood of observing the data given the model is

$$\begin{eqnarray}&&{ \mathcal L }({n}_{{\rm{s}}},\gamma | {\boldsymbol{x}})=\displaystyle \prod _{i}^{N}\left[\displaystyle \frac{{n}_{{\rm{s}}}}{{\rm{N}}}{ \mathcal S }({x}_{i}| \gamma )+\left(1-\displaystyle \frac{{n}_{{\rm{s}}}}{{\rm{N}}}\right){ \mathcal B }({x}_{i})\right],\end{eqnarray} \tag{ A2 }$$

where the likelihood is unbinned, so the product is over all N events {x_i } in the sample, and ${ \mathcal S }$ and ${ \mathcal B }$ are the signal and background probability density functions (pdf's), respectively. Both depend on the reconstructed observables for the ith event, x_i = (E_μ,i , d _μ,i , σ_i ), where E_μ is the muon energy, d _μ is the muon direction, and σ is the uncertainty on the reconstructed direction. ¹⁵ The signal pdf additionally features a conditional dependence on the spectral index, γ, of the energy distribution. Each of the two pdf's can be separated into a spatial pdf and an energy pdf. The spatial and energy distributions of background events, as well as the spatial distribution of signal events (a 2D Gaussian centered at the assumed source location and with standard deviation equal to the σ of the event), are the same used in IceCube's typical point-source analysis (e.g., Carver 2019). The energy pdf in our signal likelihood is the only term that differs from the typical IceCube analysis. In fact, while IceCube analyses usually define a signal energy pdf conditional on the reconstructed muon decl., our signal energy pdf conditionally depends on the assumed source decl. ¹⁶ This modification is needed to cope with the absence of the reconstructed muon decl. in the detector response matrix, which is used to construct the signal energy pdf as

$$\begin{eqnarray}&&P({E}_{\mu }| {\delta }_{\nu },\gamma )=\int {{dE}}_{\nu }\,P({E}_{\mu }| {E}_{\nu },{\delta }_{\nu })P({E}_{\nu }| \gamma )P({E}_{\nu }| {\delta }_{\nu }).\end{eqnarray} \tag{ A3 }$$

The terms of the energy pdf in Equation (A3) are as follows:

• 1.  
  P(E_μ ∣E_ν , δ_ν ) is the probability of the muon energy E_μ given the parent neutrino energy E_ν and decl. δ_ν , and it is obtained from the detector response matrix;
• 2.  
  P(E_ν ∣γ) is the probability of the neutrino energy given a power-law energy flux with spectral index γ;
• 3.  
  P(E_ν ∣δ_ν ) is the probability of detecting a neutrino of energy E_ν from the direction δ_ν , calculated from the tabulated effective areas.

To compare our likelihood analysis to the one in Aartsen et al. (2020), we compare the respective sensitivities to an astrophysical E⁻² neutrino flux simulated at various declinations across the sky. ¹⁷ The events for the signal injection are sampled from the detector response matrix. The analysis used in this work is sensitive to an E⁻² astrophysical flux 20%–50% higher than the one that the analysis in Aartsen et al. (2020) would detect, depending on the decl., as shown in Figure 7.

Zoom In Zoom Out Reset image size

To calculate the 10 yr neutrino L (i.e., $-{\mathrm{log}}_{10}$(p-value) map), we divide the sky into a grid of equivalent pixels of ∼(011)² constructed using a HEALPix (Górski et al. 2005) projection with Nside parameter 512. At each of them, we maximize the TS (Equation (A1)) and calculate the L value as the negative logarithm of the probability of obtaining a TS larger than the observed one from background pseudo-experiments, generated by scrambling the experimental data in R.A. multiple times (e.g., Carver 2019). Thus, the L value at each pixel represents the probability that the data belong to the background-only hypothesis. The neutrino L sky map computed in this work is shown in Figure 8.

Zoom In Zoom Out Reset image size

We check the local significance at the locations of the four most significant sources contributing to the 3.3σ excess from a binomial test on a list of candidate neutrino sources published in Aartsen et al. (2020). We find L = 4.41 for NGC 1068, L = 2.70 for TXS 0506+056, L = 3.15 for PKS 1424+240, and L = 2.52 for GB6 J1542+6129. The respective L values published by the IceCube Collaboration are L = 4.74, L = 3.72, L = 2.80, and L = 2.74. Differences are within ∼0.2 Gaussian-equivalent standard deviations, except for the blazar TXS 0506+056, which has a local p-value an order of magnitude larger than the one reported by the IceCube Collaboration using the same data set (equivalent to a difference of ∼0.66σ). In the southern sky, we find a local significance of L = 4.83 for the published hottest spot, located at R.A. 3502 and decl. −565. This is slightly lower than the value reported by IceCube of L = 5.3 but still compatible within ∼0.2 standard deviations. The differences in the observed p-values, as well as in the sensitivity fluxes (Figure 7), can have various causes, from the difference in the energy term of the spatial likelihood to the coarseness of the information provided in the public detector response matrices. Indeed, it is worth recalling here that these matrices only have three bins in neutrino decl., covering the sky from −90° to −10°, from −10° to 10°, and from 10° to 90°. Especially at the horizon (where TXS 0506+056 is located), the detector response is averaged over a decl. range where the detector behavior changes significantly: the neutrino decl. bin includes the border between southern and northern sky, at −5°, where data processing and selection cuts change (e.g., Aartsen et al. 2020, and references therein).

The significance of the well-known blazar TXS 0506+056 in the 10 yr neutrino p-value map produced in this work is smaller than the one reported by Aartsen et al. (2020). Having L = 2.7, it does not pass the threshold ${L}_{\min }=3.0$ for being selected as an interesting hot spot for the correlation analysis. As a result, even though this source is included in both the 5BZCAT and RFC catalogs, the hot spot associated with it is not correlated with the two catalogs. As an a posteriori check, we test the impact of the absence of this hot spot on our analysis. We add to the list of our neutrino hot spots in the Northern Hemisphere an additional hot spot at the pixel corresponding to the location of the blazar and with the local significance reported in Aartsen et al. (2020), i.e., L = 3.7. We then run again the cross-correlation analysis in the northern sky with both the 5BZCAT and RFC catalogs. The results of this cross-check are shown in Figure 9. The best-fit ${L}_{\min }$ changes from 4.0 to 3.5 when using the 5BZCAT, while it stays at 3.5 for the RFC catalog. The pre-trial p-value of the correlation slightly decreases in both cases, going from 2.9% to 1.3% for the 5BZCAT and from 7.9% to 1.8%, but remaining fully compatible with the background expectation.

Zoom In Zoom Out Reset image size