A search for AGN sources of the IceCube diffuse neutrino flux

The origin of the diffuse astrophysical neutrino flux measured by the IceCube Observatory remains largely unknown. Although NGC 1068 and TXS 0506+056 have been identified as potential neutrino sources, the diffuse flux of neutrinos must have additional sources that have not yet been identified. Here we investigate potential correlations between IceCube's neutrino events and the Fermi and MOJAVE source catalogs, using the publicly-available IceCube data set. We perform three separate spatially-dependent, energy-dependent, and time-dependent searches, and find no statistically significant sources outside of NGC 1068. We find that, under the most optimistic assumptions of a spectral index of 2.0 and a neutrino flux uncorrelated with the gamma ray flux, no more than 13% of IceCube's neutrino flux originates from blazars over the whole sky. Then, using an energy-dependent likelihood analysis, the limit on neutrinos originating from blazars reduces to 9% in the Northern hemisphere under the same spectral index and flux assumptions. Finally, we set limits on individual sources from the MOJAVE radio catalog after finding no statistically significant time-flaring sources.


Introduction
The spectrum of IceCube's high-energy neutrino flux, first detected in 2013 [1], is consistent with a power-law extending from tens of TeV to a few PeV, and with flavor ratios consistent with a pion decay origin [2][3][4].While searches for point sources contributing to this neutrino flux have largely produced null results, there have been significant excesses from two specific sources.Evidence for the first source, TXS 0506+056, includes a neutrino event coincident in time and space with the flaring blazar detected by Fermi and MAGIC with a significance of 3 σ [5,6].More recently, NGC 1068, a nearby non-blazar active galactic nuclei (AGN), was identified as a potential neutrino source with a significance of 4.2 σ [7].However, these sources alone do not explain the majority of the diffuse neutrino flux, indicating that neutrinos in this energy range are produced by a large number of extragalactic sources [8].
There have been various proposals for possible sources of the diffuse neutrino flux.Gamma-ray bursts (GRBs) [9][10][11][12], star-forming galaxies [13], both blazar and non-blazar AGN [14][15][16], and fast radio bursts [17][18][19] have all been considered.However, various models have been ruled out as the primary source for IceCube's neutrino flux.In particular, there is an observed lack of correlation with blazar AGN, effectively eliminating flaring blazars as the primary source of the diffuse neutrino flux [20,21].Starburst and other star-forming galaxies are also unable to account for the entirety of this signal without exceeding the measured intensity of the isotropic gamma-ray background [22,23].
Our work in this paper consists of three analyses conducted using the publicly available data set from IceCube, containing 10 years of muon track neutrino events [24].First, we present an update of an earlier analysis that used IceCube's three-year public data set [8] to investigate whether a significant fraction of the neutrino flux could originate from blazar or non-blazar AGN.In this updated analysis, we compare IceCube neutrino events from the 10 year catalog to sources in the Fourth Catalog of AGN detected by the Fermi Large Area Telescope (the 4LAC catalog) [25], and our results are consistent with both the original work and an updated work [8,26].
We then go on to incorporate neutrino energy information, new in IceCube's 10 year catalog, and perform an energy-dependent analysis with a likelihood formulation that includes the provided "energy proxy" of the neutrino, focusing on the 4LAC sources in the Northern sky using the approach outlined in [27].We restrict this analysis to the Northern sky, defined conservatively as events with a declination of > 10 • , due to the coarse binning of the smearing matrices provided by the IceCube public data set.These smearing matrices describe the spread in reconstructed energy, direction, and angular uncertainty expected after developing Monte Carlo simulations of the IceCube detector.In the Southern sky, the lack of a dedicated cosmic ray simulation makes characterizing the muon background too challenging, and thus this region is omitted from this analysis.We also choose to omit the horizon region, defined by the smearing matrices as -10 • to 10 • , to be sure any edge effects in the detector response are not impacting our analysis.Thus, only the Northern sky above 10 • declination remains in the analysis.Even still, the addition of information about the neutrino energy distribution increases the sensitivity of the search overall; previous work found that with energy information included, the total number of events needed for a 5 σ significance is reduced by a factor of two [27].The same study also outlines a process that can be repeated in the Southern sky with a more robust energy reconstruction simulation that characterizes the muon background.
The third analysis presented here is a time-dependent threshold analysis of the MOJAVE XV radio catalog [28], in which, in addition to spatial coincidences, time windows are weighted corresponding to the magnitude of flares coming from a MOJAVE source at that time.A more complete description is included in Section 5.The fluctuations seen in the radio emission from these galaxies could be correlated with high-energy neutrino emission; p − γ interactions in the sources could create both neutrinos at high energies and photons at lower energies after the initial high-energy photon is lost through a chain of pair production [29].
Radio AGN in the MOJAVE radio catalog have previously been investigated as a possible source class for the IceCube neutrino flux in a time-independent analysis [30], and our addition of a time-dependent likelihood function probes a new signal hypothesis.The multiple epochs present in the MOJAVE data make it an ideal catalog for investigating a possible correlation between neutrinos and radio emission.We apply the process outlined in [31] for an individual source to the entire MOJAVE catalog.

Method Overview
Each analysis presented here uses data from IceCube's public data release consisting of muon track events from April 2008 to July 2018 [24].The data is from the 40 string detector in 2008, the 59 string detector in 2009, the 79 string detector in 2010, and the completed 86 string detector for the remaining 7 years.The entire data set contains 1,134,450 muon track events.Each event has a reported direction, angular resolution, time of event, and "energy proxy," which is related to the energy deposited in the detector.Each year of this data set includes an effective area for the detector as determined by simulation, as a function of declination and neutrino energy.To test for evidence of a neutrino signal from an individual point source, we follow the approach outlined in [32].The likelihood that a given source results in n s events, out of a total N recorded in the detector, is given by: where S i and B i are probability distribution functions (PDFs), ⃗ x s is the direction to the source, ⃗ x i is the reported direction of the event, and δ i is the declination of the event.The signal PDF consists of individual spatial, energy, and time PDFs: where P sig i , ϵ sig i and T sig i are the respective spatial, energy, and time components for the signal PDF, σ i is the uncertainty associated with each event, given in terms of angular resolution provided by the data set, E i is the energy proxy, γ is the spectral index (set to either 2 or 2.5 for each analysis), t i is the time of the event, and ϕ s is the flux of the source as a function of time.The background PDF is defined similarly by: where again P bkg i , ϵ bkg i , and T bkg i are the respective spatial, energy, and time components for the background PDF.For the work presented in this paper, the spatial component is included in all three analyses, while the energy and time components are only used in the energy-dependent analysis and the MOJAVE catalog study, respectively.In the following sections, the PDFs used in each analysis will be described.

Catalog Descriptions
We compare the measurements from IceCube to the results from two external catalogs: the Fermi 4LAC catalog [25] and the MOJAVE radio catalog [28].
The Fermi 4LAC catalog is comprised of active galactic nuclei (AGN) sources detected in the 50 MeV to 1 TeV energy range over an eight year period from August 2008 to August 2016.Of the 2863 total AGNs, 2796 are blazars, 63 are non-blazars, and 4 are unidentified.The blazars can be further broken down into 658 flat spectrum radio quasars (FSRQs), 1067 BL Lacs, and 1071 "blazars of unknown origin."We do not consider sources within 3 • of the poles.
The MOJAVE XV catalog consists of observations in the 15 GHz band of 437 AGN over the course of 20 years, between 1996 and 2016.This list of sources includes 265 quasars, 127 BL Lacs, 27 Radio Galaxies, 5 Seyfert-1 Galaxies, and 13 unidentified objects.During this time, 5321 measurements were taken of the AGN sources, with most sources visited between 5 and 15 times.These time-dependent measurements give unique insight into whether a source is flaring, either relative to its baseline or compared to similar sources.In this study, we only consider MOJAVE data taken after 2008, to align with the operational window of the IceCube instrument.
The function P B is equal to the fraction of events in the data set averaged across a band of ±6 • in declination, δ, around a given source.We only consider sources with declination between ±87 • due to the limited amount of solid angle near the poles with which to characterize the background PDF.To find the value of n s for a given source, the likelihood function is maximized with respect to the free parameter n s , as outlined in [27].
The statistical significance of a neutrino point source over a background-only hypothesis is calculated using the following test statistic: where L(0) is the likelihood for the background-only hypothesis in which there are no signal events from a given direction.From this, the p-value can be calculated by performing an integral over a χ 2 distribution with one degree of freedom.
Here, the background is being modeled empirically from the experimental data.This is possible due to the dominance of atmospheric muons and has been validated by scrambling the Right Ascension of the dataset and achieving statistically similar maps.The background hypothesis is that there is no correlation between the Fermi 4LAC catalog and the IceCube neutrino dataset; the signal hypothesis is the existence of a statistically-significant correlation between detected neutrino output and gamma-ray emission from known classes of astrophysical objects.
We additionally constrain the fraction of IceCube's diffuse neutrino flux that originates from known classes of astrophysical objects using three different weighting hypotheses for the expected neutrino emission from a given source class.This weighting is applied to n s in the likelihood, and is determined according to the following hypotheses: 1. Gamma-Ray Scaling: The neutrino flux from a given source is proportional to the gamma-ray flux from that source.The gamma-ray flux in the 4LAC catalog is in units of photons between 1-100 GeV per area, per time.This hypothesis would be expected if the gamma-ray emission is primarily produced from hadronic interactions, resulting in a fixed ratio of photons and neutrinos.
2. Geometric Scaling: The neutrino flux of a given source is proportional to the inverse of the luminosity distance squared, 1/D 2 L .This hypothesis assumes that the neutrino luminosity of a source is only correlated to distance between the source and the detector.This calculation can only be done with sources with measured redshifts, so some sources are excluded in the analysis under this hypothesis.

Flat Scaling:
The neutrino flux from a given source is uncorrelated to any other information.
These hypotheses remain unchanged from our previous two studies [8,23].

Results
To validate our methodology, in Fig. 1 we show an all-sky map of the likelihood of a neutrino point source, √ 2∆ ln L, in terms of right ascension and declination.The scan was performed in steps of 0.2 • , and at each point we show the value of n s that maximizes the test statistic defined in Eq. 3.2.In Fig. 2 we show the distribution of this test statistic across the sky.When excluding the points containing NGC 1068, a 4.2 σ source discovery by IceCube [7], the distribution is consistent with a normal distribution.In Table 1, we present the most significant locations in the sky other than NGC 1068, along with the associated likelihoods and pre-trial p-values.The post-trial p-value for the most significant source is calculated by scrambling the right ascension of the data to create 1000 randomly distributed backgroundlike distributions, and calculating how often these background-like distributions are more extreme than the real data.There is no evidence of a statistically significant source, which is the expected result, consistent with other all-sky scans of the IceCube dataset [33].
We then set limits on the fluxes from specific source classes, calculating the 95% upper limit as done in [8,23].When considered individually, no source has a statistically significant likelihood after accounting for the appropriate trials factor.When comparing the distances, luminosities, and spectra of the sources with the highest likelihood, we again find no clear features that set these sources apart from other sources in the 4LAC catalog more generally.
Finding no statistically significant individual sources, we place an upper limit on the neutrino flux from these source classes.In Fig. 3, the limits from this analysis on the contribution of blazar and non-blazar AGN to the IceCube diffuse neutrino flux are shown.We separately show results for sources identified as "non-variable" by applying cuts to sources with a variability index below 18.48, as reported by the Fermi Collaboration.Of the original 2796 blazars and 63 non-blazars, we identify 1674 and 47 "non-variable" sources, respectively.In all cases, we apply the appropriate completeness factors to account for the incompleteness of the source catalog [8].For blazars, we apply a completeness factor of 1.4.For non-blazars, we apply a completeness factor of 50.6 (154.6) for non-blazar (non-variable) AGN.We place 95% upper limits (corresponding 2∆ ln L = -3.84)for the three hypotheses outlined in Section 3, and for two choices of spectral index.We use a spectral index of α = 2.5, based on the spectral flux measured by IceCube, and α = 2.0, the value predicted assuming Fermi acceleration.
From the least constraining upper limits presented, we conclude that blazars can produce Figure 3.The 95% confidence level upper limits on the neutrino flux for various source classes using the spatial-only likelihood, presented as a function of neutrino energy.We compare these constraints to the diffuse neutrino flux reported by IceCube [34,35].In the left frames, we consider non-blazar AGN from the Fermi 4LAC catalog, and in the right frames we consider blazar AGN.In the bottom frames we plot only sources from the catalog that we identify as "non-variable" based on the flux density measurements over time, which corresponds to a 4LAC variability index of 18.48.The solid lines are for α = 2.0, while the dashed lines are for α = 2.5.The three different neutrino flux hypotheses (gamma-ray scaling, scaling by luminosity distance, and flat scaling), are shown in green, blue, and yellow respectively.The luminosity distance-scaled upper limits for α = 2.5 are outside the axes range.
no more than 13% of the diffuse neutrino flux, which is consistent with our previous analysis [8] and a publication earlier this year [36].This is a slightly more restrictive constraint than our previous result, due to the increased livetime of the 10-year IceCube data.Additionally, we find that under the flat scaling hypothesis, non-variable, non-blazar AGN could produce the entirety of the IceCube diffuse neutrino flux, and that in the case of the other two hypotheses, non-variable, non-blazar AGN could still produce a majority of the flux.
4 An Energy-Dependent Constraint on the Neutrino Flux from Blazar and Non-Blazar AGN In addition to the spatial-only likelihood analysis described above, we have also perfomed an analysis using an energy-dependent likelihood on data from the Northern hemisphere alone, which unlike the Southern hemisphere is not background-limited by atmospheric muons.We use the "energy proxy" for the muon track events in the IceCube data set, and adopt the smearing matrices from the IceCube data release and a spectral index of α = 2 to create a function that describes the likelihood of obtaining the reconstructed neutrino energy given the assumed spectral index, as outlined in [27].This signal likelihood is normalized to one and its distribution as a function of energy can be seen in the left panel of Fig. 4. To speed up the computational process, only events with an energy proxy greater than 100 GeV were included.The likelihood distribution of the Northern hemisphere data set as a function of the given "energy proxy" (in GeV) of the event.Right: The likelihood distribution in favor of a neutrino point source from the Northern hemisphere all-sky scan in the energy-dependent analysis.The observed distribution is consistent with a background-only hypothesis and we identify no evidence of another neutrino point source population.Sky locations with ∆ ln L < 0, corresponding to a best fit with a negative point source flux, are not shown.The error bars represent the 68% Poissonian confidence interval on each bin.
Most of the events in this data set are not astrophysical neutrinos; the Northern hemisphere data set is dominated by atmospheric neutrinos, which are a background for this analysis.Therefore, when constructing the energy-dependent background likelihood distribution, we use the data set itself, dependent on the declination of the point in the sky being tested.

Results
With the test statistic based on the new likelihood formulation in hand, we find the updated distribution of the test statistic shown in the right panel Fig. 4.This distribution falls along Figure 5.The 95% confidence level upper limits on neutrino flux for various source classes using the spatial and energy-dependent likelihood, presented as a function of neutrino energy.We compare these constraints to the diffuse neutrino flux reported by IceCube [34,35].In the left frames, we consider nonblazar AGN from the Fermi 4LAC catalog, and in the right frames we consider blazar AGN.In the bottom frames we plot only sources from the catalog that we identify as "non-variable" based on the flux density measurements over time, which corresponds to a 4LAC variability index of 18.48.The solid lines are for α = 2.0, while the dashed lines are for α = 2.5.The three different neutrino flux hypotheses (gamma-ray scaling, scaling by luminosity distance, and flat scaling), are shown in green, blue, and yellow respectively.Table 2.The five most significant independent locations found in our energy-dependent all-sky scan of only the Northern Hemisphere.The most significant source has a post-trial p-value of 0.456, which is not statistically significant.

2∆ ln L
a normal distribution, consistent with no significant source excesses.The pre-and post-trial p-values are calculated in the same was as in Section 3.1.The most significant locations on the sky in this search are presented in Table 2.While not statistically significant in this search, we do note that the brightest spot in our likelihood map, with a Right Ascension=182.60 and Declination=39.40, is within two arc minutes of NGC 4151, a source with recent evidence of neutrino emission [37].NGC 1068 falls just below the lower boundary of the analysis in declination and is thus excluded.We again compare against the Fermi 4LAC catalog with the same source hypotheses outlined in Section 3.Because we only analyze half the sky, the number of sources in our catalog is decreased by roughly 50%; however, we also expect the background to decrease by more than 50%, as the atmospheric events in the Southern hemisphere contribute significantly to the overall background rate at the highest energies.
The limit set by this energy-and spatial-dependent search in the Northern hemisphere with the 10-year data set are even more restrictive than that of the spatial-only analysis presented in Section 3. We use the same completeness factor as in the spatial-only analysis with an additional factor of 2 that accounts for only looking in the Northern hemisphere.
The resulting 95% upper limits are summarized in Fig. 5 for each of the neutrino flux scaling hypotheses.At most, non-variable blazar AGN could contribute up to 9% of the neutrino flux.The improved limits for blazar AGN again suggests that the majority of the IceCube neutrino flux is likely originating from a different source class.
The limits set for non-blazar AGN are also slightly improved after folding in the event energy.Even still, non-blazar, non-variable AGN could contribute up to 95% of the IceCube diffuse neutrino flux under the flat scaling hypothesis.This suggests that non-blazar AGN are still a viable candidate source class.

A Search for Correlations with MOJAVE Radio AGN
We apply a triggered time-dependent analysis on the MOJAVE catalog, which contains measurements of the radio emission from various AGN over a period of two decades, including one decade of overlap with the IceCube observatory.The time-dependent component of the likelihood T sig is created by calculating the mean value of the flux density for all sources in the MOJAVE catalog after being scaled appropriately with respect to their luminosity distance from Earth.Then, sources are identified as flaring if their flux density at a given time is greater than 2 σ above the average flux density.This is illustrated in Figure 6.
We perform the analysis using 436 radio-loud AGN from the MOJAVE catalog.Due to detector location and limits, the catalog only consists of sources higher than -30  Probability Density Function (PDF) Figure 6.The flux density, adjusted for luminosity distance, across the MOJAVE AGN catalog as a function of Modified Julian Dat (MJD).Also shown is the extrapolated PDF.For time periods with no flaring sources, defined as 2σ above the average recorded flux, the PDF is equal to 0.
lination.Each source in the catalog has a minimum of 5 flux density recordings at different times.Two sources are eliminated from the analysis, MOJAVE source 0415+379 and source 2200+420, as all of their flux density recordings were above the 2 σ threshold and can thus be categorized as always flaring.These two sources could be analyzed independently of this temporal analysis.This methodology was chosen to fairly assess all of the sources in the MOJAVE catalog, which can individually have as few as five data points and as many as 30 over the two decades of data.
Once the threshold on the catalog is set, any flux density point measured under the threshold is given a time likelihood value of 0. Above the threshold, we bin the data into 0.1 year bins, as outlined in [29] as the smallest reasonable bin size for Radio AGN.This was selected to ensure that there is minimal likelihood overlap, as larger time bins could result in neutrinos and flaring sources being erroneously correlated.After binning the flux density, the total flux density in all the bins is normalized to one.This creates the time-dependent portion of the signal PDF, where the likelihood for each muon track neutrino event is assigned based on the trigger time of the IceCube event.This likelihood is used in combination with the spatially-dependent likelihood described in Section 3. The delay between the arrival of the neutrino and the photons from the source is on the scale of seconds and can be ignored, as it is absorbed in the 0.1 year binning.
The time-dependent component of the background PDF is independent of seasonal vari- ations [31,38], but is dependent on the declination and year of the neutrino event.An additional time-dependence arises due to the construction of IceCube itself: more atmospheric muons were recorded in the first three years of this dataset than later years due to the phased construction of the detector.This prevalence of events with declination <0 • during these early years is described in the background PDF, with unique PDFs defined as a function of declination for each 0.1 year bin.After the first 3 years, the PDF no longer changes with time, only declination.

Results
In Fig. 7, we show the likelihood distribution of the MOJAVE sources tested.After accounting for trials factors, no statistically significant sources are identified.
We also consider each source individually to set independent limits on the expected flux of neutrinos from each source.These limits are shown in Fig. 8.The list of these sources with their reference name, coordinates, and upper limit neutrino flux density can be found in Appendix A. Because we only consider individual sources, rather than source classes, no completeness factor is required.The MOJAVE catalog includes sources in four different subclasses: quasars, non-blazar AGN, blazar AGN, and Seyfert galaxies.Across each of these subclasses, there is no obvious relationship or trend between source type and neutrino flux limit.
This time-dependent method could be further refined in future work by varying the time window used to define a flaring source or by combining it with the energy-dependent method discussed in Section 4.

Conclusion
The analyses presented here utilize spatial, energy-dependent, and time-dependent likelihood functions to compare the IceCube diffuse neutrino flux to the Fermi 4LAC and MOJAVE catalogs.We confirm that blazar AGN are unlikely to explain the entirety of the IceCube neutrino flux, while non-blazar AGN are not ruled out.Adding an energy dependence to the previously described spatial analysis tightened constraints on the flux fraction for all source classes tested.The time-dependent analysis with the MOJAVE catalog sets neutrino flux limits for flaring sources broken down by source class.
While we are unable to definitively confirm a specific origin for the diffuse neutrino flux, we are hopeful that future analyses will utilize increased livetime, more complete source catalogs, and a substantive cosmic ray simulation to improve these results further.Additionally, considering specific theoretical models for source neutrino production as a function of neutrino energy could be included to further constrain possible sources.The time-dependent analysis framework could be used in tandem with other time-dependent catalogs as they become available.9.The 95% confidence level upper limits on the neutrino flux coming from the 175 MOJAVE sources that had an n s > 0 and √ 2δ ln L > 0. On the x axis is the source number and on the y axis is the neutrino flux from the source at the detector.The sources are presented in descending flux order and sorted by source type.The entirety of this list with the flux upper limits, reference name, and source class is in Appendix A.

Figure 1 .
Figure 1.An all-sky map of the likelihood distribution of neutrino point sources, √ 2∆ ln L, in Right Ascension (horizontal axis) and Declination (vertical axis) in an Aitoff projection.

Figure 2 .
Figure2.The likelihood distribution in favor of a neutrino point source from bins of our all-sky scan.The observed distribution is consistent with a background-only hypothesis when not including NGC 1068 in the best-fit function (dark blue line) and we identify no evidence of another neutrino point source population.Sky locations with ∆ ln L < 0, corresponding to a best fit with a negative point source flux, are not shown.The error bars represent the 68% Poissonian confidence interval on each bin.

Figure 4 .
Figure 4. Left: The likelihood distribution of the Northern hemisphere data set as a function of the given "energy proxy" (in GeV) of the event.Right: The likelihood distribution in favor of a neutrino point source from the Northern hemisphere all-sky scan in the energy-dependent analysis.The observed distribution is consistent with a background-only hypothesis and we identify no evidence of another neutrino point source population.Sky locations with ∆ ln L < 0, corresponding to a best fit with a negative point source flux, are not shown.The error bars represent the 68% Poissonian confidence interval on each bin.

Figure 7 .
Figure 7.The likelihood distribution in favor of a neutrino point source at the locations of the MOJAVE sources.The observed distribution is consistent with a background-only hypothesis, with a slight dip in the 1-2 √ 2∆ ln L range, and we identify no evidence of another neutrino point source population.Sky locations with ∆ ln L < 0, corresponding to a best fit with a negative point source flux, are not shown.The error bars represent the 68% Poissonian confidence interval on each bin.

Figure 8 .
Figure 8.The 95% confidence level upper limits on the neutrino flux coming from the 175 MOJAVE sources that had an n s > 0 and √ 2δ ln L > 0, as a function of declination.The shape of the limits match the declination-dependent effective volume of the IceCube instrument.The entirety of this list with the flux upper limits, reference name, and source class is in Appendix A.

Table 1 .
The five most significant independent locations found in our spatial-only all-sky scan, excluding the coordinates of NGC 1068.The most significant source has a post-trial p-value of 0.451, which is not statistically significant.
Pre-Trials p-value RA