Search for Correlations of High-energy Neutrinos Detected in IceCube with Radio-bright AGN and Gamma-Ray Emission from Blazars

The IceCube Neutrino Observatory sends realtime neutrino alerts with a high probability of being astrophysical in origin. We present a new method to correlate these events and possible candidate sources using 2089 blazars from the Fermi-LAT 4LAC-DR2 catalog and with 3413 active galactic nuclei ( AGNs ) from the Radio Fundamental Catalog. No statistically signi ﬁ cant neutrino emission was found in any of the catalog searches. The result suggests that a small fraction, < 1%, of the studied AGNs emit neutrinos that pass the alert criteria, and is compatible with prior evidence for neutrino emission presented by IceCube and other authors from sources such as TXS 0506 + 056 and PKS 1502 + 106. We also present cross-checks to other analyses that claim a signi ﬁ cant correlation using similar data samples.


INTRODUCTION
Astrophysical neutrinos are weakly-interacting, electrically neutral particles capable of revealing the origin of cosmic rays.Thought to be produced in regions near some of the most high-energy events in the Universe, astrophysical neutrinos travel undeflected by magnetic fields delivering information relating to their original production mechanism.One example of these regions is NGC 1068, an Active Galactic Nucleus (AGN), from which evidence of high-energy neutrino emission was recently found with the IceCube Neutrino Observatory (Abbasi et al. 2022a).
Despite ongoing research, the origins of astrophysical neutrinos remain largely unknown.However, AGNs are among the most promising candidate sources (Abbasi et al. 2022a,b;Murase & Stecker 2022).In particular, blazars -AGNs that have jets pointing towards Earth-are thought to be likely sources.Recently, evidence has emerged of a spatial and temporal correlation with the emission from blazar TXS 0506+056 and the high-energy neutrino event IC170922A, detected by IceCube (Aartsen et al. 2018b).An alert was sent to other observatories and a follow-up campaign in multiple wavelengths was started.At the moment of neutrino detection, the blazar was flaring in gamma rays, and it was determined that the chance of accidental coincidence was at the 3σ level.Subsequent retrospective analyses from the IceCube collaboration found an excess in neutrino events with respect to the expected atmospheric background coincident with the position of TXS 0506+056, with a significance at the 3.5σ level (Aartsen et al. 2018a).Later studies have also indicated a possible correlation between the blazar PKS 1502+106 and the IceCube neutrino event IC190730A (Rodrigues et al. 2021).† also at Institute of Physics, Sachivalaya Marg, Sainik School Post, Bhubaneswar 751005, India ‡ also at Department of Space, Earth and Environment, Chalmers University of Technology, 412 96 Gothenburg, Sweden § also at Earthquake Research Institute, University of Tokyo, Bunkyo, Tokyo 113-0032, Japan Both IC170922A and IC190730A were high-energy neutrino events with a high probability of being astrophysical in origin.For each of these events, IceCube released public realtime "neutrino alerts".For an alert to be issued, the event must have passed at least one of several filtering streams based on its reconstructed energy and direction (Aartsen et al. 2017b;Blaufuss et al. 2019).These alerts can be used to search for correlations between different candidate sources, as other studies have shown, e.g. a search for correlations with radio-selected AGNs in Plavin et al. (2020) and with blazars in Giommi et al. (2020).
In this paper we develop a method to perform correlation analyses using the latest catalog of alerts that IceCube has published, IceCat-1 (Abbasi et al. 2023) 1 , and various catalogs of candidate sources.These analyses also use different weighting schemes for the sources to account for different neutrino production mechanisms.As a use case, we present three analyses: two using blazars from Fermi 4LAC-DR2 catalog and one with AGNs from the RFC 2022a catalog.The paper is organized as follows: in Section 2 the IceCube Neutrino Observatory is introduced, in Section 3 the description of the neutrino catalog and the two source catalogs considered for this paper are discussed and in Section 4 we introduce the analysis technique.In Section 5 the results are shown.Section 6 contains the conclusion and a comparison to previous analyses that used similar samples.

THE ICECUBE NEUTRINO OBSERVATORY
The IceCube Neutrino Observatory is a cubic-kilometer detector located at the South Pole (Aartsen et al. 2017a).It is a three-dimensional array of 5, 160 Digital Optical Modules (DOMs) deployed along strings between depths of 1, 450 m and 2, 450 m below the ice.Each DOM contains a photomultiplier tube (PMT) (Abbasi et al. 2010) and readout electronics (Abbasi et al. 2009).The array of DOMs detects neutrinos indirectly via the Cherenkov radiation produced by the relativistic charged particles that result from the interactions of neutrinos in the ice or bedrock below IceCube.The analyses presented in this paper are focused on track-like events, produced by charged-current muon-neutrino interactions.

IceCube's neutrino alert catalog
The IceCube Realtime System was established in 2016 (Aartsen et al. 2017b).Recently, the selection criteria were updated (Blaufuss et al. 2019) and now consists of two alert streams, "Gold" and "Bronze" neutrino alerts, which are defined by an average astrophysical purity (signalness) of 50% and 30% respectively.The signalness is calculated as the ratio between the expected number of astrophysical signal events and the expected total number of signal and background events with energy larger than the energy of the event and at the given declination.The expected number of signal events is calculated assuming a spectral index of 2.19 ± 0.10 (Haack & Wiebusch 2017) for the signal neutrino flux.The latest diffuse astrophysical muon neutrino spectrum measured by IceCube has an index of 2.37±0.09(Abbasi et al. 2022c), but the signalness definition is kept unchanged for consistency.Typically, alert events have reconstructed neutrino energies above 100 TeV.
When IceCube detects a high-energy neutrino that passes the alert criteria, an automatic GCN2 Notice is sent out to other observatories to allow for follow-up observations in multiple wavelength bands.This first notice contains the reconstructed direction with statistical uncertainties calculated with a fast reconstruction algorithm (Abbasi et al. 2021a).Another, more sophisticated reconstruction algorithm is then run and a second, updated GCN Notice and a GCN Circular are sent out a few hours later, including the refined reconstructed direction and uncertainties that also include systematic uncertainties.
This second algorithm is a maximum likelihood reconstruction method (Aartsen et al. 2014) that determines the neutrino direction by running the track reconstruction for each direction (θ, φ) on a HEALPix (Gorski et al. 2005) grid covering the sky, with iteratively finer binning.The likelihood is calculated at every point of the grid, and the direction maximizing the likelihood represents the track whose light yield matches the observed data and therefore is used as the reconstructed direction of the neutrino.This likelihood map can be normalized to act as a spatial probability density function for the neutrino (Figure 1).From the map one can also derive uncertainty contours following the procedure explained in Abbasi et al. (2021b), in which an event is simulated and reconstructed multiple times to determine the appropriate confidence level and to account for systematic uncertainties.In practice, all the events published in the realtime alert system in IceCube are calibrated using the simulations of one archival event, IC160427A.The work presented in Abbasi et al. (2021b) indicates that this method has difficulties in correctly accounting for systematic uncertainties and leads to a possible overcoverage of the angular resolution.

Equatorial
Run: 137096 Event: 70551815 Type: HESE MJD: 59850.522665511824 Online scan best-fit (50%, 90%) Offline scan best-fit (50%, 90%) 90% bounding rectangle 0 50 100 150 200 250 2 ln(L) Figure 1.Example of a log-likelihood map of a neutrino alert reconstructed with the full-sky maximum likelihood reconstruction algorithm in right ascension and declination.The black star is the best-fit position from the likelihood scan.The solid and dashed black curves represent the 50% and 90% confidence level respectively, derived from the likelihood scan.The red dashed rectangle determines the errors quoted in the GCN Circulars and is defined as the minimum rectangle that would encapsulate the uncertainty contours.Lastly, the purple solid and dashed curves are the 50% and 90% confidence level contours calculated with the initial, less refined reconstruction method and sent out in the first GCN Notice.
The IceCube Event Catalog of Alert Tracks-1 (IceCat-1) (Abbasi et al. 2023) includes all the alerts that were issued with the revised criteria as well as archival data starting from 2011, prior to the start of the realtime stream, which would have been selected as realtime alerts and were processed in the same way.In total, it contains 275 track-like events.Abbasi et al. (2023) also shows a simple correlation analysis searching for spatial coincidences between neutrino alerts and different candidate sources.The correlation analysis presented here represents an updated method with increased sensitivity to potential signal, since it includes additional information of the alert event and the candidate counterpart.In Figure 2, the spatial distribution of the alerts is shown, with the colours representing the signalness of each event.

Fermi-LAT sources
For the first two analyses, blazars were selected from the incremental version of the Fourth Catalog of Active Galactic Nuclei detected by Fermi -LAT (4LAC-DR2, for Data Release 2) (Ajello et al. 2020;Lott et al. 2020).This catalog contains all AGNs detected by Fermi -LAT between 2008 and 2018 at high galactic latitudes b (|b| > 10 deg).We select blazars (BL-Lacs, Flat Spectrum Radio Quasars and blazar candidates of unknown types), which represent 98% of the catalog.
The catalog was further refined by studying the distribution of sources in the sky to assure that it is isotropic, which fulfills the condition of isotropy that goes into the definition of the Test Statistic, explained in Section 4. Various selection criteria with increasing limits on the brightness were applied to the gamma-ray energy flux F , which is defined as the energy flux of photons of a certain energy per unit area, F = (dN/dEdAdt)E 2 , and obtained by spectral fitting from 100 MeV to 100 GeV.The latitude and longitude distributions of the blazar data after the selection were compared to a uniform distribution, both weighted by the energy flux and unweighted.A cut of F ≥ 10 −11.6 erg cm −2 s −1 was chosen based on the results retaining 90% of the cumulative flux from all objects in the catalog.The final number of blazars used in the analysis is 2, 089 (from the original 2, 863 objects in 4LAC-DR2).
Two weighting schemes are used in analyzing the 4LAC-DR2 catalog.For the first, the normalized 10-year average energy flux in gamma rays above 100 MeV of each blazar is used as the weight, since the blazars with higher flux are expected to emit more high-energy neutrinos (Murase et al. 2014;Petropoulou et al. 2015) and, therefore, are more . Neutrino events which have passed the alert criteria in the period between 2011 and 2020, in equatorial coordinates.
The colours represent the probability that the neutrino is of astrophysical origin (signalness).The varying density in the number of events reflects the dependence of IceCube's effective area on the declination of the neutrino direction.
likely to be the source of an alert neutrino.A similar weighting scheme was followed in a previous IceCube analysis using a different neutrino catalog (Huber 2019).The second weighting scheme utilizes the light curves from the Light Curve Repository3 made available by the Fermi -LAT collaboration (Fermi Large Area Telescope Collaboration 2021).
Here we use the normalized flux in the 1-month time bin that coincides with the neutrino arrival time as the weight, and thereby use the blazar's temporal information.This procedure is similar to the one followed in Aartsen et al. (2018b).From the 2, 089 blazars selected for this analysis, 1, 161 had flux variability and were released in the light curve repository.The weights are calculated as follows: • If the light curve of a blazar is not provided in the repository, it is assumed that it is not variable and the average energy flux is considered as the flux at the time of the neutrino's arrival.
• If there is a flux measurement at the neutrino arrival time, it is used as the weight.
• If there is no flux measurement in the time bin of the neutrino arrival or it is an upper limit, we check the monthly time bin before and after the time bin of the neutrino event (neighbouring bins).If both have a flux measurement, the weight is an average of both.If only one of them has a measurement, the weight is equal to that flux.This assumes that the monthly energy flux does not vary greatly from one month to another.
• If there is no flux measurement or only an upper limit is provided in the time bin of the neutrino event and its neighbouring time bins, then a weight of zero is assigned.The weight is set to zero also when there are no measurements available in the neighbouring bins.This is to avoid mixing the hypotheses for the two weighting schemes.
For the analysis presented here, we use the flux density integrated over VLBI images at 8.6 GHz as a weight and select AGNs with flux density S 8.6 GHz ≥ 150 mJy, which leaves a total of 3, 413 AGNs (828 in common with the subsample from 4LAC-DR2 described in the previous section).We follow the selection criteria applied in Plavin et al. (2020) with the purpose of cross-checking the results obtained in their analysis.The procedure followed to compare both methods is explained in Appendix A and discussed in Section 6.

ANALYSIS TECHNIQUE
The statistical significance of the overall correlation between the neutrino alerts and various candidate source catalogs is calculated using the stacking analysis method.In this method, the hypothesis that we test is that the neutrino is correlated with a coincident source from the candidate source class."Coincident" is defined here as the source lying within the 90% error contour of the reconstructed arrival direction of the neutrino.The null hypothesis is that the neutrino is not correlated with any sources, either because it is atmospheric background or because it comes from a source not included in the catalog being tested.The likelihood L is the weighted sum of a signal Probability Density Function (PDF), S, and a background PDF, B, where n S is the number of signal events and N is the total number of signal and background neutrino candidate events.The Test Statistic (TS) for a single neutrino i (N = 1) is defined as Treating each neutrino independently, the hypothesis is tested for each coincident source b.The combination returning the highest TS is selected as the most likely counterpart.If TS < 0, the background hypothesis is chosen, i.e.S/B = 1.The signal PDF is defined as where S signalness,ν is the signalness of the neutrino alert; S spatial,ν ( x b ) is the neutrino point-spread function, calculated using the re-scaled likelihood maps from the reconstruction of the alerts that depends on the position of the source x b (see Section 3); and w b is the source weight that depends on the hypothesis of neutrino production mechanism.The most likely counterpart is then the source for which S spatial,ν ( x b ) • w b is maximum, and we denote that in Equation 2 as Ŝ.The signalness is a fixed parameter for each individual neutrino, so it does not affect the selection of the counterpart.
For the background case, in which the neutrino is not associated with any source, the probability of finding a correlation by chance can be assumed to be isotropic, and therefore the background PDF is simply defined by B = B spatial = 1 4π .Finally, the TS of the data is calculated by summing over all neutrino alerts, TS data = i TS i .To calculate the significance of the result, N = 5000 background pseudo-experiments are constructed by fixing the direction of the neutrino alerts and scrambling the candidate sources in right ascension and declination, while maintaining the high galactic latitude criterion from the Fermi 4LAC-DR2 catalog of |b| > 10 deg.Then the TS is calculated for each of the background maps and its distribution is fitted to a gamma distribution for the p-value estimation.
For the analysis using Fermi -LAT 4LAC catalog with average flux weights, the sensitivity flux, which is defined as neutrino flux needed in order to get a TS larger than the median of the background distribution in 90% of the cases, is at 5.1% of the astrophysical flux, which corresponds to 6.3 signal correlations from the 275 neutrinos in the IceCat-1 catalog, i.e. 6.3 neutrinos correlated with their sources.The sensitivity for the same catalog with the monthly energy flux weighting scheme (second analysis) is at 4.4% of the astrophysical flux, which is equivalent to 5.5 signal correlations.Using the RFC catalog and the average VLBI the sensitivity flux is 6% of the astrophysical flux, or 7.4 signal correlations.

Fermi 4LAC-DR2 with average energy flux
No significant correlation was found between the 2, 089 selected blazars from Fermi 4LAC-DR2 weighted by the average 10-year energy flux and the neutrino alerts, with a p-value of 0.248, see Figure 3 Panel a.The 10 most significant correlations (the 10 blazar-neutrino pairs with the highest TS i in the data) are listed in Table 1.The two strongest correlations have been extensively studied, (see Section 1, TXS 0506+056 -IC170922A and PKS 1502+106 -IC190730A).The third and fourth correlation are also interesting since 3C 454.3 and Mkn 421 are two of the brightest sources in the sky, which explains the high TS, but in both cases the area of the error region of the associated neutrino is very large.The areas of the 90% uncertainty contours are 26.80 deg 2 and 56.07 deg 2 respectively, compared to the median 90% contour area of all alerts which is 5.96 deg 2 , or the 90% contour area of IC170922A and IC190730A, of 1.54 deg 2 and 5.61 deg 2 , respectively.Table 1.The table shows the 10 correlations that contribute most to the TS of the data using the Fermi 4LAC-DR2 catalog and average flux weights.The first column is the neutrino using the nomenclature of the realtime program, the second column is the signalness of each neutrino, the third column is the J2000 name of the most likely counterpart, the fourth column is the angular distance between the candidate source and the neutrino (calculated as the angle between the two sightlines) and the last column is the contribution of each correlation to the final TS.

Fermi 4LAC-DR2 with monthly energy flux
The TS is again consistent with the null hypothesis when the selected blazars are weighted by the energy flux at the neutrino arrival time (see Figure 3 Panel b).The p-value of this analysis is 0.08.Table 2 shows the 10 most contributing correlations to the TS.The main contribution to TS data is given by TXS 0506+056 and IC170922A, since the blazar was flaring in gamma rays at the neutrino arrival time as mention in Section 1, leading to a higher TS value than observed in the analysis with average energy flux weights.On the other hand, one can see that the blazar 3C 454.3 was in a quiescent state at the neutrino arrival time, and therefore, the TS value of that correlation is low and is not within the 10 most contributing ones.

RFC 2022 with average VLBI flux
The result from this third analysis, shown in Figure 3 Panel c, using the RFC catalog and the average VLBI flux at 8.6 GHz as weight for the contribution of the AGNs is that the data is compatible with the null hypothesis, with a p-value of 0.57.In Table 3 the 10 correlations with the highest TS are shown.It is worth noting that the correlations of TXS 0506+056 with IC170922A and PKS 1502+106 with IC190730A appear on the list (with the nomenclature used in the RFC catalog for the AGNs), since they are both very bright in radio as well.
Table 2.The table shows the 10 correlations that contribute most to the TS of the data using the Fermi 4LAC-DR2 catalog and monthly flux weights.The first column is the neutrino using the nomenclature of the realtime program, the second column is the signalness of each neutrino, the third column is the J2000 name of the most likely counterpart, the fourth column is the angular distance between the candidate source and the neutrino (calculated as the angle between the two sightlines) and the last column is the contribution of each correlation to the final TS.In this paper we have described a new method developed within IceCube that allows us to calculate the overall correlation between neutrino alerts and a catalog of candidate sources.We conducted three analyses: two searching for correlations with blazars from Fermi 4LAC-DR2 catalog and one with AGNs from the RFC 2022a catalog.A significant correlation was not found in any of the three cases, meaning that the majority of the neutrino alerts are neither produced in bright or flaring gamma-ray blazars nor in radio-bright AGNs.

Neutrino
The analyses presented here, and the method itself, rely heavily on the robustness of the reconstruction.Due to the way the contours are calculated, which determines how broad the neutrino spatial PDF is, a better handle on the systematic uncertainties is crucial.The reconstruction method used in IceCube's realtime system is being improved and further modifications will be implemented in the future (Abbasi et al. 2021b), which might give an opportunity to revisit this analysis with improved sensitivity.
The results obtained in this paper are, however, compatible with previous hints of a signal, e.g.Aartsen et al. (2018b); Rodrigues et al. (2021).The expected number of chance correlations that have T S i > 7 is 0.077 for the 4LAC-DR2 catalog weighted by the average energy flux.In data, there are three such correlations, which point to potential signal being present but non discernible from chance coincidences.The observed T S data is compatible with the TS distribution that one would obtain if only 1% of the blazars in the selected Fermi 4LAC-DR2 catalog emits neutrinos.This is calculated by randomly choosing 1% of sources in multiple pseudo-experiments and injecting them For comparison, our definition of the TS is kept, but the spatial term in the signal PDF is calculated as in Plavin et al. (2020).The information that IceCube publishes in the GCN Circular whenever an alert is detected is the reconstructed direction along with the minimum rectangle that encapsulates the 90% error contour in right ascension and declination (see Figure 1).In Plavin et al. (2020) the coordinate-wise errors are multiplied by the ratio of the 90% quantiles of 2D and 1D Gaussian distributions (≈1.3).An extra parameter x is added linearly in all directions to account for possible unknown systematic uncertainties.It varies from 0 to 1 deg in increments of 0.1 degrees and fitted to find the minimum pre-trial p-value, which is obtained at x = 0.5 deg in Plavin et al. (2020).Finally, the error region is assumed to be composed of 4 points joint by elliptical segments.This is shown in Figure 4. IceCube already calculates the errors by assuming a 2D distribution, but the extra 1.3 factor is also applied in the comparison for consistency with Plavin et al. (2020).Systematic uncertainties are added to the contours before publishing them, although one could argue that there are unknown systematics.This is currently being studied, see e.g.Abbasi et al. (2021b), but the method to calibrate the errors is conservative and the real uncertainties could be smaller.However, it is worth investigating this scaling up of the contours since other analyses have also found interesting hints of signal following a similar procedure (although in some cases with a very different neutrino data sample focused on lower energies), e.g., Giommi et al. (2020); Buson et al. (2022).
For the comparison analysis we switch the spatial term of the signal PDF from the neutrino point-spread function (based on the reconstruction's likelihood landscape) to the simpler coincidence-based function described above, multiplying the angular errors by 1.3 and considering the added 0.5 deg as a fixed value.Then the S spatial is 1 for all AGNs inside the area of the calculated contour and 0 outside.This is performed with the datasets described in Section 3 (IceCat-1 catalog + RFC 2022a) and the original dataset from Plavin et al. (2020) (56 alert and alert-like events + RFC 2019c) with and without the scaling factors, 1.3 from the quantiles + 0.5 degrees for systematic uncertainties.Without scaling, the error bars would be similar to the ones in Figure 4 Panel b), while for the scaled-up contours the error bars would be as in Figure 4 Panel d).For the neutrinos from the 56 alert and alert-like event catalog that do not have a signalness, a signalness of 45.1% is assumed, which is the average of the IceCat-1 catalog.
The results can be seen in Table 4.The lowest p-value is achieved with the test that is the most similar to Plavin et al. (2020) (same AGN and neutrino catalogs and scaling of the angular uncertainty).However, when more events are added, the significance drops.The same happens when the angular errors are not scaled up, which indicates that the correlations which are contributing to the TS are from AGNs that are from outside the original contours, and therefore far away from the best-fit position.Only one of the four correlations that drive the result in Plavin et al. (2020) appears in Table 3, and in the two tests with the RFC 2019c + alert-like events selection, due to its high signalness and the brightness of the AGN.The four most significant correlations from Plavin et al. (2020) are found in the test of RFC 2019c + alert-like events including scaling as expected, since this is the test that is closest to their analysis.With the information available, we cannot confirm the findings of Plavin et al. (2020), and the only hint of excess disappears when we use a more sophisticated description of the spatial probability density function for the neutrino events.

Figure 3 .
Figure 3.The blue histogram shows the distribution of the TS values for N = 5000 background maps, fitted to a gamma distribution (blue curve).The orange dotted line represents the observed TS of the data, while the grey lines represent the TS data needed to reject the background hypothesis at the 3σ or 5σ level.Panel a) For the Fermi 4LAC-DR2 catalog with the average energy flux weights, TS data = 131.3,while TS = 120.4for the median background maps, leading to a p-value of 0.248.Panel b) For the Fermi 4LAC-DR2 catalog with the monthly energy flux weights, TS data = 104.6,median background maps has TS = 83.6 and the p-value is 0.08.Panel c) For the RFC 2022a catalog, TS data = 155.8,TS = 158.9for the median background maps and the p-value is 0.57.

Figure 4 .
Figure 4. Example of the simplification and scaling that is applied to the realtime contours in Plavin et al. (2020), which is also done in the cross check presented here.From the original contours in Panel a), only the minimum rectangle is published, which is translated to uncertainty in the coordinates (Panel b).Those errors are multiplied by 1.3 (Panel c) and an extra parameter is added to account for systematic uncertainties, leading to a minimum p-value at an added 0.5 degrees (Panel d).The error bars are joined by elliptical segments to obtain the simplified contours corresponding to the best-fit position from the reconstruction, represented by the black star (Panel e).

Table 3 .
The table shows the 10 correlations that contribute most to the TS of the data using the RFC 2022a catalog.The first column is the neutrino using the nomenclature of the realtime program, the second column is the signalness of each neutrino, the third column is the J2000 name of the most likely AGN counterpart, the fourth column is the angular distance between the candidate source and the neutrino and the last column is the contribution of each correlation to the final TS.