) Search for Spatial Correlations of Neutrinos with Ultra-high-energy Cosmic-Rays

For several decades, the origin of ultra-high-energy cosmic rays ( UHECRs ) has been an unsolved question of high-energy astrophysics. One approach for solving this puzzle is to correlate UHECRs with high-energy neutrinos, since neutrinos are a direct probe of hadronic interactions of cosmic rays and are not de ﬂ ected by magnetic ﬁ elds. In this paper, we present three different approaches for correlating the arrival directions of neutrinos with the arrival directions of UHECRs. The neutrino data are provided by the IceCube Neutrino Observatory and ANTARES, while the UHECR data with energies above ∼ 50 EeV are provided by the Pierre Auger Observatory and the Telescope Array. All experiments provide increased statistics and improved reconstructions with respect to our previous results reported in 2015. The ﬁ rst analysis uses a high-statistics neutrino sample optimized for point-source searches to search for excesses of neutrino clustering in the vicinity of UHECR directions. The second analysis searches for an excess of UHECRs in the direction of the highest-energy neutrinos. The third analysis searches for an excess of pairs of UHECRs and highest-energy neutrinos


(Published 2022 August 3 in The Astrophysical Journal )
ABSTRACT For several decades, the origin of ultra-high-energy cosmic rays (UHECRs) has been an unsolved question of high-energy astrophysics. One approach for solving this puzzle is to correlate UHECRs with high-energy neutrinos, since neutrinos are a direct probe of hadronic interactions of cosmic rays and are not deflected by magnetic fields. In this paper, we present three different approaches for correlating the arrival directions of neutrinos with the arrival directions of UHECRs. The neutrino data is provided by the IceCube Neutrino Observatory and ANTARES, while the UHECR data with energies above ∼50 EeV is provided by the Pierre Auger Observatory and the Telescope Array. All experiments provide increased statistics and improved reconstructions with respect to our previous results reported in 2015. The first analysis uses a high-statistics neutrino sample optimized for pointsource searches to search for excesses of neutrinos clustering in the vicinity of UHECR directions. The second analysis searches for an excess of UHECRs in the direction of the highest-energy neutrinos. The third analysis searches for an excess of pairs of UHECRs and highest-energy neutrinos on different angular scales. None of the analyses has found a significant excess, and previously reported overfluctuations are reduced in significance. Based on these results, we further constrain the neutrino flux spatially correlated with UHECRs.

INTRODUCTION
The Earth is continuously bombarded by high-energy cosmic rays, most of which are charged atomic nuclei (Particle Data Group et al. 2020). It is generally believed that cosmic rays with energies above 1 EeV (10 18 eV), known as ultra-high-energy cosmic rays (UHECRs), mostly originate from extra-galactic sources * Deceased in the nearby universe. Based on the estimated magnitude of galactic magnetic fields (Nagano & Watson 2000), cosmic rays below this energy are believed to diffuse within their host galaxy, whereas cosmic rays above this energy escape from the galaxy. These assumptions are confirmed by the observation of largescale anisotropies in the arrival directions of UHECRs above 8 EeV with the excess flux directed from outside of our Galaxy (Aab et al. 2017). The two largest observatories for UHECRs are the Pierre Auger Observatory (Auger) (The Pierre Auger Collaboration 2015) in Argentina in the Southern hemisphere and Telescope Array (TA) in the United States in the Northern hemisphere (Abu-Zayyad et al. 2012). They have performed detailed analyses of the arrival directions of detected UHECRs and have revealed interesting features, such as possible correlations with the directions of known starburst and active galaxies in the nearby universe observed by Auger (Aab et al. 2018;Abreu et al. 2021), and intermediate-scale directional clustering observed by TA (Abbasi et al. 2014;Tkachev et al. 2021). However, due to the small number of events detected at the highest energies, these indications of anisotropic arrival directions have not yet been confirmed on the 5σ level. Thus, neither Auger nor TA report an unambiguous identification of the UHECR sources to date. One of the main reasons is the deflection of cosmic rays by magnetic fields during propagation from their sources to Earth, which alters their directional information with respect to the source positions. The deflection of UHECRs increases with increasing charge, which is not well determined at the highest energies above 50 EeV considered in this work. This uncertainty additionally complicates the identification of UHECR sources. Nevertheless, UHECRs at the highest energies are deflected the least due to their extremely high magnetic rigidity (Alves Batista et al. 2019), which makes them suitable for directional correlation searches.
We search for the sources of the highest-energy cosmic rays with a multi-messenger approach using high-energy neutrinos and UHECRs with energies above ∼50 EeV. High-energy neutrinos are direct tracers of hadronic interactions of cosmic rays, and are thus expected to be produced at the acceleration site or during propagation (see e.g. Murase 2015;Meszaros 2018;Anchordoqui et al. 2008). Furthermore, the electrically-neutral neutrinos carry the directional information of their origin, since they are not deflected by magnetic fields. A combined analysis of UHECRs and high-energy neutrinos is further motivated by the observation of a diffuse flux of high-energy neutrinos of astrophysical origin above the background of atmospheric neutrinos. The flux follows a hard power law that extends to energies of multiple PeV and possibly beyond (The IceCube Collaboration 2013; Aartsen et al. 2016;Stettner 2019;Aartsen et al. 2020a;Abbasi et al. 2021). The measured flux is closely below the Waxman-Bahcall bound (Waxman & Bahcall 1999;Bahcall & Waxman 2001), which is a theoretical upper bound on the astrophysical neutrino flux derived from the observed UHECR flux, suggesting a connection between UHECRs and high-energy neutrinos (Murase et al. 2013;Ahlers & Halzen 2018). We note that the Seyfert and Starburst Galaxy NGC 1068 appears both as a neutrino source candidate in Aartsen et al. (2020b) and as a UHECR source candidate in Aab et al. (2018). Recently, the potential of high-energy neutrinos in a multi-messenger approach has been demonstrated together with observations of high-energy photons (Aartsen et al. 2018a,b), finding compelling evidence for neutrino emission from the blazar TXS 0506+056. At the same time, it can be concluded that blazars seem insufficient to saturate the total observed neutrino flux (Murase & Waxman 2016;Aartsen et al. 2017a) based on the non-observation of steady neutrino fluxes from known blazars in the universe. This result motivates different approaches to identify hadronic sources in the local universe; a promising method is to correlate the sources traced with UHECRs with high-energy neutrinos, and vice versa.
In this paper, we report results from three conceptually different approaches. All are based on correlating the arrival directions of high-energy neutrinos detected by the IceCube Neutrino Observatory (IC) (Aartsen et al. 2017b) and the ANTARES experiment (ANT) (Ageron et al. 2011) with the arrival directions of UHE-CRs with energies above ∼50 EeV measured by Auger and TA: • Analysis 1 (described in subsection 4.1) uses the measured UHECR directions as well as magnetic deflection estimates for identifying regions in which we search for point-like neutrino sources.
• Analysis 2 (described in subsection 4.2) uses the arrival directions of neutrinos with a high probability of astrophysical origin to search for a correlated clustering of UHECR arrival directions, while also accounting for the magnetic deflection.
• Analysis 3 (described in subsection 4.3) is based on a largely model-independent two-point correlation analysis of the arrival directions of UHECR and neutrinos with a high probability of astrophysical origin.
With these three approaches, we follow up on previous searches performed by the Pierre Auger, Telescope Array and IceCube collaborations (The IceCube Collaboration et al. 2016), which showed an interesting correlation at a significance level close to three standard deviations.
Since 2015, the analyses have been improved with substantially enlarged data sets, extending the Auger data from 10 yr to 13 yr (see subsection 2.3) and TA data from 6 yr to 9 yr (see subsection 2.4). The different data sets from IceCube are enlarged from 4 yr to 10.5 yr, from 4 yr to 7.5 yr, and from 2 yr to 8 yr, depending on the specific detection channel as described in subsection 2.1. Newly included are two data sets of 9 yr and 11 yr from ANTARES, as further specified in subsection 2.2. Due to the updated data sets by Auger and TA, their respective shift of the energy scale is updated from a sole shift of TA energies by −13% to a shift of +14% for Auger and −14% TA energies above the energy threshold of ∼50 EeV ). Furthermore, we include an improved magnetic deflection model that distinguishes between the Northern and Southern hemisphere for analysis 2. We report the results from the three improved correlation searches, which update the preliminary reported results in Schumacher (2019); Aublin et al. (2019); Barbano (2019). In addition, we report upper limits on the correlated fluxes of UHECRs and neutrinos based on benchmark models for the magnetic deflections.

OBSERVATORIES AND DATA SETS
All data sets used in this paper are used in previous work by the four respective collaborations. This section focuses on the main aspects relevant for our analyses.

The IceCube Neutrino Observatory
The IceCube Neutrino Observatory (Aartsen et al. 2017b) is an ice-Cherenkov detector sensitive to neutrinos with energies ≥ 5 GeV. It is located at the geographic South Pole, about 1.45 km to 2.45 km deep in the ice. Its main component consists of a volume of about 1 km 3 glacial ice instrumented with 5160 photomultipliers which are connected to the surface by 86 cable strings.
Two classes of neutrino-induced events can be phenomenologically distinguished: elongated, track-like events are produced by muons that originate mostly from charged-current ν µ interactions; and the spherical, cascade-like events that originate from charged-current ν e and ν τ interactions with hadronic and electromagnetic decays, as well as neutral-current interactions of all flavors. Typically, track-like events enable a better angular resolution than cascade-like events due to their different topologies, but provide a poorer energy resolution (Aartsen et al. 2014a;Wandkowsky 2018). One method of suppressing the dominant background of down-going muons produced by cosmic-ray interactions in the atmosphere is by selecting events with the interaction vertex within the detector (Aartsen et al. 2014b;Kopper 2017;Wandkowsky 2018). Alternatively, through-going tracks with either horizontal or up-going directions are selected, such that the atmospheric muons are blocked out by Earth (Aartsen et al. 2016;Haack & Wiebusch 2017). In the case of down-going tracks, a high energy threshold together with elaborate selection procedures are necessary to filter out atmospheric muons (Aartsen et al. 2018b(Aartsen et al. , 2017c. In all cases, the remaining event rate is dominated by atmospheric neutrinos. The selection of astrophysical neutrinos can be achieved on a statistical basis by selecting very energetic events or, in case of the very down-going region, by vetoing events where an atmospheric shower is observed in IceTop, Ice-Cube's surface detector for cosmic rays (Abbasi et al. 2013).
For the three analyses, data from multiple detection channels are used, which are: (i) a data set of throughgoing tracks from the full sky optimized for point-source searches (PS), (ii) a data set of high-energy starting events (HESE) of both topologies from the full sky, (iii) a data set of high-energy neutrinos (HENU) selected from through-going tracks with horizontal and up-going directions, and (iv) a data set of tracks from a selection of extremely high-energy events (EHE). The PS data set is used for analysis 1, while the HESE, HENU and EHE data sets are used for analyses 2 and 3. For analyses 2 and 3, track-like events from the HESE, HENU, and EHE data sets are combined, while multiple instances of identical events are removed. This results in a data set of 81 track-like events. In analyses 2 and 3, the 76 cascade-like events from the HESE data set are analyzed separately due to their larger directional uncertainty. The sky distribution of selected events is shown in Figure 1 and an overview of the nomenclature is presented in Table 1.
The PS data set consists of a combination of data collected from 7 years of operation between 2008 and 2015 that was used for point-source searches (Aartsen et al. 2017d) and data from 3.5 years of operation between 2015 and 2018 that was selected for the real-time gamma-ray follow-up (GFU) program of IceCube (Aartsen et al. 2018b(Aartsen et al. , 2017c. The combined data set consists of about 1.4 million track-like events above ∼ 100 GeV. It is dominated by atmospheric neutrinos in the Northern hemisphere and by atmospheric muons in the Southern hemisphere. The median of the angular resolution (Ψ) is better than 0.5°above energies of a few TeV. Figure 2 shows the median of the angular resolution for the different detector configurations of IceCube, and for data from the ANTARES detector (see subsection 2.2). Note that analysis 1 is not based on the median of the angular resolution, but on the estimator for the angular resolution on an event-by-event basis (σ).
The (2018). It consists of neutrinos of all flavors that interacted inside the detection volume, called starting events, with deposited energies ranging between about 20 TeV and 2 PeV. Integrated above 60 TeV, which corresponds to 60 events in total, the percentage of events of astrophysical origin, i.e. the astrophysical purity, is larger than 75 % (Wandkowsky 2018), while the percentage is lower below 60 TeV. The angular resolution is about 1°for track-like events and 15°for cascade-like events above 100 TeV. The resolution of the deposited energy for tracks and cascades is around 10 % (Aartsen et al. 2014a) without accounting for systematic uncertainties, but the cascades have a better correlation with the primary neutrino energy since they deposit most of their energy inside the detector, while tracks do not. The HENU data set consists of the 35 highest energy track-like events with a reconstructed declination ≥ −5 • , which have been collected between 2009 and 2016 (Haack & Wiebusch 2017) for the measurement of the diffuse muon-neutrino flux (Aartsen et al. 2016). From the original data set starting at ∼ 100 GeV, only events of high probability of non-atmospheric origin have been selected by applying an energy threshold of ≥ 200 TeV on the reconstructed muon energy. This corresponds to an astrophysical purity of more than 50 %. The astrophysical purity here is defined as the flux ratio of the astrophysical to the sum of atmospheric and astrophysical differential fluxes and thus depends on the assumed astrophysical flux. For the estimation of the astrophysical purity, the best fit astrophysical flux from Aartsen et al. (2016) dφ/dE = 1.01 × 10 −18 GeV −1 cm −2 s −1 sr −1 (E/100 TeV) −2. 19 is used as well as the best fit atmospheric flux from the same reference. The EHE data set consists of 20 events that have been collected between 2008 and 2017 (Aartsen et al. 2018c(Aartsen et al. , 2017c. The selection is targeting high-energy track-like events of good angular resolution ≤ 1°. The selection has been optimized to be sensitive to events in the energy range of 0.5 PeV to 10.0 PeV. The integrated astrophysical purity depends on the assumption on the spectrum of astrophysical events at the highest energies, but can be estimated as approximately 60 % purity (see Table 2 in Aartsen et al. (2017c)).

The ANTARES Neutrino Telescope
The ANTARES telescope (Ageron et al. 2011) is a deep-sea Cherenkov neutrino detector located 40 km offshore from Toulon, France, in the Mediterranean Sea. The detector is composed of 12 vertical strings anchored at the sea floor at a depth of 2475 m. The strings are spaced at distances of about 70 m from each other, instrumenting a total volume of ∼0.01 km 3 . Each string is equipped with 25 storeys of 3 optical modules (ANTARES Collaboration et al. 2002), vertically spaced by 14.5 m, for a total of 885 optical modules. Each optical module houses a 10-inch photomultiplier tube facing 45 • downward to optimize the detection of light from upward going charged particles. The detector was completed in 2008.
The ANTARES and IceCube detection principles are very similar (see subsection 2.1). Particles above the Cherenkov threshold induce coherent radiation emitted in a cone with a characteristic angle of 42 • in water. Median of the angular resolution, Ψ, which is the angle between the true neutrino direction and reconstructed muon direction for IceCube and ANTARES point-source data as a function of the true neutrino energy. This is calculated based on simulation data sets. The data from partial detector configurations of IceCube are denoted by the number of operating strings, (IC40, IC59, IC79; Aartsen et al. 2014c) and are shown in blue. The data of the full detector configuration with 86 strings including the GFU data set (Aartsen et al. 2017c(Aartsen et al. ,d, 2018b are identified by the years of data taking and are shown in orange. The angular resolution of ANTARES data for the 11-year data set (Albert et al. 2018;Illuminati et al. 2019) is shown in black.

Detector
Analysis 1  The position, time, and collected charge of the signals induced in the photomultiplier tubes by detected photons are used to reconstruct the direction and energy of events induced by neutrino interactions and atmospheric muons. Trigger conditions based on combinations of local coincidences are applied to identify signals due to physics events over the environmental light background due to 40 K decays and bioluminescence (Aguilar et al. 2007). ANTARES is thus also capable of detecting neutrino charged-and neutral-current interactions of all flavors.
For analysis 1, we use all track-like events from the 11-year data set used for point-source searches (PS) recorded between 2007 and 2017 (Albert et al. 2018;Illuminati et al. 2019). The high-energy events (HENU) for analyses 2 and 3 are selected from an earlier data set of track-like and cascade-like events collected between 2007 and 2015 (Albert et al. 2017). In order to ensure a high probability of astrophysical origin for analyses 2 and 3, we require an astrophysical purity ≥ 40 % based on the same definition as used for the HENU data set of IceCube. This selection results in a total of three tracklike events and no cascade-like events, of which the arrival directions are shown in Figure 1. An overview of the nomenclature is presented in Table 1. The track-like events are combined with the respective IceCube data sets. All events have a good angular resolution that is below 0.4°above 10 TeV (Albert et al. 2017). Overall, the angular resolution of ANTARES tracks is better than for IceCube tracks for energies below 100 TeV, and comparable around 100 TeV and above. The median angular resolution, Ψ, of the 11-year data set compared to the IceCube PS data set is shown in Figure 2.
Despite the smaller detection volume, the inclusion of ANTARES data significantly improves the all-sky coverage of the neutrino data set used in analysis 1. This is shown in subsection 4.1 and Figure 5 through the relative contribution to the expected number of signal events for the individual PS data sets of ANTARES and IceCube. Particularly in the Southern hemisphere, where the background from atmospheric muons results in a higher energy threshold in IceCube, the ANTARES data contributes substantially to the signal acceptance for soft source spectra.

The Pierre Auger Observatory
The Pierre Auger Observatory is located in Argentina (32°S, 69°W) at a mean altitude of 1.4 km above sea level (The Pierre Auger Collaboration 2015). The observatory has a hybrid design combining an array of particle detectors at ground (surface detector, SD) and an atmospheric fluorescence detector (FD) for detecting the air showers caused by UHECRs interacting with the atmosphere. The SD array is composed of 1660 water-Cherenkov detectors spread over an area of 3000 km 2 . The SD area is overlooked by the FD, which consists of 27 wide-angle optical telescopes located at four sites. The reconstruction of SD events is described in detail in . The energy estimate is based on the signal at 1000 m from the reconstructed intersection of the shower axis with the ground, which is extracted with a fit of the lateral distribution of signals in the individual detectors. This value is then corrected to take into account the different absorption suffered by showers coming at different angles. Due to the hybrid design of the observatory, this energy estimator is calibrated via the correlation with the near-calorimetric energy measured by the FD with the events observed by both SD and FD. Through this calibration, the energy estimate for SD is done without relying on Monte Carlo. At the energies considered in this work, the systematic uncertainty on the energy scale is 14 % (Dawson 2019), the statistical uncertainty on the energy due to the number of triggered SD stations and uncertainties of their response is smaller than 12 % (Abreu et al. 2011;, and the angular uncertainty is less than 0.9°( Bonifazi & The Pierre Auger Collaboration 2009). The UHECR data set used here consists of 324 events observed with the SD between January 2004 and April 2017 (Aab et al. 2018) with reconstructed energies ≥ 52 EeV and zenith angles ≤ 80°. The data statistics is enlarged with respect to the previously used data set (Aab et al. 2015) by 93 events and includes an updated energy calibration based on a larger number of hybrid events and an improved absorption correction. The angular acceptance translates into a field of view ranging from −90°to 45°in declination. All energies are scaled up by a constant factor of 14 % to match both Auger and TA energies on a common energy scale (see subsection 2.4), following the recommendation of the Auger-TA joint working group ). This scaling factor was determined by cross-calibrating the measured UHECR fluxes in the declination band around the celestial equator covered by both observatories. It is chosen such that the fluxes in the common declination band match at 40 EeV of the Auger data set and at 53 EeV of the TA data set, respectively. The celestial distribution of the UHECR arrival directions is shown in Figure 1. The directional exposure as function of declination can be found in Figure 3, which amounts to an integrated geometric exposure of 91 300 km 2 yr sr.  The Telescope Array (TA) is located in Millard County, Utah, USA (39.3°N, 112.9°W) at an altitude of about 1.4 km (Kawai et al. 2008). It consists of a surface detector (SD) array, composed of 507 plastic scintillation detectors of 3 m 2 each. The SD stations are located on a square grid with 1.2 km separation, which extends over an area of 700 km 2 (Abu-Zayyad et al. 2012). The TA SD array observes UHECRs with a duty cycle near 100 %. With its wide field of view, the SD array covers a range from −15°to 90°in declination. In addition to the SD, there are three fluorescence telescope stations, instrumented with 12 to 14 telescopes each (Tokuno et al. 2012). The telescope stations observe the sky above the SD array and measure the longitudinal development of the air showers as they traverse the atmosphere.
The data set for this analysis follows the selection in Abbasi et al. (2014), which has been updated to 9 years of data in Abbasi et al. (2018). The data set is identical to the data set used for anisotropy analyses presented in Troitsky et al. (2017). It consists of 143 events observed with the SD between May 2008 and May 2017 with reconstructed energies ≥ 57 EeV and zenith angles ≤ 55°. At these energies, the angular uncertainty is about 1.5°. The statistical uncertainty on the reconstructed energy is about 15 % to 20 %, with an additional systematic uncertainty on the energy scale of 21 % (Abbasi et al. 2018). All energies are scaled down by −14 % to match both Auger and TA energies on a common energy scale (see subsection 2.3), following the recommendation of the Auger-TA joint working group .
The celestial distribution of selected events is shown in Figure 1. The right-hand side of the figure is a histogram of the sine of declination of all UHECR events, showing that TA data substantially contributes to the full-sky exposure of the combined UHECR data set. Figure 3 shows the directional exposure as a function of declination, with an integrated total exposure of 11 600 km 2 yr sr ).

MAGNETIC DEFLECTIONS
Despite their extremely high rigidity, UHECRs are deflected in galactic and in extragalactic magnetic fields (GMF/EGMF) by a non-negligible amount. Neither the strength and correlation length of the extragalactic magnetic fields (Durrer & Neronov 2013;Kronberg 1994) nor the distance of the UHECR sources are known well. Measurements of the Faraday rotation of extragalactic sources indicate that the extragalactic magnetic fields are weaker than 1 nG (Pshirkov et al. 2016). Assuming a correlation length of ∼1 Mpc, this results in a deflection less than 2°for protons of 100 EeV, even at source distances of 50 Mpc. In line with the previously reported results (The IceCube Collaboration et al. 2016), the deflection outside of our galaxy is assumed to be generally weaker than within our galaxy, and it is not modeled explicitly but benchmarked within the uncertainty of the deflection by GMF and the uncertainty of the rigidity due to the unknown composition.
The measurement and modeling of the GMF is a complex task and subject of on-going discussion 1 . Among different proposed models, we use the JF2012 (Jansson & Farrar 2012) model and the PT2011 (Pshirkov et al. 2011) model to estimate the deflection of the UHE-CRs in the GMF. Both models consist of a disc and a halo component, while the JF2012 model has an additional x-shaped field component perpendicular to the disc. A propagation simulation has been conducted with a Monte-Carlo approach to estimate the deflection of protons with an energy of 100 EeV that are distributed isotropically outside of the GMF. The resulting deflections are shown in Figure 4, split into the Galactic Northern and Southern hemisphere based on the arrival direction of the proton. The estimated deflection shows a complex structure which differs considerably for the two models. However, the mean deflection is about 3°i n both cases. Due to the heavy tails of the distributions, the mean is consistently larger than the median, thus we choose the mean as a conservative deflection estimate. The split into the two hemispheres shows a considerably smaller deflection in the Galactic North than in the Galactic South due to North-South asymmetries present in both models of the GMF. For this work we have chosen a robust benchmark modeling of the effect of the GMF, thus avoiding biases of detailed model uncertainties. The deflection process is assumed to be random, resulting in a symmetric 2D Gaussian distribution of UHECR arrival directions around a given source direction. In reverse, the source position is assumed to be located within the 2D Gaussian distribution around the arrival direction of a UHECR event. The standard deviation of the Gaussian deflection, σ MD , depends on a scaling factor, C, and inversely on the energy of the UHECR event, E CR The default deflection, D 0 , for C = 1 and E CR = 100 EeV is derived from the mean of the deflection values obtained with the simulation, which are shown in Figure 4. Note that the deflection is usually a 2D quantity with (x, y) coordinates, which in Figure 4 is projected to the absolute value, i.e. x 2 + y 2 . Furthermore, larger and thus more conservative values of the deflection will be tested to include an uncertainty of the GMF model. The scaling factor C thus accounts for uncertainties of both GMF model and UHECR charge, which is not known on an event-by-event basis at highest energies. Analysis 1 uses a default deflection of D 0 = 3°over the whole sky, while analysis 2 uses D North/South = (2.4°, 3.7°) for UHECRs with arrival directions from the Galactic Northern and Southern hemisphere, respectively. The different choices have been made based  on the respectively better sensitivity for the two analyses. Analysis 1 and 2 implement the deflection into their methods as described in the respective sections 4.1 and 4.2. Analysis 3 employs a model-independent approach with respect to the UHECR deflection.

Unbinned neutrino point-source search with UHECR information
The goal of analysis 1 is to find point-like neutrino sources that are spatially correlated with UHECR arrival directions within a region derived from their magnetic deflection estimate. The search for neutrino sources utilizes the unbinned maximum-likelihood analysis commonly used in IceCube (Aartsen et al. 2017d(Aartsen et al. , 2019(Aartsen et al. , 2020b. This enables an easy combination of the high-statistics, full-sky neutrino data sets, which are the PS data sets of IceCube and ANTARES described in subsection 2.1 and subsection 2.2, respectively. In addition to this standard method, the magnetic deflection regions, defined by the UHECR arrival directions, energy, and scaling factor, are used for constraining the possible source regions. The unbinned neutrino likelihood in the source direction Ω = (α, δ) in right ascension and declination consists of the sum of a signal PDF, S, and a background PDF, B The likelihood product with index i runs over all neutrino events within each data set with index j (see subsection 2.1 and subsection 2.2). The likelihood product with index j runs over all data sets to yield the final likelihood function. The likelihood combination of the seven data sets by ANTARES and IceCube is thus handled in the same way as described in Aartsen et al. (2017d). The signal and background PDFs, S and B, are evaluated for four observables per neutrino event, weighted with the number of signal events over total number of events per data set, n s /N j . The observables are the reconstructed right ascension and declination, summarized as Ω i = (α i , δ i ), reconstructed energy, E i , and angular error estimator, σ i . The signal PDF consists of two terms: the first term is a declination-dependent reconstructed energy distribution, where the underlying neutrino flux is modeled as a power law Here, Φ ν 0 is the flux normalization, 1 GeV is the corresponding pivot energy, and γ is the spectral index of the power law. The second term of the signal PDF is a spatial term modeled as a Gaussian with the width, σ i , given by the angular error estimator of each neutrino candidate on an event-by-event basis. The background PDF is determined as function of reconstructed energy, E i , and declination, δ i , from randomized experimental data. A full description of the signal and background PDFs, S and B, can be found in Aartsen et al. (2017d). The proportionality factor, f j (sin δ, γ), is the relative signal acceptance per neutrino data set calculated from the expected number of signal events via This factor is used to correctly weight the signal contribution of each data set j for a given source declination and spectral index. Per data set j, the livetime is denoted with T j and the effective area is denoted with A j eff . Figure 5 shows the proportionality factor of each data set as a function of declination for two spectral indices, γ = 2.0 and 2.5. Note that all data sets are evaluated with the same formulation of the likelihood function as used by IceCube, instead of combining different formulations as described in Albert et al. (2020).
The best-fit signal parameters,n s andγ, at a given source position, Ω, are obtained with the maximumlikelihood method. The number of events and the proportionality factor are related to the neutrino flux using the respective livetime and effective area of each data set via The corresponding significance of a source at position Ω is evaluated using a likelihood ratio test, which yields the test statistic (TS) TS ν ( Ω) = 2 log L (n s ,γ) L(n s = 0) Ω .
Here, the null hypothesis is defined via n s = 0, representing the case of no neutrino sources in the vicinity of UHECR events. Instead of searching for a single neutrino source, the whole sky is searched on a HEALpix (Górski et al. 2005) grid with a resolution of approximately 0.2°. This results in a skymap of the TS values at each grid center, Ω grid . The second step is to combine the neutrino likelihood with the information provided by UHECR events and their deflection estimate. The deflection estimate for one UHECR with index k and arrival direction Ω CR k , which possibly originated in the direction of Ω grid , is defined as a 2D Gaussian. Its width is the quadratic sum of the magnetic deflection estimate, σ MD of Equation 1, and the UHECR angular reconstruction error, σ 2 CR , which is 0.9°for Auger events and 1.5°for TA events The term P k acts as a spatial constraint, which is multiplied to the neutrino likelihood function defined in Equation 2 via L → L·P k . Effectively, the constraint is added as a logarithmic term to the neutrino TS defined for each grid point in Equation 6 via TS → TS + 2 · log(P k ). Finally, the maximum of the combined UHECR-neutrino TS is found at the best-fit grid point,ˆ Ω grid , via The maximum marks the neutrino-source candidate which is the most-likely counterpart to the respective UHECR event used to calculate the spatial constraint. The normalization in Equation 7 adds only a constant term to the TS, which can be omitted when calculating significances and p-values based on pseudo-experiments (see below). As a third, final step, this procedure is repeated for all UHECRs and all obtained TS values are added up yielding the final test statistic. This search strategy was first developed in the context of this pointsource correlation analysis and first outlined in sec. 4 of Schumacher (2019). It has already been applied also to other neutrino correlation searches with spatially extended source localization, namely the correlation with ANITA events (Aartsen et al. 2020c) and with gravitational wave events (Aartsen et al. 2020d). Note that the analysis described here is improved with respect to the previous search (The IceCube Collaboration et al. 2016, sec. 5), as it models the displacement of a point-like neutrino source with respect to the UHECR arrival direction based on the assumed magnetic deflection. Here, the position and flux of a point-like source are fit in the vicinity of the UHECR direction, while in the previous search, a spatially-extended neutrino emission around the UHECR direction was fit.
Six different signal hypotheses are tested, which are characterized by a lower cut on the UHECR energy, E CR > E cut with E cut ∈ {70, 85, 100} EeV, and the scaling factor of the deflection estimate, C ∈ {1, 2}. The lower energy cut is a threshold for selecting only the highest-energy UHECRs with the lowest deflection for this analysis. No analogous energy cut is applied to the neutrino data. The scaling factor C is a model parameter for scaling the baseline magnetic deflection, and it is not derived from the UHECR data. The baseline magnetic deflection is D 0 = 3°for all UHECRs, which is then converted into the spatial constraint of an individual UHECR event using Equation 1 and Equation 7. The corresponding sensitivity and 3σ-discovery potential are evaluated based on the normalization of the neutrino flux per source, Φ ν 0 (see Equation 3). The 3σ threshold is chosen, since the background expectation needs to be calculated based on pseudo-experiments, as described in the next paragraph. A 5σ-discovery potential is computationally too expensive to be calculated for the various hypotheses. Sensitivity is defined as the expected median upper limit at 90% confidence level on Φ ν 0 in case of a null measurement, while the discovery potential is defined as the median of Φ ν 0 required for a rejection of the background hypothesis with 3σ significance.
Both sensitivity and discovery potential are calculated based on data challenges, for which signal-like and background-like pseudo-experiments (PEs) are generated. Since the calculation is based on PEs, the constant term in Equation 7 and Equation 8, i.e. the normalization of the constraint term, can be omitted. For all signal and background PEs, the UHECR arrival directions are kept fixed. The background PEs are obtained by randomizing the experimental neutrino data in right ascension. The signal PEs reflect the six different signal hypotheses; one signal PE is based on experimental neutrino data randomized in right ascension, to which Monte-Carlo neutrino events, representing a neutrino source, are added. The location of this neutrino source is chosen randomly within the spatial deflection constraint of one UHECR event. This way, we mimic a neutrino source that is displaced with respect to the UHECR direction due to the UHECR deflection. In the baseline model of the signal PEs, all UHECR events have such an artificial neutrino source in their vicinity. For one PE, all neutrino sources have an E −2 powerlaw spectrum and the same flux normalization. Note that the UHECR energy cut and scaling factor used to generate the spatial constraints for the likelihood function are the same as for determining the neutrino source location.
We generate additional signal PEs by varying the fraction of UHECR events, f corr , with a neutrino source in their vicinity. This mimics signal models where not all UHECRs have a corresponding neutrino source in the vicinity of their arrival direction. Thus, the number of neutrino sources is determined by the roundedup product of the correlation fraction with the numbers of UHECRs, N src ν = f corr N CR , where N CR is the number of UHECRs passing the energy cut. We choose the correlation fraction from three discrete values, f corr ∈ {1, 0.5, 0.1}, where f corr = 1 represents the baseline model. Note that the correlation fraction is only used to choose the number of neutrino sources in the signal PEs, while in a real experimental data, it is not known a priori which UHECRs have a correlated neutrino source and which do not. For f corr < 1, only a number of N src ν < N CR of randomly selected UHECRs will have a correlated neutrino source in the signal PEs. Independent of the correlation fraction, the TS values of all UHECR constraints are added up, since in real experimental data it would not be known which UHECRs have a correlated neutrino source. The sensitivity and 3σ-discovery potential in terms of the flux normalization per neutrino source for all signal models and hypotheses are reported in Table 2, which are also shown in Figure 6. We see that a smaller correlation fraction increases the flux normalization per neutrino source for both sensitivity and 3σ-discovery potential. This is expected since a smaller correlation fraction corresponds to fewer neu-trino sources, which in turn need to have a higher flux normalization to be detected by the analysis.
It is noticeable that for a small correlation fraction of f corr = 0.1, the signal model with a larger scaling factor of C = 2 yields a better 3σ-discovery potential than the smaller scaling factor of C = 1. This is unexpected, as a smaller UHECR deflection and thus a smaller spatial constraint should yield more accuracy in detecting the corresponding neutrino sources. However, the amount of sky covered by the constraints does not grow uniformly with the scaling factor and numbers of UHECRs, since the UHECR spatial constraints can overlap. In the case of a small correlation fraction, a high neutrino flux per source is required to reach the 3σ-discovery potential due to the small number of sources present 2 . A closer study showed that in the case of few, but strong neutrino sources, it is beneficial to search a larger area of the sky for sources so that fewer of them are missed. This effect is thus caused by an interplay between deflection size, number of UHECRs and relatively high neutrino flux per source in this particular case.

Unbinned likelihood-stacking analysis of UHECRs and high-energy neutrinos
The second analysis is based on using the neutrinos with a high probability of astrophysical origin as markers for the location of the sources of UHECRs and neutrinos. The UHECR events are stacked using an unbinned likelihood analysis with the arrival directions of the highenergy neutrinos as the source positions. The signal hypothesis is defined by the number of UHECR events, which are clustered around the neutrino arrival directions, as well as the size of the magnetic deflection. The background hypothesis is defined by an isotropic flux of UHECRs. This approach is thus complementary to the point-source correlation analysis, where the UHECR events are used as source markers and an isotropic flux of neutrinos defines the background hypothesis. The logarithm of the likelihood function is defined as where the free parameter is the total number of UHECR signal events, n s . The sums run over all UHECR events in the respective data set, i.e. N Auger and N TA , where the total number of UHECR events is N CR = N Auger + N TA . In contrast to the likelihood function of analysis 1, the signal PDFs, S i Auger and S i TA , and the background PDFs, B i Auger and B i TA , are stacked for all high-energy neutrinos such that n s is a global parameter. The background PDFs per UHECR detector, B i det , are the normalized expectations for an isotropic UHECR flux as function of declination, which correspond to the normalized detector exposures (see Figure 3). The signal PDF per UHECR detector is defined as as a function of the UHECR arrival direction, Ω i , and energy, E i . The term R det is the relative exposure of each detector as a function of the declination per UHECR event, δ i . Figure 3 shows the absolute directional exposures, which are each normalized to 1 over the whole sky to obtain R det . The sum runs over all neutrino events, N ν , where the terms of S j ( Ω i , σ(E i )) are the spatial likelihood maps obtained from the neutrino directional reconstruction smeared with the UHECR uncertainty with width σ(E i ) (see Equation 7). Thus, these terms are PDFs representing the total uncertainty on the common source position of UHECR event i and neutrino event j, evaluated at the UHECR arrival direction, Ω i . Similar to analysis 1, three different deflections corresponding to scaling factors of C ∈ {1, 2, 3} are represented in the signal terms of the likelihood function. Again, C is a model parameter that scales the magnitude of the deflection and it is not derived from UHECR data. Depending on whether the UHECR arrival direction is in the Northern or Southern Galactic Hemisphere (see Figure 4), the baseline magnetic deflection is thus The final value of the deflection, σ MD , is then calculated based on the UHECR energy using Equation 1. The best-fit of the number of signal events,n s , is found with a maximum likelihood estimation. The resulting test statistic is defined as the likelihood ratio of the maximum likelihood over the background likelihood with n s = 0 TS = 2 log L(n s ) L(n s = 0) .
The significance is then determined based on the assumption that the background expectation is distributed Sensitivity, ν flux per source Number of CR events Figure 6. Sensitivity (left) and 3σ-discovery potential (right) of the neutrino flux per source as function of the UHECR energy cut in EeV, for the three correlation fractions and for the two scaling factors we used.
according to a χ 2 function with one degree of freedom, which has been verified using background-like PEs. The test statistic is calculated separately for all track-like and all cascade-like neutrino events, and separately for the three different deflections. The analysis approach is essentially the same as published in The IceCube Collaboration et al. (2016), except for the updated magnetic deflection values, which are split into the Galactic Northern and Southern Hemisphere as described in section 3.
The sensitivity and the 3σ-discovery potential in terms of n s are obtained with data challenges based on PEs. Here, the neutrino arrival directions are kept fixed per PE, while the UHECR arrival directions and energies are generated resembling the signal and background hypothesis. For the background PEs, all UHECR arrival directions are derived from an isotropic flux and thus according to the exposure of the respective UHECR experiments. The energies of the UHECR events are sampled from a power-law proportional to E −4.2 for the Auger events and to E −4.5 for the TA events, same as in The IceCube Collaboration et al. (2016). For the signal PEs, a number of n s UHECR arrival directions are distributed randomly based on their respective spatial signal terms in Equation 10, S j ( Ω i , σ(E i )). Note that all UHECRs have the same scaling factor, corresponding to the respective signal hypothesis as implemented in the likelihood function. The sensitivity and 3σ-discovery potential for the three different benchmark values of C are presented in Table 3, separately for the track-like and cascade-like neutrino events. We find that the sen-sitivity and 3σ-discovery potential using neutrino tracks requires much fewer UHECRs than when using neutrino cascades, which is expected due to the large differences in angular reconstruction uncertainty of the two event topologies.

Two-point angular correlation of UHECRs and high-energy neutrinos
The 2pt-correlation analysis relies on counting pairs of UHECRs and the high-energy neutrinos, where the angular separation between their arrival directions is within a maximum angular separation, α. The observed number of pairs within this radius, n obs (α), is compared to the mean number of pairs expected from pure background, i.e. uncorrelated arrival directions, n bckg (α) . The relative excess of pairs is defined as As the separation angle is not known a-priori, angles between 1°and 30°in steps of 1°are tested. The angle with the largest deviation from the background expectation is chosen after the scan. The significance of the experimental result is calculated with respect to background PEs. Similar to analysis 2, the background PEs are generated with a fixed set of neutrino arrival directions, and uncorrelated UHECR arrival directions are generated according to the exposure of the respective UHECR experiments. The energies are sampled from the same spectra as described in subsection 4.2. As a cross-check, additional background PEs are generated with the fixed set of UHECR arrival directions, and neutrino arrival directions are randomized in right ascension. The different types of background PEs approximate an isotropic flux of UHECR and high-energy neutrinos, respectively.
Note that this analysis does not include an assumption of the magnetic deflection of UHECRs, which makes it a robust and model-independent approach. The analysis approach is the same as published in The IceCube Collaboration et al. (2016).

Unbinned neutrino point-source search with UHECR information
The test statistic of experimental data is obtained with the neutrino point-source correlation analysis as described in subsection 4.1, assuming the six different combinations of signal parameters, i.e. combinations of scaling factor, C, and lower energy cut, E cut . Each of the six experimental TS values is evaluated with respect to the corresponding distribution of TS obtained from background-like PEs. The resulting p-values for all six tests are listed in Table 2. The smallest p-value (9.7%) is found for the scaling factor of C = 2 and an UHECR energy cut of E cut = 85 EeV. The p-value increases to 24% when correcting for the trials due to the six correlated tests. This correction is based on the combined distribution of all minimum p-values of the background PEs. All p-values are fully compatible with the background hypothesis of no resolved neutrino sources in spatial correlation with the UHECRs considering the assumed signal models. Based on the experimental TS value, 90% C.L. upper limits on the normalization of the flux of neutrino sources correlated with UHECR arrival directions are reported in the last three rows of Table 2. In the case of the single p-value larger than 50% found for C = 1 and E cut = 100 EeV), the limits are set equal to the sensitivity in order to not over-estimate the limits based on an under-fluctuation in data. In addition to the baseline correlation fraction of 100%, the limits are calculated for two smaller correlation fractions of 50% and 10%. In these cases, the limits are relaxed by about 40 % to 60 % for f corr = 0.5 compared to f corr = 1, and by about a factor of 3 to almost 5 for f corr = 0.1 compared to f corr = 1.

Unbinned likelihood-stacking analysis of UHECRs and high-energy neutrinos
The UHECR stacking analysis is performed for tracklike and cascade-like high-energy neutrinos separately for three different scaling factors of C ∈ {1, 2, 3}. This results in six p-values with respect to the background hypothesis of an isotropic flux of UHECRs, of which none is significant. The smallest p-value is 0.26, which is found for cascades and the largest scaling factor of C = 3, corresponding to the benchmark deflection in the Northern and Southern Galactic hemisphere of (7.2 • , 11.1 • ). Based on the observed TS values, 90% C.L. limits on the number of UHECR events correlated to the high-energy neutrinos are calculated using the PEs of the corresponding signal hypothesis for each of the six tests. Since all p-values using the track-like neutrinos are larger than 0.5, the limits are set equal to the sensitivity in order to not over-estimate the limits based on an under-fluctuation in data. All p-values and the corresponding limits are reported in Table 3 together with the sensitivity and discovery potential.

Two-point angular correlation of UHECRs and high-energy neutrinos
The 2pt-correlation analysis is applied to the sets of neutrino tracks and cascades separately, as well as for the background hypothesis of an isotropic UHECRs flux and an isotropic high-energy neutrino flux. The significance of the result with respect to the background hypothesis is corrected for the scan over the separation angles. This results in four p-values and four respective best-fit angular separation, as listed in Table 4. The largest deviation from an isotropic flux of high-energy neutrinos (p-value=15%) is found using cascades and a maximum angular separation of 16°. None of the pvalues show a significant result, thus the results are all compatible with their respective background hypothesis. The relative excesses of pairs with respect to an isotropic distribution of neutrinos as a function of the separation angle are shown in Figure 7, separately for neutrino tracks and cascades.

CONCLUSIONS
Three complementary analyses have been performed on the UHECR data sets provided by the Pierre Auger and the Telescope Array Collaborations combined with the high-energy and full-sky neutrino data sets provided by the IceCube and ANTARES Collaborations. For both the UHECR and neutrino data sets, the combination of data from the two respective observatories provides a field of view over the entire sky. None of the analyses found a result incompatible with the assumed background hypotheses of either an isotropic neutrino flux or an isotropic UHECR flux. This is an important update on the results reported in The IceCube Collaboration et al. (2016), where p-values close to the 3σ level were reported when applying analyses 2 and 3 to cascade-like neutrino events.  Table 2. Sensitivity, discovery potential and 90% C.L. upper limit for the different analysis parameters, which are the UHECR energy cut, Ecut, magnetic deflection, D0 · C, and correlation fraction, fcorr, for the point-source correlation analysis. The values are given as normalization of the neutrino flux per source in the unit of 10 −10 GeV −1 cm −2 s −1 . The neutrino flux per source is modeled with a power-law of the form dΦ/dE = Φ0 · (E/1 GeV) −2 . The last row states all six experimentally obtained pre-trial p-values with respect to the null hypothesis of an isotropic neutrino flux.    Table 4. Post-trial p-value and best-fit angular separation for the 2pt-correlation analysis obtained with the neutrino data sets of high-energy tracks and cascades, as stated in the first column. The p-values are corrected for choosing the largest deviation out of all maximum angular separations. The respective background hypotheses are stated in the second column.

Analysis parameters
Based on the results, 90% C.L. upper limits are calculated for analysis 1 on the flux of neutrinos from pointlike sources correlated with UHECRs, and for analysis 2 on the number of UHECR events correlated with high-energy neutrinos. Analysis 1 reports upper limits on the correlated neutrino flux per source given different sets of parameters that make assumptions about the lower cut on the UHECR energy and the fraction of UHECRs with a correlated neutrino source, as well as about a scaling factor of the magnetic deflection (see Table 2). Analysis 2 reports upper limits on the number of correlated UHECR events with either track-like or cascade-like neutrinos depending on the scaling factor of the magnetic deflection (see Table 3). This is the first time that upper limits from a direct correlation search of UHECRs and neutrinos are reported.
These limits are calculated based on the assumed signal and background hypotheses and are thus based on the parameters of the respective signal models. In the following, we discuss the limitations of these signal models in the light of underlying assumptions made for the neutrino-UHECR correlation. There are several plausible explanations why we have not observed a significant correlation of UHECRs and neutrinos.
Neutrino sources are presumably distributed over the whole universe, and the emitted neutrinos are able to reach Earth without deflection or absorption. For example, the first neutrino source identified with compelling evidence, TXS 0506+056 (Aartsen et al. 2018a,b), is located at a redshift of z = 0.34 (Paiano et al. 2018), corresponding to a comoving distance of about 1.3 Gpc. This is beyond the assumed horizon for UHECR sources in the local universe of up to about 200 Mpc. Therefore, the fraction of correlated sources within the UHECR horizon is small compared to the total number of neutrino sources. This was quantified by Palladino et al. (2020) for muon neutrinos above 200 TeV based on the non-observation of neutrino multiplets above this energy. They concluded that the high-energy flux measured with IceCube must come from a large number of sources such that no multiplets are observed. Necessarily, most of the neutrino sources lie beyond the UHECR horizon, which is used in Palladino et al. (2020) to explain the lack of UHECR-neutrino correlations found in The IceCube Collaboration et al. (2016) and Schumacher (2019), and the same applies to the current study. In these studies, however, the energy threshold for neutrinos lies around 20 TeV for the combined highenergy data sets and as low as 100 GeV for the PS data set. The non-observation of UHECR-neutrino correlations thus extends to even lower energies than considered by Palladino et al. (2020).
Even if there were local sources of both UHECRs and neutrinos, they could be transient phenomena, such that the UHECRs emitted over a short period of time arrive much later due to their deflection, whereas the neutrinos travel on a straight path. It has been estimated by Davoudifar (2011) that the deflection in the EGMF causes a typical time delay on the order of 10 5 years considering a source distance of 50 Mpc and an EGMF strength of 2 nG. A delay of a couple of decades as expected from propagation in the GMF is already sufficient to de-correlate the observed UHECRs and neutrinos, as the livetime of the neutrino and UHECR data sets ranges between 7 and 13 years.
Another source of uncertainty is the mass composition and thus the charge of UHECRs at the highest energies. The measurements of the UHECR composition above ∼50 EeV are not yet conclusive. However, several composition constraints at the highest en-ergies (Aab et al. 2014a,b) and the lack of significant, magnetically-induced structures in the UHECR arrival directions  suggest that less than 10% of UHECRs above ∼50 EeV might be light elements. Only the light UHECRs are expected to have their source in the vicinity of their arrival direction as derived from the benchmark deflection model. The sources of heavier UHECRs, e.g. iron nuclei, are spatially almost uncorrelated to the UHECR arrival directions due to the 26times larger magnetic deflection compared to protons. We see in analysis 1 that the limits on the neutrino flux are significantly relaxed when assuming a correlation fraction of 50% or 10%. This correlation fraction, defined as the number of UHECRs with a neutrino source in their vicinity, is a simplified model of the UHECR composition.
The modeling of the magnetic deflection as a 2D Gaussian as assumed for analysis 1 and 2 is a simplification with respect to the PT11 and JF12 GMF models (Pshirkov et al. 2011;Jansson & Farrar 2012), which themselves are subject to large uncertainties. Especially in the region of the Galactic Plane, we expect deflections that are larger than the assumed mean value of around 3°. In addition, coherent deflection of UHECRs in the GMF or the deflection of UHECRs in the IGMF are not explicitly accounted for. Overall, a coherent shift of UHECRs depending on their arrival direction or an overall significantly larger deflection in the GMF and IGMF can dilute the spatial correlation of UHECRs and neutrinos. Nevertheless, a scaling of the overall deflection is covered partially with the scaling factor applied in analyses 1 and 2, since the overall strength of the deflection and the UHECR charge are largely degenerate parameters.
From a theoretical perspective, the connection of UHECRs and high-energy neutrinos is plausible based on their similar energy budget (Waxman & Bahcall 1999;Bahcall & Waxman 2001;Murase et al. 2013;Ahlers & Halzen 2018). However, this connection could not be proven with the current data and analysis approaches. It is to note that neutrinos with energies in the TeV-PeV range originate most likely from cosmic rays with energies below ∼100 PeV (Murase & Waxman 2016), which is below the energy threshold of > 50 EeV of the UHECR data sets used here. A direct connection to UHECRs can only be proven with ultra-high-energy neutrinos in the EeV range that have not been discovered yet. The non-observation of a UHECR-neutrino correlation rather points to the possibility that efficient UHECR sources are less efficient neutrino sources and vice versa: Sources with efficient UHECR acceleration and emission require an optically thin proton and ra-diation environment, while sources with dense proton and radiation targets are efficient in neutrino production, but not in UHECR emission (Murase et al. 2014;Rodrigues et al. 2021). This argument, however, can be relaxed when considering models where cosmic rays below ∼100 PeV are confined in a calorimetric environment and subsequently produce TeV-PeV neutrinos, while a fraction of the cosmic rays are accelerated to the highest energies and escape the source before interacting (Murase & Waxman 2016;Ahlers & Halzen 2018). Although no indication of such a scenario has been found in this analysis, the first indication of such a connection could be the observation of the Seyfert Galaxy NGC 1068 as a potential neutrino source candidate (Aartsen et al. 2020b;Inoue et al. 2020) and UHECR source candidate (Aab et al. 2018). As such sources are numerous also within the UHECR horizon, dedicated future searches correlating UHECR, photon, and neutrino observations might be able to set constraints on specific source candidates, instead of relying on UHECR and neutrino data alone. NGC 1068 as a potential neutrino source candidate (Aartsen et al. 2020b;Inoue et al. 2020) and UHECR source candidate (Aab et al. 2018) could serve as blueprint for a catalog of potential common sources to be constrained with dedicated analyses.
In summary, the three analyses presented here reflect complementary approaches for tackling the question of a common origin of UHECRs and high-energy neutrinos. Despite substantially enlarged data sets and improved methods, the initially reported results in The IceCube Collaboration et al. (2016) could not be strengthened. For future analyses, we expect substantial gains in sensitivity if the charge of the UHECRs could be estimated on an event-by-event basis, as it is expected for future measurements by Auger Prime (Aab et al. 2016).