Multiwavelength search for the origin of IceCube's neutrinos

The origin of astrophysical high-energy neutrinos detected by the IceCube Neutrino Observatory remains a mystery to be solved. In this paper we search for neutrino source candidates within the $90$% containment area of $70$ track-type neutrino events recorded by the IceCube Neutrino Observatory. By employing the Fermi-LAT 4FGL-DR2, the Swift-XRT 2SXPS and the CRATES catalogs, we identify possible gamma, X-ray and flat-spectrum radio candidate sources of track-type neutrinos. We find that based on the brightness of sources and their spatial correlation with the track-type IceCube neutrinos, the constructed neutrino samples represent special populations of sources taken from the full Fermi-LAT 4FGL-DR2/Swift-XRT 2SXPS/CRATES catalogs with similar significance ($2.1\sigma$, $1.2\sigma$, $2\sigma$ at $4.8~\mathrm{GHz}$, $2.1\sigma$ at $8.4~\mathrm{GHz}$, respectively, assuming 50% astrophysical signalness). After collecting redshifts and deriving sub-samples of the CRATES catalog complete in the redshift--luminosity plane, we find that the 4.8 GHz ($8.4$~GHz) sub-sample can explain between 4% and 53% ($3$% and $42$%) of the neutrinos (90% C.L.), when the probability to detect a neutrino is proportional to the ($k$-corrected) radio flux. The overfluctuations indicate that a part of the sample is likely to contribute and that more sophisticated schemes in the source catalog selection are necessary to identify the neutrino sources at the $5\sigma$ level. Our selection serves as a starting point to further select the correct sources.


INTRODUCTION
High-energy cosmic rays (HECRs) are energetic subatomic particles with a non-thermal energy spectrum that ranges from GeV up to at most ZeV energies (e.g., Abraham et al. 2010;Abu-Zayyad et al. alternative is the detection of high-energy gamma-rays arising from hadronic interactions. These are, however, absorbed in interactions with the CMB photons, so the high-energy gamma-ray horizon at TeV energies is constrained to the universe at z < 1, and at PeV energies it is constrained to our Galaxy. In addition, the leptonic production of high-energy gamma-rays makes an interpretation of the photon spectrum alone much more difficult.
Neutrinos are produced in high-energy (photo-)hadronic interaction processes via the production of pions and kaons, which subsequently decay to produce high-energy neutrinos and muons, where the latter further contribute to the neutrino flux (Gaisser et al. 1995;Becker 2008). High-energy gamma-rays are always co-produced via the decay of neutral pions. The IceCube Neutrino Observatory detects cosmic high-energy neutrinos with energies between ∼ 100 GeV and ∼ 10 PeV. Among millions of atmospheric neutrinos and hundreds of billions of cosmic-ray muons, about 200/year cosmic neutrinos for an E −2.5 spectrum (Aartsen et al. 2015) are detected via data reduction by introducing intelligent physicsbased cuts to remove a large portion of the atmospheric particles (IceCube Collaboration 2013a,b, 2014Icecube Collaboration 2016;IceCube Collaboration 2017, 2020. After their first detection by the IceCube Neutrino Detector (IceCube Collaboration 2013a), high-energy astrophysical neutrinos quickly became powerful probes of astrophysics (e.g. Ahlers et al. 2018). In 2018, multimessenger observations of the blazar TXS 0506+056 were reported (IceCube Collaboration et al. 2018). The high-energy neutrino event IceCube-170922A detected at 22/9/2017 was in temporal and positional coincidence with a strong gamma-ray flare shown by TXS 0506+056 detected by Fermi-LAT at GeV energies. Knowing where to look, the IceCube Collaboration analyzed archival track-type data and found 3.5σ evidence for a neutrino flare in the direction of TXS 0506+056 (IceCube Collaboration 2018) during a low gamma state of TXS 0506+06 (Garrappa 2019), that lasted about half a year and was of significantly lower energy as compared to the IceCube event IceCube-170922A. Since 2019 June, neutrino events with high chance to be astrophysical and detected in the framework of the IceCube Realtime Alert System (Aartsen et al. 2017b) are distributed to the Gold (signalness is above 50%) or Bronze channel (signalness is below 50% but above 30%). The signalness is the probability that an observed neutrino is of astrophysical (s = 1) or atmospheric (s = 0) origin (see its definition in Aartsen et al. 2017b).
In general, radio-loud active galactic nuclei (AGN) are among the primary candidate sources of IceCube neutrinos. A prime candidate is their blazar sub-class. These sources usually show low-energy and high-energy bumps coexisting in their spectral energy distribution (SED). Throughout the paper we refer to blazars as the subclass of AGN that point their jet close to our line-ofsight, according to the unification scheme of radio-loud AGN (Urry & Padovani 1995). The flat spectrum at the 5 GHz range of radio frequencies is indicative of these sources pointing to us and boosting the emission (Biermann et al. 1981;Gregorini et al. 1984).
High-energy neutrinos can be produced in the interaction of cosmic rays with either a photon or a gas target. Here, threshold of the kinematic energy of a cosmic-ray proton for interactions with the gas is at E p 280 MeV. Interactions with a photon target, on the other hand, is typically constrained to higher energies, as its threshold is dominated by the production of the delta-resonance at E 1.3 GeV. As photon targets are typically most dense for thermal emission spectra in the infrared to optical energy range, these interactions might be very effective in certain regions like close to the accretion disk of the supermassive black hole, but the energy spectrum will be limited to energies above ∼ 10 18 eV for the Bethe-Heitler pair production and above ∼ 5 × 10 19 eV for the photopion production.
Two photohadronic processes are of astrophysical interest as mentioned, the Bethe-Heitler pair production and the photopion production, such that neutrinos are produced only by the latter process. Enhanced gamma-ray emission from the decay of neutral pions is bound to accompany the emission of neutrinos born in charged current pion production. However, no statistically significant association of neutrino events with an ensemble gamma-ray loud blazars has been found yet (e.g. Aartsen et al. 2017a;Neronov et al. 2017;O'Sullivan & Finley 2019;IceCube Collaboration 2020). Moreover, at lower energies there is more flux in the IceCube diffuse neutrino component than that could be expected from the Fermi diffuse gamma-ray sky if they have the same sources (IceCube Collaboration 2021).
Recently Bartos et al. (2021) found blazars (mainly FS-RQs, based on their assumption of the blazar density) cannot contribute more than 11% to the total diffuse neutrino flux (90% confidence level, C.L.). We note these results are model dependent as regards the diffuse backgrounds.
Hadronic gamma-rays produced in photonrich environments, such as AGN, may lose energy through electromagnetic cascading, which results in a high flux emitted in the hard Xray to MeV range (e.g. Murase et al. 2018;Rodrigues et al. 2019;Halzen & Kheirandish 2020).  had shown that if relativistic protons interact with the synchrotron blazar photons producing gamma-rays through photopion processes, despite being a subdominant proton cooling process (for protons satisfying the threshold condition for photopion production), the typical blazar SED should have an emission feature due to the synchrotron emission of Bethe-Heitler pairs in the energy range 40 keV-40 MeV. Therefore, typical neutrino emitters would show this third, although less energetic bump between the typical low-energy and high-energy bumps of blazar SEDs (e.g. Abdo et al. 2010).
It is a nontrivial task to find a well-established connection between the radio flux and neutrino flux in AGN, mainly because neutrino production zones may be realized in different parts of an AGN (central or extended parts of the galaxy), and the resolution of the instruments, especially at high energies, is somewhat poor. While a significant connection is not widely concluded yet, catalog searches and studies of particular individual radio-loud sources associated with neutrinos (e.g. Kadler 2016; Kun et al. 2019;Britzen et al. 2019;Ros et al. 2020;Britzen et al. 2021;Kun et al. 2021;Rodrigues et al. 2021) deliver more and more results that suggest radio observations play an important role in the understanding of the sources of astrophysical neutrinos.
One way to study the connection between radio observations of AGN and neutrino data observed with Ice-Cube is the search for correlations between neutrinos that have a high probability of being astrophysical and individual sources or source catalogs. By analyzing public muon-track data from the 7 years of IceCube observations against the 8 GHz VLBI data of 3, 411 radio-loud AGN in the Radio Fundamental Catalog, Plavin et al. (2021) found a connection between the neutrinos and the radio sources at a 4.1σ significance level, combining the post-trial significance of their study and the results in Plavin et al. (2020). Hovatta et al. (2021) studied association of IceCube neutrinos with radio sources observed at Owens Valley (at 15 GHz) and Metsähovi Radio Observatories (at 36.8 GHz) by identifying sources in their radio monitoring sample that are positionally consistent with IceCube high-energy neutrino events. According to their results, observations of large-amplitude radio flares in a blazar at the same time as a neutrino event are unlikely to be a random coincidence at 2σ level.
Another way to study the connection between radio observations of AGN and neutrino data taken with Ice-Cube is the concept of stacking, see Achterberg et al. (2006) for the introduction of the method of source stacking, for the first time in neutrino astronomy. Zhou et al. (2021) studied the connection between 10 year of IceCube muon-track data (with 1, 134, 450 tracks, taken between April 2008 and July 2018), and 3, 388 radio-bright AGN at 8 GHz selected from the Radio Fundamental Catalog. Using the unbinned maximum likelihood method they analysed the neutrino data in the position of radio-bright AGN. Their results imply the 3, 388 radio-bright AGN can account for at most 30% (95% C.L.) of the flux of neutrino tracks in the 10 years catalog of IceCube. Desai et al. (2021) presented a correlation analysis between 15 GHz radio observations of AGN reported in the MOJAVE XV. Catalog and 10 years of IceCube detector data. They did not find a significant correlation between them.
In this paper we carry out a multiwavelength search of sources within the 90% containment area of 70 IceCube track-type neutrino events. After the Introduction in Section 1., in Section 2. we design an analysis to match the 70 track-type neutrino events published by IceCube and the Fermi-LAT 4FGL-DR2 catalog, the Swift -XRT 2SXPS catalog and the CRATES Catalog. We perform a randomness test of the average brightness of the candidate samples in Section 3. We investigate neutrino-radio statistical correlations in Section 4. by deriving number counts of the CRATES catalog after collecting redshifts of more then 3, 000 radio sources, and by defining a 4.8 GHz and a 8.4 GHz sub-sample that are complete in radio luminosities. Then we estimate how many neutrinos can be explained by the above sub-samples, assuming the probability to detect a neutrino is assumed to be proportional to the radio flux. We summarize our results in Section 5. and discuss them in Section 6. We give our conclusions in Section 7.

POINT SOURCES WITHIN 90% CONTAINMENT OF ICECUBE NEUTRINOS
The neutrino list compiled by Giommi et al. (2020) summarizes 70 public IceCube high-energy cosmic neutrino tracks (recorded between 2009-08-13 and 2019-07-30) that are well reconstructed and off the Galactic plane. We use this sample to find neutrino candidate

TXS 0506+056
Fermi/LAT 4FGL-DR2 point sources and track-type IceCube neutrinos  Giommi et al. (2020). The map is shown in Mollweide projection. The equator is plotted by a purple continuous line. TXS 0506+056 and PKS 1502+106 are marked on the map. The 7 brightest sources (with F0.1−100 GeV > 10 −7 ph cm −2 s −1 ) are not shown to enhance the visibility of the flux-diversity of the sources. Table 1. Excerpt of the neutrino source candidates from 10-year Fermi -LAT Point Source Catalogue (4FGL-DR2) (n = 29, out of 5, 787). Neutrino ID (1), detection time (2), arrival direction of neutrino (90% containment, 3-4), 4FGL-DR2 source ID (5), equatorial coordinates (J2000) (6-7), predicted number of photons in the 50 MeV-1 TeV energy range (8), flux between 0.1-100 GeV (9). The full table is available in electronic format here. sources. Given that the errors on the neutrino sky positions are asymmetric, the center of their error-ellipses were shifted in the direction of the most significant errors by an amount equal to half the difference between the larger and smaller errors, and by setting the major and minor axes equal to the sum of the two asymmetrical errors.
2.1. Gamma-ray neutrino source candidates in the Fermi-LAT 4FGL-DR2 catalog The Fermi Gamma-ray Space Telescope is one of the most successful telescopes in the gamma range. Its main instrument, the Large Area Telescope (LAT) is sensitive between 20 MeV and 300 GeV energies, though measurements are made up to 1 TeV. The 4FGL-DR2 catalog includes 5, 787 gamma sources based on 10 years of survey data in the 50 MeV-1 TeV energy range (Ballet et al. 2020), being an incremental version to the 8 years 4FGL catalog (Abdollahi & et al. 2020). More than 3,413 of the identified or associated sources are blazars in the Fermi-LAT 4FGL-DR2 catalog (Abdollahi & et al. 2020;Ballet et al. 2020). We plot this catalog in Fig.  1 overlaid on the 90% containment area of the Ice-Cube track-type neutrino events (Giommi et al. 2020). Table 2. The number of sources (N), the average (¯), standard deviation (σ), median (med), minimum (min) and maximum (max) values of the flux the sources in the total (T) Fermi -LAT 4FGL-DR2 (F0.1−100 GeV ), of the mean rate of sources in the full Swift-XRT 2SXPS catalog (Rm), of the flux density at 4.8 GHz (S4.8) and 8.4 GHz (S8.4) of sources in the full CRATES catalog, as well as the same properties of the sub-samples of neutrino candidate sources (S) derived from the above catalogs by the list of IceCube track-type neutrino events compiled by Giommi et al. (2020). We found 29 gamma-ray sources from the Fermi-LAT 4FGL-DR2 catalog within the 90% containment area of track-type neutrinos, as summarized in Table 1, and from now on we refer to them as the gamma "ν-sample". In Table 2 we summarize the average, standard deviation, median, minimum and maximum values of the gamma flux in the 0.1-100 GeV energy range of the sources in the full Fermi-LAT 4FGL-DR2 catalog, and separately such values for the gamma ν-sample. The average flux of the Fermi-LAT ν-sample emerged as 2.4 × 10 −9 ph cm −2 s −1 .
Similar studies were already conducted in the gamma regime, but either with different neutrino lists or with different source catalogs. Franckowiak et al. (2020) studied the IceCube real-time alerts, that selects highenergy ( 100 TeV) starting and through-going muon track events (Aartsen et al. 2017b), and archival neutrino events that would have passed the same selection criteria. We used the neutrino list by Giommi et al. (2020) that is partially different. When searching for the neutrino source-candidates, Giommi et al. (2020) used the Fermi-LAT 3FHL (that reports hard Fermi-LAT sources significantly detected in the 10 GeV-2 TeV energy range during the first 7 years of the Fermi mission using the Pass 8 event-level analysis, Ajello et al. 2017), 4LAC (the Fourth Catalog of Active Galactic Nuclei detected by the LAT based on 8 years of data, Ajello et al. 2020) and 3HSP catalogs (a sample of extreme and high-synchrotron peaked blazars and blazar candidates, Chang et al. 2019), while we employed the whole 10 years catalog of Fermi-LAT data.

X-ray neutrino source candidates in the Swift-XRT 2SXPS catalog
The Swift -XRT Point Source (2SXPS) catalog covers 3, 790 square degrees on the sky and contains position, fluxes, spectral details and variability information for 206, 335 X-ray point sources. Observations were made with the X-ray telescope on-board the Neil Gehrels Swift Observatory between 2005 Jan 01 and 2018 August 01 (Evans et al. 2020). Each source has a detection flag which indicates how likely it is to be a real astrophysical object. We applied several selection criteria: DetFlag=0 (good), SNR(=Count rate/BG rate)≧ 10 (in the total 0.3-10 keV energy band). We plot the catalog in the upper panel of Fig. 2 overlaid on the 90% containment area of the IceCube track-type neutrino events. Due to the usually wide diversity of sources with different nature in an X-ray catalog such like the 2SXPS (with typically AGN and stars as the most dominant populations), we applied an additional selection criterion that allows only AGN being in our sample. Tranin et al. (2022) applied a robust probabilistic classification of X-ray sources in the Swift -XRT and XMM-Newton catalogs. We selected such X-ray sources, for which the posterior probability of being an AGN is > 0.99 according to their probabilistic classification. The resulting sub-sample contains 9, 079 sources having high chance of being an AGN, from which we found 61 sources to be located within the 90% containment area of track-type neutrinos, as summarized in Table 3.
In Table 2 we summarize the average, standard deviation, median, minimum and maximum values of the X-ray full band rate in the 0.3 keV-10 keV energy range of the sources in the full 2SXPS catalog, and separately such values for the X-ray ν-sample. The the average flux of the Swift -XRT ν-sample emerged as 0.049ctss −1 .

Radio neutrino source candidates in the CRATES catalog
We searched radio sources within the 90% containment area of neutrinos in the Combined Radio All-Sky Targeted Eight GHz Survey (CRATES, Healey et al. 2007). The authors have assembled an 8.4 GHz survey of bright, flat-spectrum (α > −0.5, with the convention S ν ∼ ν α , Gregorini et al. (1984)) radio sources with nearly uniform extragalactic (|b| > 10 degrees) coverage for sources brighter than a 4.8 GHz flux density of S 4.8 GHz = 65 mJy. The catalog is assembled from existing observations, especially the Green Bank 6cm Radio Source Catalog (GB6, 4.8 GHz Gregory et al. 1996), the Parkes-MIT-NRAO Survey (PMN, 4.8 GHz Wright et al. 1994;Griffith et al. 1994), and the Cosmic Lens All-Sky Survey (CLASS, 8.4 GHz Myers et al. 2003). The existing observations were augmented by reprocessing of archival VLA and ATCA data and by new observations to fill in coverage gaps. The resulting catalog provides precise positions, sub-arcsecond structures, and spectral indices for 14, 467 sources at 8.4 GHz. The number of sources exceeds the number of sources at 4.8 GHz since in many cases there are multiple 8.4 GHz counterparts to a single 4.8-GHz source due to the higher frequency and therefore better resolution. There are also 762 sources for which no 8.4 GHz counterpart was detected. We plot the catalog in Fig. 3 overlaid on the 90% containment area of the IceCube track-type neu- CRATES flat spectrum radio sources and track-type IceCube neutrinos  Table 4. Excerpt of the neutrino source candidates from CRATES catalog of flat spectrum radio sources (n = 100, out of 14,467 at 8.4 GHz). Neutrino ID (1), detection time (2), arrival direction of neutrino (90% containment, 3-4), ID of the CRATES radio source (5), its position (6-7), its detected flux at 4.8 GHz (8) and 8.4 GHz (9) observing frequencies, morphology type (10) as N = No detection at 8.4 GHz, P = point source, S = short jet (<= 1" for VLA maps), L = long jet (> 1"). The full table is available in electronic format here. GHz, and separately such values for the radio ν-samples.
The the average flux of the CRATES ν-samples emerged as 544 mJy at 4.8 GHz, and 423 mJy at 8.4 GHz.

RANDOMNESS TEST OF THE ν-SAMPLES
In this section we investigate whether the selection method based on track-type IceCube neutrinos chooses random ν-samples, or ν-samples with flux properties significantly different from the full catalogs.
The null-hypothesis states a ν-sample represents a sample randomly chosen from the full catalog, while the alternate hypothesis states the average observed flux of a ν-sample is rather special value characteristic to the selection based on the neutrinos.
First we took the angular error of the 70 neutrinos, randomly generated their RA coordinates (with 0.1 decimal precision), repeated the source search described in Section 2.1, calculated the average flux contained by the gamma sources, and repeated this process 10, 000 times. We did not scramble the declination of the neutrinos because of the strong variation of the sensitivity of IceCube with declination. Integrating the probability density functions (PDFs) from the observed averages to infinity we get the probability of the corresponding ν-sample to be random under the null hypothesis. We repeated the same process for the Swift -XRT 2SXPS catalog, and for the CRATES catalog at 4.8 GHz and 8.4 GHz.
We repeated the tests for different values of the signalness. For example, when we set the signalness to be 0.6, each of the neutrinos has 60% probability to be of astrophysical origin relative to the total background rate. Then the chance of a neutrino to enter a cross-match equals the actual signalness of the test. We plot the resulting p-values as function of the signalness in Fig. 4. It seems that increasing the signalness of the neutrinos, the p-value of the Fermi-LAT 4FGL-DR2 and CRATES catalogs first decrease, then the p-values fluctuate about constant values. Increasing the signalness, however, the p-value of the Swift-XRT 2SXPS continuously decrease without reaching a constant significance level.
We found that there is a reasonable, though not conclusive connection between the neutrinos and the Fermi-LAT 4FGL-DR2 gamma sources (p ≈ 0.016 for 50% signalness, with corresponding significance ∼ 2.1σ), as well as between the neutrinos and the CRATES radio sources at 4.8 GHz and at 8.4 GHz observational frequencies (p ≈ 0.025, p ≈ 0.019, respectively, corresponding significances ∼ 2.0σ, ∼ 2.1σ both for 50% signalness). We note, the 50% signalness (s=0.5) represents the cut for which neutrino events detected before 2019 June (since when events with high signalness are unified to Gold and Bronze channels) probably would end up in the Gold channel of IceCube's Real Time Alert System, which is why we chose this value in accordance to previous work (Giommi et al. 2020). A less significant connection emerged between the track-type neutrinos and the Swift -XRT 2SXPS catalog (p ≈ 0.11 for 50% signalness, with corresponding significance ∼ 1.2σ). Increasing the signalness of the neutrinos to 100%, however, the significance of the subsample derived from the Swift-XRT 2SXPS catalog to be the origin of the neutrinos reaches of a p-value of 0.04 (∼ 1.8σ).
In the rest of the paper we focus on the radio-neutrino argument. We will present our study focusing on the Xray and gamma regimes in a subsequent paper.

INVESTIGATION OF NEUTRINO-RADIO CORRELATIONS BASED ON THE CRATES CATALOG
In this section we estimate how many neutrinos can be explained with the CRATES catalog, assuming the neutrino flux is proportional to the radio flux.

Error estimates of the flux density and position of CRATES sources
We needed the flux density uncertainties; however, this information is not directly accessible through the CRATES catalog 1 . The 4.8 GHz part of CRATES is composed mainly from the PMN and GB6 surveys. The single-dish observations of the southern sky at 4.85 GHz are available from the PMN catalog, which covers the region −87 • < δ < +10 • . Observations were made with the Parkes 64m Radio Telescope. The standard error (dS) of the flux densities (S) in the PMN catalog can be calculated as : where δ is the declination. We used this equation to estimate the flux density uncertainties of the radio sources in CRATES south of δ < 0 • , at 4.8 GHz.
The GB6 catalog of sources covers the declination range 0 • < δ < +75 • , and the observations were made with the 91 m telescope at Green Bank. Since the noises and calibration errors of the map are unknown (to calculate the rms uncertainty in the corrected peak flux density and the source integrated flux density), we fitted a linear function to estimate the error in the flux density of sources in the region 0 • < δ < +75 • at 4.8 GHz, where CRATES is based on mostly GB6 survey. The maximum declination of neutrinos ∼ 67.40 • (IC-140216A) is well south of the maximal declination covered by the GB6 catalog, we do not need to know flux density in the far north region δ > 75 • . The CLASS survey contains observations made with the VLA at A configuration at 8.4 GHz. Augusto et al. (1998) mentioned the typical error on the CLASS flux densities as 5%. Since CLASS is a parent catalog of CRATES, and most of the sources in CRATES were observed by the VLA (11, 333 out of 14, 467), we estimated the errors as 5% of the flux densities.
In the north region (0 • < δ < +75 • ) the pointing errors of GB6 source positions are estimated as (Gregory et al. 1996): At 8.4 GHz, CRATES is largely based on two surveys, the CLASS survey in north, and the (unpublished) PMN-CA survey in the south. The equatorial south was surveyed for gravitational lenses, and X-band data are available in this region. In the CLASS region 0 < δ < +75 • , |b| < 10 • the typical angular resolution error is about 0.2" (Myers et al. 2003). In the south the PMN-CA survey quotes a typical position error of 0.6" in each coordinate for the PMN-CA data.

Number counts of CRATES sources in the z-samples
Our first goal was to find the luminosity down to which our sample is complete at a given redshift, which we call "z-sample". To construct a complete z-sample, we defined the parameter C(L, z) = N o (L, z)/N p L, z), where N o (L, z) (N p (L, z)) is the observed (predicted) number of sources in a given luminosity and redshift bin.
We found the redshift of 3, 467 CRATES sources in the Simbad database using the positional uncertainties summarized in Section 4.1, and by running Astropy queries (Astropy Collaboration et al. 2013. When more objects were found within the errorbars, the object closest to the CRATES position was assumed (less then 0.5% of the cases). The majority of these sources are quasars (2, 262 sources), BL Lacs (486) and blazars (97). The maximum of their observed number counts happens between z ≈ 0.75 and z ≈ 2 (see Fig. 5). It is important to note that nearby blazars like the Markarian sources, galaxies with extraordinary UV brightness compared to normal galaxies, are much closer to us: their redshift ranges from 0.0009 to 1.912 with median value 0.0238, situating them as residents of typically the Local Universe. The First Byurakan Survey (Markarian) Catalog of UV-Excess Galaxies (first paper, and extended list, see Markarian 1967;Markarian et al. 1989) includes the next classes of objects (with their abundance in brackets): galaxy (1, 251), starburst galaxy (94), Seyfert type 2 (65), galaxy HII region (26), LINER (17), QSO (13), Seyfert type 1.5 (12), BL Lac (3). As it can be seen from this list, only 16 objects out of 1469, 13 QSOs and 3 BL Lacs are classic radio-loud objects. Therefore, despite their closeness, one does not expect Markarian sources to be a dominant component of neutrino sources observable at radio frequencies, simply because the Markarian galaxies are typically radio-quiet objects (e.g. Biermann et al. 1980).
The intrinsic luminosity of a source with flux density S ν at observing frequency ν, spectral index α (with the convention S ν ∼ ν α ), seen at luminosity distance D L , redshift z, is: after the k-correction K(z) = (1 + z) −(1+α) . In Fig. 5 we show the number counts within our z-sample with ∆ log L = 0.25, ∆z = 0.25 binning. There is luminosity and redshift evolution in the sample, such that the more luminous sources were more abundant in the past. We found the most sources between log L(W Hz −1 ) ≈ 26 and log L(W Hz −1 ) ≈ 27.5, between the redshifts z ≈ 0.5 and z ≈ 2, such that the peak falls within the bin centered at log L p (W Hz −1 ) = 26.5, z p ≈ 0.75.
The space density of sources as function of the luminosity and redshift can be given in the form (e.g. Ajello et al. 2012) where φ(L, z) is the radio luminosity function (LF). In this equation dV /dz is the co-moving volume element per unit redshift and unit solid angle (e.g. Porciani & Madau 2001) which is in a flat universe: where c is the speed of the light, H 0 = 67.66 km Mpc −1 s −1 is the Hubble parameter at z = 0, ∆Ω s is the solid angle covered by the survey (= 4π for the all-sky CRATES), D L is the luminosity distance, Ω g = 5.4 × 10 −5 , Ω ν = 1.440 × 10 −3 , Ω m = 0.310, Ω Λ = 0.689 is the density parameter for the radiation, neutrinos, matter, and dark energy, respectively, taken from Planck Collaboration et al. (2020). First we estimated the LF by employing the 1/V a method of Avni & Bahcall (1980), which is the generalization of the 1/V max method of Schmidt (1968) for cases when sufficient number of sources is available in a given survey. If N object appears in the ∆ log L∆z bin element (L 1 < L < L 2 , z 1 < z < z 2 ), then the LF is where is the available volume of the ith source, when it is pushed to z i max , where the source would be on the edge of visibility given the lower flux limit of the survey (which is S min = 65 mJy at 4.8 GHz for CRATES). The LF φ(L, z) estimates the space density if all objects in a given population would be detectable by the survey with a flux limit. The dimension of φ(L, z) in our study is Mpc −3 / log L. The rms error estimate of the LF in each luminosity and redshift bin is calculated by weighting each sources by its contribution (e.g. Marshall 1985): However, if there are ten sources or less in a ∆ log L∆z bin, we used the tabulated upper and lower 84% confidence intervals given by Gehrels (1986), who calculated the Poissonian statistical asymmetrical errors when event rates are calculated from small numbers of observed events. These intervals correspond to the Gaussian 1σ errors such that σ φ = φ × σ N /N , where σ N is the Poissonian statistical asymmetrical error on the measured number of sources. The limiting luminosity distance D L,max of the ith source taking into account the CRATES flux limit at 4.8 GHz can be calculated by keeping L and setting S 4.8GHz = 65 mJy. Then the maximal redshift is calculated by solving for z max . Now we estimate the completeness of the z-sample derived from the CRATES catalog that comes from our redshift search method. We compared the observed number density to the number density predicted by the 1/V a method, which gives the space density if all sources could have been observed in the survey (with 65 mJy flux density limit). If the observed number density of sources in the jth ∆ log L∆z bin is n j (= N j ∆ log L −1 , ∆z −1 , the luminosity function is φ j (L, z), then The completeness parameter C(L, z) ∈ [0 ÷ 1] takes into account every geometrical and statistical correction factors (e.g. the effective area of the survey) that lead to N j,o = N j,p , where N j,o (N j,p ) is the observed (predicted) number of sources in the jth bin. In Fig 8 (see Appendix) we show the resulting completeness parameter C(L, z) in a grid with ∆ log L = 0.25, ∆z = 0.25. In Fig. 6 we show the luminosity of sources as function of the redshift, over which we capture all sources predicted by the 1/V a method and C(L, z) = 1. It seems at small redshift the completeness reaches 1 for smaller luminosities, while at larger redshifts we miss the sources with smaller luminosities. We fitted a linear-exponential model in the form of log L c 4.8 GHz (z) = az exp bz + c to set the completeness limit in the luminosity redshift plane. The resulting parameters are a = 1.97 ± 0.17, b = −0.309 ± 0.020, c = 25.248 ± 0.094 (χ 2 red = 0.67). Sources satisfying log L 4.8 GHz ≧ log L c 4.8 GHz (z) constitute our complete z-sample. The final, complete zsample at 4.8 GHz contains 1, 423 sources. Next we transformed the 8.4 GHz flux densities to 4.8 flux densities, such as then calculated the luminosities by using Eq. 6 and constrained the complete sample at 8.4 GHz based on Eq. 15. The complete z-sample at 8.4 GHz contains 1, 428 sources. We also calculated the 8.4 GHz luminosities.
We repeated the catalog-matching taking into account only the complete 4.8 GHz and 8.4 GHz z-samples derived from the CRATES catalog, and the 70 Ice-Cube track neutrinos from Giommi et al. (2020). We found that 13 neutrinos contain 14 CRATES sources at 4.8 GHz, and 14 neutrinos contain 15 CRATES sources from the z-sample. The average flux density of the neutrino-bright CRATES sources is 925 mJy (768 mJy) at 4.8 GHz (8.4 GHz) when we assume 100% signalness, and 956 mJy (762 mJy) at 4.8 GHz (8.4 GHz) when we assume 50% signalness of the neutrinos. The 50% signalness was taken into account as described in Section 3.

4.3.
Does the probability to detect a neutrino relate to the radio flux?
After constraining a sub-sample of the CRATES catalog complete in radio luminosities, we estimated that how many neutrinos could be explained with this subsample at 4.8 GHz and 8.4 GHz, if the probability to detect a neutrino is assumed to be proportional to the radio flux at these observational frequencies. With this we assume that neutrinos are emitted whenever a source is in an active stage, and assume that this probability scales with flux at radio frequencies. We also assume once a source is in an active stage, the probability to actually detect it runs with the k-corrected (radio) flux F = L/(4πD L 2 ). We assigned the sky coordinates of N CRATES sources from the 4.8 GHz complete z-sample to N neutrinos keeping the (randomly assigned) angular error of the neutrinos. The probability of the assigned direction depends on the radio flux of the given radio source. The right-ascension of the remaining 70 − N neutrinos was scrambled (with 0.1 decimal precision). We did the source finding as described in Section 2, and calculated the average flux density contained within the containment area of 70 neutrinos. We repeated the same process 10, 000 times, and then calculated the mean of the 10, 000 sample average flux densities and its 90% con- fidence interval. We then run N from N min = 1 to N max = 70. We then repeated the whole process employing the 8.4 GHz complete z-sample.
We plot the mean of the 4.8 GHz and 8.4 GHz sample average flux densities and their 90% confidence intervals in Figure 7, where we run N from 1 to 70. We see that higher the number of neutrinos with flux density-driven directions higher the sample average, as more and more neutrino directions equal the direction of the brightest CRATES sources and the catalog matching indeed finds these CRATES sources.
Taking into account the observed average flux densities (925 mJy at 4.8 GHz and 768 mJy at 8.4 GHz), we find that the CRATES catalog can explain 15.7 +12.6 −10.5 neutrinos at 4.8 GHz and 10.0 +11.3 −7 neutrinos at 8.4 GHz (90% C.L.), when the probability to detect a neutrino is assumed to be proportional to the radio flux. We repeated the test by assuming 50% astrophysical signalness of the neutrinos. In that case we found the CRATES catalog can explain 6.6 +11.9 −5.25 neutrinos at 4.8 GHz and 4.4 +10.3 −3.5 neutrinos at 8.4 GHz (90% C.L.).

SUMMARY OF RESULTS
5.1. Gamma-ray, X-ray and radio position matches with IceCube neutrino directions and probability of the ν-samples In Section 2. we found 29 Fermi-LAT 4FGL DR-2 gamma-ray sources (out of 5, 787), 61 Swift -XRT X-ray point sources (out of 9, 097, with SN R ≧ 10 and considering only AGN-type objects), and 87 (96) CRATES radio sources at 4.8 GHz (8.4 GHz) located within the 90% containment area of the 70 IceCube track-type neutrino events listed by Giommi et al. (2020).
We found that out of the 70 IceCube track-type neutrino events, 53 cover at least one gamma, or one Xray, or one radio source. Out of these 53 neutrinos, 15 contain at least one Fermi, one Swift -XRT and one CRATES source. One neutrino contains at least one Fermi and one Swift source. Five neutrinos contain at least one Fermi and one CRATES sources, and also five neutrinos contain at least one Swift and one CRATES sources. One neutrino contains only Fermi-LAT, 7 contain only Swift, and 18 neutrinos contain only CRATES sources.
With hypothesis tests we concluded the 70 IceCube track-type neutrinos select a special sample from the Swift -XRT 2SXPS catalog at a 1.2σ significance level, when we assume 50% signalness of the neutrinos. We found that there is a reasonable connection between the neutrinos and the Fermi-LAT 4FGL-DR2 gamma-ray sources (2.1σ), as well as between the neutrinos and the CRATES radio sources at 4.8 GHz (2σ), and at 8.4 GHz (2.1σ) assuming 50% signalness.

Test of radio flux dependence of neutrino directions based on complete z-samples
After constraining z-samples complete in 4.8 GHz and 8.4 GHz radio luminosities, we found that 13 neutrinos, out of the 70 IceCube track neutrinos from Giommi et al. (2020), contain 14 CRATES sources at 4.8 GHz, and 14 neutrinos contain 15 CRATES sources from the z-sample. Their observed average flux density is 925 mJy (768 mJy) at 4.8 GHz (8.4 GHz). We concluded that the complete z-samples derived from the CRATES catalog can explain 15.7 +12.6 −10.5 neutrinos at 4.8 GHz and 10.0 +11.3 −7 neutrinos at 8.4 GHz (90% C.L.), when the probability to detect a neutrino is assumed to be proportional to the radio flux. We found the complete z-samples derived from the CRATES catalog can explain 6.6 +11.9 −5.25 neutrinos at 4.8 GHz (between 4% and 53% of the neutrinos) and 4.4 +10.3 −3.5 (between 3% and 42% of the neutrinos) at 8.4 GHz (90% C.L.), when assuming 50% astrophysical signalness of the high-energy neutrinos.  (IC40, IC59, IC70, IC86-I, IC86-II, IC86-III, IC86-IV). They found the five maps with simulated source extension of 1 • , 2 • , 3 • , 4 • , 5 • , all are consistent with background only hypothesis. Hooper et al. (2019) studied one year of track data (IC86-2011) applying spatial match between this data set, the 3LAC and the Roma BZCAT catalogs. They also searched spatial and temporal match between the above neutrino data set and the Fermi Collaboration's All-Sky Variability Analysis. They found no evidence that blazars generate a significant flux of high-energy neutrinos, and found that no more than 5%-15% of the diffuse neutrino flux (IC86-2011) can originate from this class of objects. Smith et al. (2021) used three years of IceCube data (IC79, IC86-2011, IC86-2012 to search for evidence of neutrino emission from the AGN in the fourth catalog of AGNs detected by the Fermi-LAT (4LAC). They found at the 95% confidence level, that blazars can produce no more than 15% of IceCube's observed flux (2010)(2011)(2012), and that for non-blazar AGN it remains possible that this class of sources could produce the entirety of the diffuse neutrino flux observed by IceCube in three years.
We employed the full 4FGL-DR2 catalog (ten years of Fermi point-like sources) to search the sources of the 70 neutrinos. We found a ∼ 2σ connection between this catalog and the 70 track-type neutrino events. It has a similar meaning to the above findings: blazars, dominant population of the 4FGL catalogs, contribute to the IceCube neutrino flux but either the contribution is subdominant, or possible in-source gamma attenuation of neutrino sources complicates the picture.
What we can learn from these results is that there are indications for a correlation between the high-energy neutrinos measured with IceCube and the radio, X-ray and GeV gamma-ray catalogs under investigation in this paper. The question now is how to proceed in order to strengthen these potential correlations to move toward a 5σ confirmation. Knowledge on the sources and from other multimessenger channels needs to be exploited: what we need to understand is how the different messengers and energy channels are correlated. As an example, when looking at the broader picture of gammaray vs. neutrino emission, we are in this paper testing the correlation between GeV gamma-rays and neutrinos. However, it is clear from different observation channels (e.g. Kun et al. 2021) that the GeV gammaray intensity seems to be significantly lower than what is predicted in the typical correlated intensities of highenergy gamma-rays and neutrinos that arise from the decay of pions. This problem prevails both when looking at the diffuse signal (Ahlers & Halzen 2015), but also when looking at potential sources like NGC 1068 with a 3σ evidence for neutrino emission in the 10 year point source sample (IceCube Collaboration 2020). One solution to this problem is that the secondaries are produced in an environment of high photon or gas densities. In the case of NGC 1068, it is likely that the emission happens near the Corona of the accretion disk of this active galaxy (see e.g. Eichmann et al. 2021). As for blazars, a time-dependent modeling is certainly necessary in order to understand the complex multimessenger lightcurves, where the evidence of two very different flaring behaviors of TXS0506+056 is a first indication that no flare is like another (IceCube Collaboration et al. 2018;IceCube Collaboration 2018). Different AGN propagation codes have been developed in recent years in order to tackle this problem of time dependence, like ATHEνA (Dimitrakoudis et al. 2012), Böttcher (Böttcher et al. 2013), Paris (Cerruti et al. 2015), AM 3 (Gao et al. 2017), CRPropa (Hoerbe et al. 2020). Most importantly, progress in the understanding of the astrophysical properties of AGN is necessary to really be able to produce quantitative models. Knowing things like the gas and photon field distributions, magnetic field structure and strength from astrophysical observations help to tighten the parameter space and thus make the results of the modeling more reliable. Here, a crucial role lies in the understanding of the plasma, and acceleration and transport process of cosmic rays (Sironi & Spitkovsky 2011;Becker Tjus et al. 2022).

CONCLUSIONS
We examined if there is a correlation between the list of the 70 IceCube track neutrinos compiled by Giommi et al. (2020), the gamma-ray Fermi-LAT 4FGL-DR2, the X-ray Swift -XRT 2SXPS, and the radio CRATES catalogs. By collecting redshifts of sources in the CRATES catalog and deriving sub-samples complete in luminosities, we investigated whether the neutrino flux depends on the radio flux of the complete subsample. In the followings we list the most important takeaways from our recent work.
• We found similar levels of correlation between the 70 IceCube neutrinos and gamma-ray, X-ray and radio source catalogs. Each of these catalogs correlated with neutrino sources at the 1.2-2.1σ.
• The sub-sample of CRATES radio sources complete in radio luminosity can explain between 4% and 53% of the neutrinos at 4.8 GHz and between 3% and 42% of the neutrinos at 8.4 GHz (90% C.L. and 50% signalness), when the probability to detect a neutrino is assumed to be proportional to the (k-corrected) radio flux. This result is consistent with the contribution of AGNs to IceCube's neutrinos based on individually identified sources .
The next step is to disentangle the correlation between the three photon energies to understand which one is most relevant to neutrinos.