Search for neutrino emission from cores of active galactic nuclei

The sources of the majority of the high-energy astrophysical neutrinos observed with the IceCube neutrino telescope at the South Pole are unknown. So far, only a flaring gamma-ray blazar was compellingly associated with the emission of high-energy neutrinos. However, several studies suggest that the neutrino emission from the gamma-ray blazar population only accounts for a small fraction of the total astrophysical neutrino flux. In this work we probe the production of high-energy neutrinos in the cores of active galactic nuclei (AGN), induced by accelerated cosmic rays in the accretion disk region. We present a likelihood analysis based on eight years of IceCube data, searching for a cumulative neutrino signal from three AGN samples created for this work. The neutrino emission is assumed to be proportional to the accretion disk luminosity estimated from the soft x-ray flux. Next to the observed soft x-ray flux, the objects for the three samples have been selected based on their radio emission and infrared color properties. For the largest sample in this search, an excess of high-energy neutrino events with respect to an isotropic background of atmospheric and astrophysical neutrinos is found, corresponding to a post-trial significance of 2 . 60 σ . If interpreted as a genuine signal with the assumptions of a proportionality of x-ray and neutrino fluxes and a model for the subthreshold flux distribution, then this observation implies that at 100 TeV, 27% – 100% of the observed neutrinos arise from particle acceleration in the core of AGN at 1 σ confidence interval.


1
The sources of the majority of the high-energy astrophysical neutrinos observed with the IceCube neutrino telescope at the South Pole are unknown. So far, only a flaring gamma-ray blazar was compellingly associated with the emission of high-energy neutrinos. However, several studies suggest that the neutrino emission from the gamma-ray blazar population only accounts for a small fraction of the total astrophysical neutrino flux. In this work we probe the production of high-energy neutrinos in the cores of active galactic nuclei (AGN), induced by accelerated cosmic rays in the accretion disk region. We present a likelihood analysis based on eight years of IceCube data, searching for a cumulative neutrino signal from three AGN samples created for this work. The neutrino emission is assumed to be proportional to the accretion disk luminosity estimated from the soft x-ray flux. Next to the observed soft x-ray flux, the objects for the three samples have been selected based on their radio emission and infrared color properties. For the largest sample in this search, an excess of high-energy neutrino events with respect to an isotropic background of atmospheric and astrophysical neutrinos is found, corresponding to a post-trial significance of 2.60σ. If interpreted as a genuine signal with the assumptions of a proportionality of x-ray and neutrino fluxes and a model for the subthreshold flux distribution, then this observation implies that at 100 TeV, 27%-100% of the observed neutrinos arise from particle acceleration in the core of AGN at 1σ confidence interval. DOI: 10.1103/PhysRevD.106.022005

I. INTRODUCTION
IceCube is a cubic-kilometer neutrino telescope operating at the South Pole, detecting neutrinos of all flavors with energies from tens of GeV to several PeV, based on the Cherenkov radiation of charged secondaries produced in neutrino interactions in the ice and bedrock [1]. In 2013, IceCube discovered high-energy astrophysical neutrinos [2] and has more recently found compelling evidence for the gamma-ray blazar TXS 0506 þ 056 during its flaring state being a source of high-energy neutrinos [3,4]. However, the sources responsible for the emission of the majority of the detected astrophysical neutrinos are still unknown. Furthermore, several studies imply that the population of gamma-ray blazars detected by the Fermi Large Area Telescope (LAT) can only account for a small fraction of the observed astrophysical neutrinos (e.g., [5][6][7][8][9][10][11][12]).
Active galactic nuclei (AGN) have been considered promising potential sites for high-energy neutrino production, as they are among the most powerful emitters of radiation in the known Universe [13][14][15][16][17][18][19][20][21][22][23][24]. They have the potential to accelerate protons up to the highest observed cosmic-ray (CR) energies and are surrounded by highintensity radiation fields where photo-nuclear reactions with subsequent neutrino production can occur. A common scenario, especially in blazars, assumes neutrino production in the relativistic jet. However, in this case, the production of neutrinos is boosted due to the relativistic motion of the jet, and therefore dominated by the blazar subclass with their jets closely aligned towards our line of sight.
Alternatively, CRs can be accelerated in the core regions of AGN, inside the accretion disk or the corona surrounding it [17,[25][26][27][28][29][30][31][32][33][34]. Neutrinos produced in the interactions of these CRs would escape freely, while coproduced GeV and TeV gamma rays would be reprocessed to lower energies in the intense radiation field of the accretion disk before escaping. Indeed, the measured intensity of the extragalactic gamma-ray background (EGB) [35,36] favors such scenarios over gamma-ray transparent neutrino sources [37]. Production of the observed astrophysical neutrino flux in the cores of AGN is compatible with the measured EGB level for both pp [29] and pγ [17,25,26] production scenarios.
A contribution of AGN cores to the flux of high-energy neutrinos could also be consistent with the most recent allsky search for the sources of high-energy neutrinos using 10 years of IceCube data [38]. The nearby, Seyfert 2 Galaxy NGC 1068 has been found as the most significant point in the sky at a post-trial significance of 2.9σ. The best-fit neutrino flux measured by IceCube for this source exceeds the observed gamma-ray emission by Fermi-LAT [39], giving support to the hypothesis that the sources of highenergy neutrinos could be opaque to high-energy gamma rays in the GeV-TeV range [40].
In this work, we probe two scenarios of high-energy neutrino production in AGN cores: the first scenario assumes that the neutrino production is proportional to the accretion disk luminosity of the AGN [17,25,26], e.g., due to the acceleration of cosmic rays in the disk's corona [27,30]. The accretion disk emits predominantly in the UV wavelength. However, a strong correlation between UV and x-ray luminosities has been observed in AGN, indicating a connection between the primary radiation from the disk and the soft x-ray emission from the hot electrons in the corona [41,42]. Therefore, the soft x-ray flux measured in two allsky surveys (see below) is used as a proxy for the accretion disk luminosity and the expected neutrino flux in this work.
The second scenario assumes that CR acceleration is suppressed in the AGN with the highest luminosities, favoring instead low-luminosity AGN (LLAGN) for the neutrino production [28,40]. LLAGN are commonly attributed with radiative inefficient acceleration flows (RIAFs) which are formed instead of an optically thick Shakura-Sunyaev accretion disk [43], when the mass accretion rate into the super-massive black hole is relatively small, less than ≈1% of the Eddington accretion rate. For such low accretion rates, the nuclei in the accretion flow do not thermalize and therefore provide a natural seed population for particle acceleration by, e.g., stochastic processes or magnetic reconnection. The breakdown of this seed population for higher accretion rates would favor particle acceleration in such LLAGN over their more powerful siblings, and the neutrino flux would only be proportional to the accretion disk luminosity until the accretion rate limit for RIAFs is reached [40]. In either scenario, pp interactions and pγ interactions can be relevant for generating neutrinos, and gamma rays would not escape because of the large optical depths of the accretion disk for photon pair production, thus both models should be regarded as gamma-ray opaque neutrino source models [40].

II. ACTIVE GALACTIC NUCLEI SELECTION
In order to test the two AGN core models, three samples of AGN have been defined, aiming at a minimal contamination from stellar sources and other galaxy populations.
The soft x-ray emission is an excellent probe of accretion in AGN, being produced by processes related to the accretion disk. However, based on x-ray data alone, it is challenging to eliminate the x-ray binaries and galaxies with x-ray emission associated with star formation rather than an AGN from the sample [44]. For this reason, the soft x-ray data is combined with radio and mid-IR wavelengths for the creation of the AGN samples. On one hand, by crossmatching x-ray catalogs with radio catalogs, it is possible to remove non-AGN sources and to selected only luminous AGN, since luminous radio sources are mostly AGN [45,46]. Radio emission in AGN is produced by synchrotron emission of electrons, which can be due to processes related to the accretion disk and/or large-scale radio jets [47,48]. At the 1.4 GHz, for flux densities of ∼10 mJy and above, AGN are strongly dominating over star-forming galaxies [48]. Moreover, AGN radio emission is unaffected by dust obscuration, and thus relatively unbiased with respect to orientation.
On the other hand, the IR wavelength combined with the x-ray data allows us to select AGN with emission produced in the core, since the dust surrounding the accretion disk reprocesses the emission of the accretion disk into the IR. The IR waveband is relevant for the identification of many AGN in the Universe that have remained hidden from short-wavelength surveys because of reddening and obscuration by dust in and around their nuclei. This emission dominates the AGN spectrum from wavelengths longer than ∼1 μm up to a few tens of microns [47]. It is particularly prominent in the mid-IR (MIR; 3-50 μm) regime, which also contains less stellar contamination. Therefore, IR colors, the ratios of intensities between several mid-IR bands, as collected by WISE [49], can be used to separate AGN from other source types [49,50]. Accordingly, we define two samples of luminous AGN, denoted as radio-selected AGN and infrared-selected AGN. Since they are selected using different criteria (i.e., radio vs IR), they allow us to test the same hypothesis (i.e., neutrino emission from the core of luminous AGN) and to make sure the analysis provides a coherent result once it is corrected for the different sky coverage and x-ray contribution of each of the two samples. IR colors can further be used to distinguish between luminous AGN and LLAGN. We use the ratio between the W1 (∼3.4 μm) and W2 (∼4.6 μm) WISE color bands to define a subset of the infrared-selected AGN that are likely of the low-luminosity type and form our third sample, the LLAGN. Additional details about the selection criteria are presented in Appendix A.
The cross-matching between sources in the different catalogs used is performed using the EXTCAT code [51]. The primary x-rays catalogs used are the ROSAT All-sky Survey (2RXS; [52]), and the second release of the XMM-Newton Slew Survey (XMMSL2; [53]). They have been cross-matched to ALLWISE counterparts in [50] and provide 106,573 (17,665) x-ray sources from the 2RXS (XMMSL2) surveys with ALLWISE IR counterparts [49], covering ∼95% of the extragalactic sky (jbj > 15°). The radio-selected AGN sample was compiled by crossmatching the 2RXS and XMMSL2 sources in this catalog with the NRAO VLA Sky Survey (NVSS; [54]). To avoid biases from the potential neutrino emission of gamma-ray blazars for this analysis, the three obtained AGN samples are further cross-matched with the 3LAC Fermi-LAT catalog [55] to remove all known gamma-ray blazars from the final samples. Finally, all sources below a declination of δ < −5°are discarded, as this part of the sky is not covered by the sample of IceCube events used in this analysis, and IceCube's sensitivity weakens rapidly towards the southern hemisphere.
The radio-selected AGN sample consists of 9749 sources with an estimated contamination from non-AGN sources of only ∼5% and an efficiency of selecting AGN of ∼94% (see Appendix A 1 for more details). It covers ∼55% of the sky. The IR-selected AGN sample is the largest sample in this analysis, and consists of 32249 sources shown in Fig. 1. The contamination from non-AGN sources here is ∼6%, for an efficiency of selecting AGN of ∼89%. The LLAGN sample is a subset of the IR-selected AGN sample. A normalized parameter has been defined based on the IR intensity ratios in the WISE W1 and W2 bands, named Seyfertness, to distinguish Seyfert-type galaxies which are commonly attributed as LLAGN from their more luminous counterparts (see Appendix A 3 for details). Only AGN with a Seyfertness ≥ 0.5 are accepted for the LLAGN sample, resulting in a total number of 15,887 sources for this sample. All three AGN samples are distributed over ∼53% of the sky.
The selection of the sources based on IR color ratios, in particular efficiency, contamination and the Seyfertness parameter, has been cross-validated using the 20% of the sources in the 2RXS catalog that also have counterparts in the VERONCAT [56] catalog, where spectroscopic classifications for each object can be found.
There is, expectedly, significant overlap between the three AGN samples, as shown in Fig. 2. About 82% of the radioselected AGN sources are also found in the IR-selected AGN sample. The LLAGN sample, which itself is a subset of the IR-selected AGN sample (i.e., 100% overlap), has about ∼23% of its sources in common with the radio-selected AGN sample. Table I summarizes the properties of the three AGN samples created for this work, including the cumulative x-ray flux from all sources in the respective sample and the completeness. Completeness is defined here as the ratio between the cumulative x-ray flux included in the sample and the total x-ray flux expected from all AGN in the Universe, estimated using their x-ray luminosity function (luminosity-dependent density evolution model [57][58][59]). The completeness allows an estimation of the contribution from sources not included in the sample, and to extrapolate the analysis results below to the full AGN population. See also Appendix C for details on the calculation of the completeness factors.

III. ANALYSIS
A stacking analysis is performed to search for the cumulative signal from each of the defined AGN samples [60], using a neutrino event sample of about 497,000 upward-going neutrinos, collected in eight years of IceCube operations. Details about this sample are given in [61]. The sample includes only muon-neutrinos to obtain the necessary pointing accuracy and from declinations δ > −5°i n order to reduce the background of atmospheric muons from cosmic-ray air showers. An unbinned maximum likelihood ratio test is performed, to obtain the best fit for n s , the number of signal events, and γ, the index of the energy spectrum of the signal events assuming a single power-law shape. Both a signal and a background probability density function (PDF) enter into the likelihood function [Eqs. (3) and (4) in [61]], and are constructed from Monte Carlo simulations as in [61].
In a stacking analysis, the total signal PDF is given by the weighted sum of the signal PDFs for the individual AGN. They enter into the signal PDF weighted with their expected relative contribution to the neutrino flux [6]. By weighting the sources and normalizing the PDF distribution, we make sure to correlate neutrinos to each AGN with negligible risk of source confusion that could be caused by the almost isotropic distribution of the AGN in the three samples. As described above, the soft x-ray flux reported in the catalogs summarized in Table I is used as a proxy for the accretion disk luminosity and expected neutrino flux.
For each of the three AGN samples, the likelihood function is maximized with respect to n s and γ. The log of the likelihood ratio between the best-fit hypothesis and the null hypothesis (n s ¼ 0) forms the test statistic (TS). The obtained TS value is compared to the TS distribution generated entirely by background-only simulations to derive the p value for the IceCube observations to arise from background fluctuations alone. The p value is defined as the fraction of the generated background-only samples that has an equal or larger TS value than that obtained from the real data. Simulations with different injected signal strengths are then used to calculate the confidence intervals for n s and γ.
In order to minimize statistical bias, the analysis is performed blindly. Therefore, the event selection and the analysis method procedures are fixed before looking at the data. Table II shows the obtained TS and the corresponding p values for the three AGN samples, along with the best-fit parameters obtained from the fit: the number of neutrinosn s and the spectral indexγ. All three tested samples show mild excesses pre-trials relative to background expectations, the largest corresponding to a 1.59 × 10 −3 (2.95σ) p value for the IR-selected AGN sample. The pretrial p values obtained for the radio-selected and LLAGN samples are 0.03 and 0.08, respectively. For each AGN sample, the constraints on the best-fit parameters are shown as profile likelihood contours TABLE I. Properties of the AGN samples created for the analysis. The surveys used for the cross-match to derive each sample, the final number of selected sources, cumulative x-ray flux in the 0.5-2 keV energy range from the selected sources [50] and the completeness (fraction of total x-ray flux from all AGN in the Universe contained in the sample) are listed. To account for testing three alternative hypotheses, the p values have to be corrected by trial factors to obtain the final significance of the observation. Since the three AGN samples are not independent but have significant overlap (see Fig. 2), simulations are necessary to determine these factors instead of a simple analytic correction. In order to correctly take into account the correlation between the three AGN samples, Oð10 4 Þ background-only simulations have been performed. For each simulation, three p values are calculated, one for each AGN sample, and the minimum p value is selected. The distribution of the minimal p values from each background-only simulation is then compared with the pretrial p values measured in the analysis. The resulting post-trial p values are shown in Table II. The best-fit spectrum for the IR-selected AGN sample, for which the background-only hypothesis is rejected at 2.60σ significance after accounting for multiple trials, is shown in Fig. 6 as a quasidiffuse flux (i.e., the cumulative flux from all objects in the sample divided by 2.2π to account for the sky coverage), in units of intensity. Excesses observed in the other two AGN samples are not significant after trials correction. Therefore, upper limits are set for the LLAGN and radio-selected AGN samples, as shown in Fig. 7. The upper limits are calculated for an E −2 energy spectrum between 30 TeV and 10 PeV,  [62,63] are plotted as a differential flux unfolding using 95% C.L. The best-fit 1σ contour is scaled by a correction factor that takes into account the flux from unresolved sources (completeness of the sample). Systematic uncertainties and the error on the completeness factor are not included. The models from [26] (dashed, gray line) and [25] (dotted, gray line) are overlaid for comparison. which was found to be the energy range to which this analysis was sensitive. For each upper limit and best-fit flux we determine the valid energy range ([30 TeV, 10 PeV]) according to the procedure in [6] (see Appendix B for details). This energy range specifies where IceCube has exclusion power for a particular model. However, given the different spectral index and that the energy range has been calculated based on Monte Carlo simulations, it does not coincide with that of the diffuse flux limits, since in real data no neutrinos with energy greater than few PeV are detected.

IV. RESULTS
The best-fit flux and the upper limits shown in Figs. 6 and 7 are scaled from the results for the individual samples using the completeness of the catalogs to display the respective values for the total population of AGN. This allows a comparison to the observed astrophysical diffuse neutrino flux measured by IceCube in [62,63]. If the result obtained for the IR-selected AGN sample is interpreted as an upper-limit, a neutrino flux up to 1.97 × 10 −18 GeV −1 cm −2 s −1 sr −1 is expected at 100 TeV from the cores of the entire population of luminous AGN, assuming an E −2 energy spectrum (see Sec. V for a comparison to the diffuse neutrino flux observed by IceCube).
The figures showing best-fit flux and upper limits do not include systematic uncertainties. Following [61], the total systematic uncertainty on the neutrino flux normalization is estimated to be 11%. The major contribution to the systematic error comes from uncertainties on the optical efficiency of the IceCube sensors and on the optical properties of the Antarctic ice.

V. SUMMARY AND DISCUSSION
We have presented an analysis probing the origin of astrophysical neutrinos by searching for a correlation between the cores of AGN and eight years of IceCube neutrino data. Two complementary models for neutrino production have been tested in this paper: one that favors neutrinos to be produced in the geometrically thin, optically thick accretion disks of luminous AGN, and one that predicts the bulk of the neutrino emission from the RIAF of LLAGN. In total, three AGN samples, each one consisting of Oð10 4 Þ sources, have been compiled using radio and IR survey data to identify AGN, and distinguish low-luminosity from high-luminosity objects.
The soft x-ray flux obtained from the 2RXS and XMMSL2 catalogs is used as a proxy for the accretion disk luminosity and expected neutrino emission. Each one of the (statistically not independent) AGN samples shows a positive correlation to the neutrino data, however for the LLAGN it is weak and compatible with no correlation within 1 standard deviation. The IR-selected AGN sample shows the strongest indication for a correlation, with a significance corresponding to 2.60 standard deviations after accounting for trial factors from studying more than one sample. The best-fit spectrum of the correlated events, assuming a power-law shape, has a spectral index close to 2 for all studied samples, as expected for particle acceleration scenarios in cosmic environments, and much harder than the background of atmospheric neutrinos. The fact that the three best-fit spectral indexes are so close to each other is not unexpected, since the samples are highly correlated (see Fig. 2). However, this spectral index is significantly harder than the index seen from IceCube diffuse flux measurements [62,63]. This implies that the IceCube diffuse flux might arise from multiple populations of sources with different spectra and that the AGN cores would be responsible for the majority of the emission at the highest energies (> 1 PeV). In this scenario, the other populations contributing to the diffuse flux would have softer spectra [37,[63][64][65][66].
Within the framework of the tested model, i.e., a linear proportionality between accretion disk luminosity (estimated from soft x-rays) and the neutrino flux, the total contribution of AGN to the astrophysical neutrino flux can be extrapolated using x-ray luminosity functions to estimate the contribution of sources not selected in the source samples. The contribution of the IR-selected AGN themselves to the diffuse flux at 100 TeV measured by IceCube [63] amounts to 10 þ5 −4 %. The associated population's total contribution can be 27%-100% at 1σ confidence interval after completeness correction, assuming soft x-ray and neutrino luminosities are correlated. The error on this fraction also includes the error on the completeness, which has been combined with the flux error by a bootstrapping method. This is consistent with a predominant origin of neutrinos at this energy from the cores of AGN, while potentially accommodating subdominant contributions from blazar jets [4] and potentially tidal disruption events [67]. It is also consistent with the contribution extrapolated from the best fit to the radio-selected AGN sample, which tests the same hypothesis, albeit for this sample the correlation is statistically less significant.
The results reported here are the first constraint on the AGN contribution to the IceCube diffuse flux. If confirmed by future analyses, our findings might point to cosmic rays accelerated in the AGN core regions being responsible for the bulk of the astrophysical neutrino flux observed by IceCube above 100 TeV. The significance of this finding is 2.60σ post-trial.
Assuming a genuine signal, the results can be interpreted as the sources of astrophysical neutrinos being dominated by AGN and that they are opaque to GeV-TeV gamma rays. Several qualitative aspects of the present results, the hard spectrum of the neutrinos attributed to the signal, the consistency between AGN samples and with the total neutrino flux observed, create an intriguing scenario. Therefore, it is essential to follow-up on these first indications in the future, using significantly larger IceCube datasets, an improved IceCube angular resolution [68], the next-generation IceCube-Gen2 detector [69,70], and, potentially, deeper soft x-ray surveys, such as the one currently conducted by eROSITA [71].
The radio-selected AGN sample is obtained by correlating the NVSS with the 2RXS and the XMMSL2 catalogs, keeping all sources whose radio and x-ray positions differ by less than 60 arcsec. The search radius has been chosen based on the x-ray source positional errors, which are energy dependent and on average ∼20 arcsec (∼10 arcsec) for the 2RXS (XMMSL2) sources [50]. About 6.12% (5.93%) of the 2RXS (XMMSL2) objects have more than one radio counterpart, and for this case the closer NVSS source to the x-ray counterpart is chosen.
The x-ray sources with a radio counterpart are then crossmatched with the 3LAC catalog [55] in order to remove blazars, using the 95% source position error of the gammaray sources as search radius. The largest group of objects in the x-ray/NVSS-3LAC correlation are the BL Lacs, followed by FSRQs and blazars of uncertain type. Only a few percent of the sources are nonblazar AGN. Before combining the NVSS/2RXS and NVSS/XMMSL2 samples, we remove duplicated x-ray sources (i.e., XMMSL2 sources already included in the 2RXS catalog) and convert the x-ray fluxes to the common 0.5-2 keV energy range. The final radio-selected AGN sample contains 13,927 sources. A further geometrical cut is applied in the end to select only sources in the northern hemisphere (declination > −5 deg), resulting in the final 9749 sources (see Fig. 8). The contamination from non-AGN sources and efficiency of the selection have also been calculated by cross-matching the NVSS catalog with two randomized x-ray catalogs reproducing the density distribution of the 2RXS and XMMSL2 catalogs. It results in a contamination of 5% and an efficiency of ∼94%.

IR-selected AGN sample
In this paragraph, we describe in detail the cuts applied on the original 2RXS and XMMSL2 catalogs from [50] to obtain the IR-selected AGN sample. They are based on the following IR-color magnitudes: W1 is defined in the 3.4 μm band, W2 at 4.6 μm, and W3 and W4 at 12 and 22 μm, respectively. The first cut we apply is based on the x-ray/MIR relation suggested by [50]. It is seen that the x-ray sources that verify the empirical relation are mainly AGN, which can thus be separated from galaxies and stars over six orders of magnitude. It is also seen that most of the sources below the relation are also stars based on their ALLWISE colors. We confirm the validity of this cut by overlapping the position of the AGN counterparts from the VERONCAT catalog [56] to the 2RXS sources in the same plane: most of these AGN verify Eq. (A1), suggesting that most of the sources above this value are indeed AGN. After applying the [W1] cut of Eq. (A1), we can use the ALLWISE [W1-W2] and [W2-W3] of the sources for their qualitative characterization as in [49] and to isolate the position of the sources we are interested in, removing blazars, starburst, and normal galaxies. Furthermore, we can look at how the VERONCAT counterparts of the x-ray sources are distributed in the same color-color diagram to validate these cuts. The VERONCAT counterparts of the 2RXS sources lay in the area in between the black lines in Fig. 9, that represent the cuts used to isolate the AGN in the [W1-W2] and [W2-W3] plane. This validates thus the 3 color-color cuts and we can therefore apply them on the 2RXS and XMMSL2 samples (see Fig. 9, left panel).
The final cut is applied by looking at the [W1-W2] vs [W3-W4] plot in the right panel of Fig. 9. Also in this case, by overlapping the position of the VERONCAT AGN we discover that the sources above the black solid line are AGN, while below that curve the sources are mostly normal galaxies. Therefore, we apply by-eye one more cut to [W1-W2] as a function of [W3-W4].
After applying all the color cuts, a cross-match with the 3LAC catalog is performed to remove the blazars from the final sample.

LLAGN sample
The LLAGN sample is a subsample of the IR-selected AGN one (see Fig. 2). Only the IR-selected sources with a 2RXS counterpart are selected, since most of them have a counterpart in the VERONCAT catalog. Using the VERONCAT classification, it is possible to separate the bright AGN from the Seyfert galaxies in the [W1-W2] space. We can use these two distributions (after normalizing them) to define a "Seyfertness" PDF as where PðSÞ and PðBÞ are the probability of being a Seyfert or a bright galaxy, respectively, approximated from the W1 − W2 histograms. Since LLAGN are mostly Seyfert galaxies [72], we can use this PDF to give a weight to our sources based on how likely they are to be LLAGN. We assign to each source a Seyfertness PDF between 0 and 1.
In the final sample we only include sources with Seyfertness ≥ 0.5, since at this value we have the best trade-off between efficiency (77%) and contamination (21%) of the selection. The final LLAGN sample contains 25,648 sources and 15,887 sources are at a declination larger than −5 deg (see Fig. 10).

APPENDIX B: ENERGY RANGE CALCULATION
The sensitivity of IceCube depends on the neutrino energy. For this reason, it is important to determine, for a given spectral index, the relevant energy range ½E min ; E max where this analysis is mainly sensitive. Outside this energy range there are few events that contribute little to the TS. In order to estimate it, we progressively change the energy range of the injected neutrino signal. The point at which the sensitivity decreases by 5% indicates the limit for the analysis. The lower energy bound E min is determined by varying E min and only injecting pseudosignal events with an energy E > E min , while E max ¼ 10 PeV remains unchanged. Analogously, we determine the upper energy bound by varying E max and only injecting pseudosignal events with an energy E < E max , while E min ¼ 100 GeV remains unchanged. Figure 11 shows the results for γ ¼ 2, using 100 sources of the radio-selected AGN sample. The energy range where we are more sensitive is [30 TeV,10 PeV]. We expect the same results to be valid also for the other two AGN samples and when stacking more sources.

APPENDIX C: CORRECTION FACTORS FOR UNRESOLVED SOURCES
In order to extend the results of the analysis of each sample to the total radio galaxy and LLAGN populations, we need to take into account also those sources not included in our samples. In fact, some sources can be cut out because of the selections applied for the creation of the samples, either because they are too faint to be detected, or they fall outside the field of view of the x-ray telescopes, or the wavelength of their radiation falls outside the windows investigated here, due to redshift effects.
The correction factor is called "completeness" and is defined as the fraction of the total flux from all sources in the observable Universe that is resolved into individual point sources in the catalogs used for the analysis. The total x-ray flux expected from all AGN has been estimated using the soft x-ray luminosity function in the energy range 0.5-2 keV. Recent x-ray surveys have found that the soft xray luminosity function of AGN is best described by a luminosity-dependent density evolution model, rather than the classical pure luminosity evolution or pure density evolution models, which tend to overestimate the cosmic x-ray background [42]. According to the luminositydependent density evolution model, the luminosity and the number density distributions change simultaneously as a function of redshift.
By integrating the luminosity functions from [57][58][59] over the luminosity and the comoving volume, we derive the source count distribution Nð> SÞ, that is the number of sources detectable at Earth above a minimum photon flux. Figure 12 shows the cumulative number of AGN Nð> SÞ as a function of the x-ray flux S in the 0.5-2 keV energy range. The expected number of sources derived from the luminosity functions (solid lines) is compared to the number of AGN in the three samples used in this analysis. All samples are scaled by their sky coverage, and the LLAGN sample is also scaled by the efficiency of the selection. The completeness factor is obtained by taking the ratio between the area below the three theoretical curves and the area below the distribution of each AGN sample. For the radio-selected and IR-selected AGN samples, the model from [57] is used (left panel of Fig. 12), while the LLAGN sample is better described by the model of [59] (right panel of Fig. 12). The final values for the completeness are shown in Table I.