{'fraction_non_alphanumeric': 0.029978962131837306, 'fraction_numerical': 0.016693937977247934, 'mean_word_length': 4.589830139372823, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Search for correlations between the arrival directions of IceCube neutrino events and ultrahigh-energy cosmic rays detected by the Pierre Auger Observatory and the Telescope Array

This paper presents the results of different searches for correlations between very high-energy neutrino candidates detected by IceCube and the highest-energy cosmic rays measured by the Pierre Auger Observatory and the Telescope Array. We first consider samples of cascade neutrino events and of high-energy neutrino-induced muon tracks, which provided evidence for a neutrino flux of astrophysical origin, and study their cross-correlation with the ultrahigh-energy cosmic ray (UHECR) samples as a function of angular separation. We also study their possible directional correlations using a likelihood method stacking the neutrino arrival directions and adopting different assumptions on the size of the UHECR magnetic deflections. Finally, we perform another likelihood analysis stacking the UHECR directions and using a sample of through-going muon tracks optimized for neutrino point-source searches with sub-degree angular resolution. No indications of correlations at discovery level are obtained for any of the searches performed. The smallest of the p-values comes from the search for correlation between UHECRs with IceCube high-energy cascades, a result that should continue to be monitored.


Introduction
The origin of the ultrahigh-energy cosmic rays is one of the main open questions in high-energy astrophysics. It is likely that they are accelerated in astrophysical objects, since more exotic scenarios for their origin, such as the top-down models, are strongly constrained by stringent limits on the flux of primary photons at such energies [1,2]. However, their sources still remain undiscovered after five decades of experimental efforts that have culminated in the two largest UHECR experiments ever built, the Pierre Auger Observatory [3,4] and the Telescope Array [5], in the Southern and Northern hemispheres, respectively. The two experiments, which measure UHECRs through the detection of extensive air showers produced in the Earth's atmosphere, have collected in total more than 300 events with energies in excess of ∼ 50 EeV (where 1 EeV ≡ 10 18 eV). These energies are well within the region of very significant cosmic-ray flux suppression observed by both experiments [6,7]. If this suppression were due to the degradation of the cosmic-ray (CR) energies by their interactions with the cosmic radiation backgrounds [8], one would expect that the sources of such energetic CRs are relatively nearby, at distances shorter than ∼ 200 Mpc. However, the results of several studies of UHECR arrival direction distributions [9,10] have shown no significant small-scale anisotropy and have not allowed the establishment of statistically-significant correlations with nearby astrophysical sources.
The identification of UHECR sources by correlation analyses is a challenging task. Since UHECRs are mostly charged particles, their trajectories from their sources to the Earth are deviated by magnetic fields, most notably by the Galactic field. The deviations depend on the energy and charge of the cosmic ray: while at energies above 50 EeV they may be of order of a few degrees in the case of protons, such deviations would be much larger for nuclei of the same energy but having a larger charge. This situation is further complicated by the facts that there are still large uncertainties in the modeling of magnetic fields (see [11,12] and references therein) and that the UHECR mass composition at energies above 50 EeV is still largely unexplored [13,14].
Further, UHECRs accelerated in astrophysical sources are naturally expected to produce high-energy photons and neutrinos in interactions with the ambient matter and radiation. In contrast to UHECRs, photons and neutrinos are not subject to magnetic deflections in their propagation to the Earth and hence they point directly to their sources. However, high-energy photons may also be produced in purely leptonic processes, and hence may be unrelated to hadronic CRs. Moreover, photons at energies exceeding a few TeV can be significantly attenuated by e + e − -pair production in interactions with the Cosmic Microwave Background (CMB) and other radiation backgrounds as they propagate from faraway sources. On the other hand, neutrinos interact with matter very weakly, so that they travel essentially unattenuated, and hence are excellent tracers of the sources of cosmic rays. Due to the challenge of neutrino detection, neutrino astronomy has however become viable only very recently, thanks to the discovery by the IceCube neutrino telescope [15] of a flux of extraterrestrial high-energy neutrinos extending up to at least 2 PeV.
The observations of the highest-energy events, on the one hand in neutrinos by IceCube and on the other, UHECRs, by the Auger Observatory and the Telescope Array, set the stage for a joint analysis of the three Collaborations presented in this paper, namely the search for a possible association between the IceCube neutrinos and the UHECRs. Some analyses searching for possible correlations have been presented in [16][17][18], but using smaller data sets.
Attempting to exploit the IceCube neutrino observations to identify the sources of UHECRs is motivated by several reasons. First, the arrival directions of the IceCube high-energy sample, in which a cosmic component has been identified [15], is consistent with an isotropic distribution and therefore with the hypothesis of an extragalactic origin for the bulk of the events, although the presence of a Galactic component is possible. Second, the inferred diffuse flux of astrophysical neutrinos, with energies between 30 TeV and 2 PeV, is at the level of the Waxman-Bahcall flux [19]. This is the upper limit for the extragalactic neutrino flux calculated assuming that UHECRs are protons escaping from optically thin sources with an injection spectrum scaling as E −2 and normalized to the measured UHECR flux above EeV energies. Clearly, if the sources are not opticallythin, or if the injection spectrum differs from E −2 , or the UHECRs are not protons, the upper limit on the flux would be considerably modified [20].
If neutrinos result from the decays of pions produced in pγ or pp processes, they would carry about 3-5% of the proton energy. Hence, the neutrinos observed by IceCube with energies of 30 TeV to 2 PeV would have been produced in this case by protons with PeV-100 PeV energies. Although these energies are much smaller than those of the UHECRs considered here, the third argument for our analysis is based on the fact that it is only at the highest energies that CR trajectories could point to the sources, while at much lower energies CRs diffuse and may not even arrive to the Earth from faraway sources. Since UHECR sources will also produce lower energy CRs, it is possible that they are the same sources that produced the observed neutrinos. One should keep in mind, however, that there may also be sources able to produce the PeV neutrinos but which do not accelerate CRs up to ultrahigh-energies. We also note that electron antineutrinos produced in the β decay of neutrons (or of radioactive nuclei produced in photo-disintegration processes) would carry about 4 × 10 −4 of the energy of the parent nucleon, and hence PeV neutrinos could in principle also arise from EeV neutrons.
In this paper we present several searches for correlations between the arrival directions of UHECRs detected by the Auger Observatory and by the Telescope Array, and of different sets of events observed by the IceCube Neutrino Observatory. The main characteristics of the three detectors, and their respective data sets, are described in Section 2. The Auger and Telescope Array samples include 231 events with E ≥ 52 EeV and 87 with E ≥ 57 EeV, respectively. The IceCube data sets that we consider, obtained using different selection criteria, comprise 39 cascades (signatures of charged-current ν e interactions as well as neutral-current interactions of all flavors) and 7 high-energy tracks (signatures of charged-current ν µ interactions) of the so-called 'High Energy Starting Events' (HESE) [15,21,22], as well as 9 high-energy muon tracks whose extraterrestrial origin is highly probable [23]. We also exploit the so-called point-source sample of about 400,000 tracks with a sub-degree angular resolution [24].
The search for correlations between neutrinos and UHECRs is not straightforward. On the one hand, as noted above, due to the magnetic deviations discussed in Section 3, the UHECR arrival directions do not exactly correspond to the source positions. On the other hand, the neutrino HESE cascade events are characterized by a large angular uncertainty. Moreover, the IceCube samples are not exempt from backgrounds due to atmospheric secondary particles. In particular, the point-source sample is largely dominated by atmospheric muons in the Southern hemisphere and by atmospheric neutrinos in the Northern one. To account for these aspects, three different analyses have been devised. Two analyses use the HESE and high-energy muon samples (Section 4). The first one is a cross-correlation analysis where the number of pairs between CRs and the highenergy neutrino sample is evaluated within a large range of angular windows, from 1 • up to 30 • . The second one adopts a stacking likelihood analysis, summing the contributions from the different sources. In this case, it is the arrival directions of the neutrinos that are being stacked. Plausible UHECR magnetic deflections and the neutrino point-spread functions are accounted for, together with the possible backgrounds. A third analysis uses a likelihood method applied to the point-source sample (Section 5). We note that the details of the analyses (angular scales explored, data samples considered, null hypotheses) were selected, as discussed in Sections 4.1, 4.2 and 5.1, before unblinding the analyses (see also [25]), while some a posteriori cross checks are also mentioned. Finally, we conclude in Section 6.
2 The observatories and data sets

The IceCube Neutrino Observatory
The IceCube Observatory is a km 3 -sized Cherenkov detector embedded in the ice at the geographic South Pole [26]. Optimized to detect neutrinos above ∼ 100 GeV energies, it consists of 5160 photomultiplier tubes (PMTs) along 86 cables (called 'strings') instrumented between depths of 1450 and 2450 m beneath the surface of the ice sheet. Each PMT is housed in a digital optical module (DOM), consisting of a pressure-resistant sphere with on-board digitization and Light Emitting Diodes (LEDs) for calibration [27]. The DOMs detect Cherenkov photons emitted by charged leptons that traverse the ice volume [28]. Starting in the year 2005, the detector ran in partial configurations consisting of 9, 22, 40, 59 and 79 strings until its completion in December 2010.

High-energy starting events (HESE)
In 2013, IceCube reported the first evidence for a high-energy neutrino flux of extraterrestrial origin from a search carried out on data collected between May 2010 and May 2012 in the 79-string and full 86-string configurations of the detector [15]. That search, targeting high-energy neutrinos interacting within the IceCube detector, was later updated to include two more years of data taken with the full 86-string configuration, extending the observation period to May 2014 [21,22]. These searches have consistently demonstrated the existence of an astrophysical flux of neutrinos at the level of E 2 ν φ(E ν ) ∼ 10 −8 GeV cm −2 s −1 sr −1 per flavor, emerging above atmospheric neutrinos starting at around 100 TeV and extending to the PeV range. The hypothesis of a purely atmospheric explanation of the neutrino events collected during 4 years has been rejected at around 6.5σ. However, the 54 events consisting of 39 cascades and 15 track-like events have shown no significant directional clustering, leaving unanswered the question of their origin [21,22].
The median angular resolution of the events containing muon tracks is around 1 • , while the resolution of cascade-like events is around 15 • [15]. The expected atmospheric background for the HESE sample is predominantly track-like atmospheric muons not vetoed by external layers of DOMs and sneaking in through the region called the 'dust layer' (where dust obscures the transparency of ice) and from atmospheric ν µ chargedcurrent interactions from the Northern hemisphere. We have therefore not used all the detected track-like events, but only the 7 that are more likely to be of extraterrestrial origin due to their high energies and/or arrival directions. Atmospheric neutrinos from the Southern hemisphere are accompanied by muons from the same shower and this background is greatly reduced by the event selection using a muon veto [21].
Within this paper these 39 cascade events and 7 track-like events are referred to as 'high-energy cascades' and 'HESE tracks', respectively, and their properties are listed in Table 4 and Table 5 in the Appendix.

High-energy through-going muons
A data sample of 35,000 through-going muons from the Northern sky was selected from events recorded between May 2010 and May 2012 in a search for a diffuse up-going ν µ flux [23]. While the majority of these events originate from neutrinos produced by cosmic-ray interactions in the Earth's atmosphere, at the highest energies a clear excess was observed. The hypothesis that the excess is of atmospheric origin was rejected at the 3.7σ confidence level [23]. This excess is compatible with the sum of the predicted muon neutrinos of atmospheric origin and an astrophysical flux consistent with the flux resulting from the HESE analysis. This deviation from the expectation that all events are of atmospheric origin occurs at the highest energies, where 9 events have been detected. These events (which are available in [29]) are listed in Table 5 and have a median angular resolution of around 1 • . The combination of these events with the 7 selected HESE tracks is referred to as 'high-energy tracks'.

4-year point-source sample
High-energy muons have long tracks within IceCube and they can be reconstructed with a median angular resolution smaller than 1 • , and hence they have good sensitivity for point-source searches. The corresponding effective volume to neutrino-induced muon events can be larger than the instrumented volume since the interaction vertex of the neutrino does not have to be contained within the detector. Four dedicated samples of such events have been isolated in three years of operation of IceCube in the incomplete 40-, 59-and 79-string configurations [30,31] as well as in one year of operation of the completed 86-string detector [24]. Note that 7 events out of the 9 high-energy throughgoing muons described in Section 2.1.2 are contained in the point-source data set.
The angular resolution of this sample has improved over the years, benefiting from the larger instrumented volumes achieved as the number of strings was increasing as well as from improvements in the reconstruction methods (see Figure 1). These samples are dominated by muons from charged-current interactions of atmospheric ν µ (andν µ ) in the Northern hemisphere (up-going data) and muons produced by interactions of cosmic rays with the Earth's atmosphere in the Southern hemisphere (down-going data). In the Northern sky, the effective area of the detector, determined by the analysis cuts and the opacity of the Earth for neutrinos with energies above ∼ 100 TeV, leads to a sample of events peaked in the 1 TeV -1 PeV range. In the Southern sky the sensitivity is in the 100 TeV -100 PeV range (see Figure 2) due to the overwhelming background of atmospheric muons that needs to be reduced by imposing a higher threshold cut on the energy proxy.

The Pierre Auger Observatory and its data sample
The Pierre Auger Observatory [3,4] (Malargüe, Argentina, 35.2 • S, 69.5 • W, 1400 m a.s.l.) combines a large surface detector (SD) with an air-fluorescence detector (FD) to observe, in a complementary way, extensive air showers generated by ultrahigh-energy cosmic rays. The SD is composed of an array of 1660 water-Cherenkov stations spread over an area of about 3000 km 2 . It measures, with a duty cycle of nearly 100%, the particles (mainly muons, electrons and photons) reaching ground level. The FD, comprising 27 telescopes at four sites, overlooks the surface array and observes, with a duty cycle close to 15%, the fluorescence light emitted by nitrogen molecules excited by the particles from the air showers. The data set used in the present work consists of 231 cosmic-ray events with zenith angle θ ≤ 80 • and energy E CR ≥ 52 EeV recorded by the SD from January 1, 2004 up to March 31, 2014 1 . These events satisfy a fiducial cut requiring their impact position to be well-contained in the SD array. This ensures an accurate reconstruction of the shower geometry and energy as well as a robust estimation of the exposure.
The arrival directions of these cosmic rays are determined from the relative arrival times of the shower front in the triggered stations. The angular resolution, defined as the radius around the true cosmic-ray direction that would contain 68% of the reconstructed shower directions, is better than 0.9 • for the energies considered here [32]. The energy estimate is obtained from the signals recorded by the SD stations [33,34] and is calibrated using 'hybrid' events (i.e., detected simultaneously by SD and FD) and using the quasi-calorimetric energy determination obtained with the FD [35,36]. The statistical uncertainty in the energy determination is smaller than 12% for the energies used here [34,36]. The systematic uncertainty on the absolute energy scale is 14%.
At energies above 4 EeV the SD trigger is fully efficient, and the determination of the exposure is purely geometrical [37]. For the time period considered and for the applied energy and zenith angle selection, it amounts to 66,452 km 2 sr yr. The adopted zenith-angle selection allows the field of view of the Observatory to extend from −90 • to +45 • in declination.

Telescope Array and its data sample
The Telescope Array (TA) is situated in Utah, USA (39.3 • N, 112.9 • W, 1400 m a.s.l.). It consists of 507 plastic scintillator detectors, each of 3 m 2 in area, located on a 1.2 km square grid and covering an area of approximately 700 km 2 (for further details see [5]). The atmosphere over the surface array is viewed by 38 fluorescence telescopes arranged in 3 stations. The detector has been fully operational since March 2008.
In this analysis, we use surface-detector data recorded between May 11, 2008 and May 4, 2014. For the reconstructed events, the energies determined by the SD array were renormalized by 1/1.27 to match the SD energy scale to that of the fluorescence detector, which was determined calorimetrically [7]. Of these events, 87 met the following criteria: (1) each event triggered at least four SD counters; (2) the zenith angle of the event arrival direction was less than 55 • ; and (3) the reconstructed energy was greater than 57 EeV. The event-selection criteria described above are the same as in Ref. [38] and were optimized to increase the observed number of cosmic rays. The angular resolution of these events is about 1.5 • , while the energy resolution is ∼ 20%. The systematic uncertainty on the energy scale is 22% [7].
To conclude the description of the data, a remark is in order. As already mentioned, the absolute energy scale of UHECRs may contain a systematic error which may be, in principle, different for the two experiments. In accordance with the findings of the TA-Auger Energy Spectrum working group [39], the UHECR spectra measured by Auger and TA may be made coincident in the region around 10 19 eV (the 'ankle' region) by up-shifting the Auger energies (equivalently, down-shifting the TA energies) by ∼ 13%. Assuming that the relative shift is energy-independent, and to make the energy scales compatible, we choose to down-shift the TA energies by 13% in the likelihood analyses presented here. Note that assuming lower energies is a conservative choice with respect to the assumed magnetic deflections for the analyses since the deflections grow as the energy decreases. In view of its relatively small magnitude, the overall effect of this shift on the analyses is minor.

Magnetic deflections of UHECRs
The deflections in the magnetic fields to which, unlike neutrinos, cosmic rays are subject, are crucial in determining the sensitivity of the analyses presented in this paper and any interpretation of the results requires their understanding. The cosmic magnetic fields that deflect the UHECR trajectories are naturally separated into the Galactic and the extragalactic fields. Of these, the extragalactic fields are the least known. The Faraday rotation measurements of extragalactic sources indicate that the extragalactic fields are smaller than ∼ 10 −9 G [40][41][42]. In this scenario, for a correlation length smaller than 1 Mpc, the deflections of protons of energy 10 20 eV over a distance of 50 Mpc are smaller than 2 • .
The Galactic magnetic field can be further divided into regular and turbulent components. The regular part of the Galactic magnetic field is expected to give a dominant contribution to the UHECR deflections. There are several models of the regular Galactic field in the literature (see, e.g., [12] for a recent review). While the averages of the magnitudes of the deflections across the sky predicted in these models are similar, individual deflections in a given direction of the sky vary substantially, so that reliable predictions have not yet been achieved. To get an idea of the expected magnitude of the deflections, a set of 10 5 random cosmic-ray events was simulated following the Auger and TA exposures, weighted in a proportion corresponding to the actual number of events in the data sets of each experiment that we consider in this work. For each of these events we computed the deflection angle according to two recent Galactic magnetic field models [43,44], assuming the primary CR to be a proton and its energy to be 10 20 eV. The corresponding distributions are shown in . Distribution of UHECR deflections in two Galactic magnetic field models marked PT2011 [43] and JF2012 [44] for the regular component. The energies of actual UHECRs are renormalized to show the distributions for E/Z = 100 EeV. The double-peak structure is mostly due to the fact that UHECRs from different Galactic hemispheres undergo different deflections. models were found to be of the same order, ∼ 3 • in both cases and having a rather wide distribution.
Additionally, the cosmic rays are also deflected in the turbulent Galactic magnetic field. There are several estimates of these deflections in the literature [45]. In all cases the contribution due to the turbulent fields was found to be sub-dominant as compared to that due to the regular field.
Apart from the magnetic field, the magnitude of the deflections also depends on the composition of UHECRs. The current measurements of the UHECR composition at the highest energies are still uncertain [13,14], so that the assumed CR deflections are uncertain too.
In summary, it appears unlikely that typical UHECR deflections at E = 10 20 eV are smaller than 2 − 3 • , but may be significantly larger in the case of a heavier composition. In particular, if the deflections are not too large, they scale approximately as Z/E, where Z is the atomic number. In view of these large uncertainties, when considering likelihood tests we perform the analyses for benchmark deflection values, parameterized by the deflection D of a proton at E = 10 20 eV. For simplicity, we model individual deflections as a random variable with a 2-dimensional Gaussian distribution with the energy-dependent width We consider test values of D = 3 • , 6 • and when possible, also of 9 • , which are compat-ible with the distributions for protons displayed in Figure 3 and may also account for somewhat heavier cosmic rays.

UHECR correlation analyses with high-energy cascades and highenergy tracks
Two different analyses, a cross-correlation test and a stacking likelihood analysis, are performed with the neutrino data sample of 39 high-energy cascade events (described in Section 2.1.1) and the 16 high-energy tracks (7 from the sample described in Section 2.1.1 and the 9 described in Section 2.1.2). The total number of UHECRs is N CR = 318 (N Auger = 231 and N TA = 87). For the two analyses, cascades and tracks are considered separately since, due to their different angular resolutions, the angular distance at which a signal (if any) can be observed would be different.

Cross-correlation analysis with high-energy cascades and high-energy tracks
If cosmic rays and neutrinos come from the same sources, then the angular separation between their arrival directions would be related to the magnetic deflections suffered by the cosmic rays, convoluted with the resolution in the determination of their respective arrival directions. This angular distance is unknown a priori due to the poor knowledge of the intervening magnetic fields and the uncertain value of the charges of cosmic rays. The cross-correlation method consists of computing the number of UHECR-neutrino pairs as a function of their angular separation α, n p (α), and comparing it to the expectation from an isotropic distribution of UHECR arrival directions. The angular scan performed in this case is between 1 • and 30 • with a step of 1 • . Due to this scan, the cross-correlation method does not rely on any assumption about the exact value of the strength of the magnetic deflections, unlike the likelihood method described in the following subsection (though there is a large trial factor).
To estimate the significance of any excess, the pre-trial p-value is first computed by evaluating for each angular distance scanned the fraction of isotropic simulations having more pairs than the data. This is done by keeping the neutrino positions fixed and simulating an isotropic distribution of the arrival directions of UHECRs arriving at the Earth according to the corresponding geometric exposures of the Observatories. Then, we find the angular scale for which this fraction is minimized, the minimum value of this fraction being the pre-trial p-value. The post-trial p-value is calculated as the fraction of simulations of isotropic arrival directions of cosmic rays that under a similar analysis would give rise to a smaller pre-trial p-value than the one observed in the data.
Before applying the cross-correlation method to the data, its sensitivity and discovery potential are computed. The sensitivity is defined as the number of signal events, n s , for which one would obtain a p-value of less than 50%, with respect to the isotropic expectation, in 90% of the simulations. The discovery potential at, say, 3σ, is the number of signal events for which one would obtain a p-value corresponding to more than 3σ with 50% probability in repeated experiments. To obtain these quantities we perform simulations in which a number n s of cosmic-ray events are assumed to come from sources. They are generated by smearing the positions of the neutrinos with a Gaussian distribution with a standard deviation obtained from the median angular resolution of each neutrino event, while the remaining N CR − n s events are assumed to be isotropic. The sources are assumed to have equal apparent luminosity at Earth. To generate the n s events, a decision is first made as to whether each event is a cosmic ray measured by Auger or by TA, with a probability given by the fraction of the number of events measured by each Observatory. A source position is then chosen with a probability given by the relative exposure of the corresponding Observatory (shown in Figure 4).
The arrival direction of the i th UHECR event is obtained by smearing the source position with a Gaussian with standard deviation σ(E i ) calculated as: where σ Auger = 0.9 • and σ TA = 1.5 • are the angular resolutions of Auger and TA, respectively, and σ MD (E i ) is the magnetic deflection defined in Eq. (3.1). The UHECR energy E CR is sampled according to the spectrum determined by Auger or TA in the energy range of the UHECRs considered in this analysis, i.e., scaling as E −4.2 [46] and E −4.5 [47], respectively. The values obtained for the sensitivity (90% CL) and for the 3σ discovery potential are listed in Table 1 in terms of n s , the total number of cosmic-ray signal events. They imply for instance that for D = 3 • a total number of around 15.9 cosmic-ray events should originate on average from the directions of the neutrino tracks (or 93.9 from the directions of the neutrino cascades) in order to reach a p-value corresponding to 3σ in 50% of the cases. Note that since in this case we are considering a total of 318 CR directions, this would correspond to a situation in which ∼ 5% (∼ 30%) of the UHECRs are required to have common sources with the neutrino tracks (cascades) in order to reach a 3σ result if D = 3 • (somewhat larger fractions are required for the D = 6 • case). Taking into account that there are 16 tracks and 39 cascades, we see that the average number of required events per source is smaller for the tracks than for the cascades. This is more evident for D = 3 • since in this case the signal events happen to lie quite close to the neutrino directions. Table 1. Sensitivities and post-trial 3σ discovery potentials DP(3σ) for the cross-correlation method in terms of n s , the total number of cosmic-ray signal events.

Stacking analysis with high-energy cascades and high-energy tracks
Stacking a set of sources is a well-known way of accumulating multiple weaker signals to enhance the discovery potential. Even if the background is also stacked, a better sensitivity can be achieved with this procedure than when looking at the same number of sources separately [31]. Since neutrinos are neither deflected nor absorbed on their way to Earth, here we stack the arrival directions of the neutrinos and search for coincident sources of CRs. The unbinned likelihood method is used in this analysis, with the log of the likelihood function defined as: where n s , the number of signal events, is the only free parameter, N CR = N Auger +N TA is the total number of UHECRs, S i Auger and S i TA are the signal probability density functions (PDFs) for Auger and for TA, respectively, and B i Auger and B i TA are the corresponding background PDFs. The Auger signal PDF has the following form: where r i is the angular position of the i th UHECR event, R Auger (δ i ) takes into account the Auger detector response, e.g., the relative exposure for given event declination δ i (see Figure 4) and N src is the number of stacked sources, 39 for the cascades and 16 for the tracks. The last term, S j ( r i , σ(E i )) is the value of the normalized directional likelihood map for the j th source (i.e., the j th neutrino) taken at r i and smeared with a Gaussian with standard deviation σ(E i ), defined in Eq. 4.1, with σ Auger = 0.9 • . The signal PDF for Telescope Array S i TA has the same form as Eq. 4.3, but the relevant parts are replaced with the Telescope Array equivalents, namely R Auger (δ i ) is replaced by the Telescope Array relative exposure R TA (δ i ) (see Figure 4) and the angular resolution is σ TA = 1.5 • .
The background PDFs, B i Auger and B i TA , represent the probabilities of observing a cosmic ray from a given direction assuming an isotropic flux. Therefore they are taken to be the Auger and TA normalized exposures (see Figure 4).
In Figure 5 we show the normalized neutrino likelihood maps convoluted with the exposures of the two CR Observatories, to demonstrate the spread of the different observed neutrino arrival directions visible from each site. The test statistic T S is defined as: T S = 2 log L(n s ) L(n s = 0) .

(4.4)
In accordance with Wilk's theorem [48], the test statistic is expected to follow a distribution close to χ 2 of one degree of freedom, as n s is the only free parameter. This assumption has been explicitly verified and is used to calculate the pre-trial p-values.
To evaluate the sensitivity and discovery potential of this method we perform the same simulation of cosmic-ray signal events as in Section 4.1. The only difference is that in this case the sources are sampled from each neutrino likelihood map. In Table 2 the discovery potential and sensitivity calculations are shown in terms of the total number of cosmic-ray signal events. We note that in these simulations we adopted the same value of deflection parameter D in the log L as the one used to simulate the CR deflections. If the actual deflections were different than the values of D considered in log L the discovery potential would become worse. The comparison of the discovery potentials obtained for this method with those obtained for the cross-correlation in the previous subsection is non-trivial, because in the cross-correlation case they are already penalized for the trial factor of the scan in the angular scale. On the other hand, for the likelihood method a fixed deflection parameter D is adopted, and hence smaller required values of n s may result with the likelihood method. Table 2. Sensitivities and pre-trial 3σ and 5σ discovery potentials (DP(3σ) and DP(5σ)) for the likelihood stacking method in terms of the total number of cosmic-ray signal events.

Results
The results of the cross-correlation method applied to the data are shown in Figures 6a  and 6b for the high-energy tracks and high-energy cascades, respectively. Here the fractional excess relative to the expectation for an isotropic CR distribution, [n p (α)/ n iso p (α) ]− 1, is plotted as a function of the angular separation between the neutrino and cosmicray pairs. The dots represent the data, and the ranges corresponding to the 1σ, 2σ and 3σ confidence bands obtained in isotropic realizations of CR arrival directions are also shown.
For the cross-correlation with the sample of neutrino tracks, the maximum departure from the expectation of an isotropic CR flux occurs at an angular distance of 1 • (Figure 6a), where 0.38 pairs were expected on average and 2 pairs are detected. The post-trial p-value is 28%. For the analysis done using the high-energy cascades (Figure 6b), the smallest pre-trial p-value is obtained at an angular distance of 22 • , for which 575 pairs are observed while 490.3 were expected on average. The post-trial p-value is 5.0 × 10 −4 assuming an isotropic flux of CRs arriving at the Earth.
By looking at Figure 7 one can infer that most of the excess of pairs in this case is due to the fact that there are several high-energy cascades in regions with large densities of UHECRs, i.e., near the Super-Galactic plane and at the TA 'hot spot' [38], which is a 20 • radius region centered at the Galactic coordinates (ℓ, b) ≃ (177 • , 50 • ) 2 .
We also perform some a posteriori tests of the cross-correlation results, considering separately the data of Auger and TA. It was observed that both samples lead to a minimum at 22 • , with the TA post-trial p-value being 9.3 × 10 −4 and that of Auger being 4.1 × 10 −2 . Thus, when considering the entire UHECR data set, this minimum gets reinforced.
The results of the likelihood stacking method are summarized in Table 3. The most significant deviation from the isotropic flux is found for the magnetic deflection parameter D = 6 • for the cascade sample. The observed pre-trial p-value is 2.7 × 10 −4 . Due to this rather small value the post-trial p-value calculation based on generating background-only samples and counting the fraction of those more significant than the result is not feasible. We then conservatively apply a trial factor of 3 to account for the 3 values of the magnetic deflection parameter D used in the analysis 3 . The obtained post-trial p-value is 8.0 × 10 −4 . It is important to stress that all the p-values quoted for both analyses above are with respect to the null hypothesis of an isotropic UHECR flux, as analyses of the distributions of their arrival directions yielded no evidence of anisotropy at discovery level. However, directions with higher densities of UHECRs, such as the TA 'hot spot' [38] and the direction of Cen A [9], have been reported. Hence, as an additional a posteriori study for both analyses, we have also evaluated the significance under the hypothesis of an isotropic distribution of neutrinos. In this case, the UHECR positions have been kept  Table 3. Results for the likelihood stacking analyses with the high-energy tracks and high-energy cascades.
fixed, thus preserving the degree of anisotropy in the arrival directions of CRs. The arrival directions of neutrinos have been simulated producing random right ascensions, while keeping their declination fixed and thus accounting for the declination dependence in the IceCube acceptance. This random-neutrino test assesses the extent to which the p-values previously obtained for the clustering of cosmic rays around the neutrino directions could be affected by chance alignments of the neutrinos with known clustering in the cosmic-ray sky.
For the cross-correlation analysis the post-trial p-value obtained under the hypothesis of an isotropic distribution of neutrinos is 8.5 × 10 −3 . A similar post-trial p-value can be calculated for the likelihood stacking analysis by applying the analysis for all three angular deflection hypothesis D = 3 • , 6 • , and 9 • and selecting the most significant result. This can then be compared to the significance obtained for real data and it was found that in four cases out of 3000 the significance of the generated samples was higher. Thus, for the likelihood stacking analysis the post-trial p-value with respect to the isotropic neutrino flux hypothesis is 1.3 × 10 −3 (i.e., of about 3σ). We see that for both the cross-correlation and the likelihood stacking analyses, the p-values obtained under the null hypothesis of isotropic neutrinos turn out to be larger than the ones obtained under the null hypothesis of isotropic CRs, the differences reflecting the extent to which the original p-values, from the isotropic cosmic-ray hypothesis, are due to an accidental alignment of the neutrinos with the known clustering of the cosmic rays.
An a posteriori scan over the values of the assumed deflection D for the stacking analysis with cascades, shown in Figure 8, reveals that the minimum p-value happens very close to one of the values adopted for the analysis, i.e., D = 6 • .
The angular distance at which an excess would occur in the case of the crosscorrelation includes not only the magnetic deflections at the corresponding CR energies but also the experimental angular uncertainties. In the case of cascades, the angular uncertainty is ∼ 15 • and it is ∼ 1 • for the CRs. Since most CRs in the data set have E CR ∼ 60 EeV, the assumed magnetic deflection, where the smallest p-value is found in the case of the likelihood analysis with the cascades (σ MD (E CR ) = 6 • × 100 EeV/E CR ), is ∼ 10 • in most cases. To translate this into an angular scale where one would find an excess in the cross-correlation analysis (if there were a signal), we add the different angular scales in quadrature, obtaining This scale is comparable to the 22 • where the smallest p-value is found for the cross-correlation performed with the cascades. Hence, the magnetic deflection of the CRs one would infer from the cross-correlation analysis with the cascades is comparable to the value leading to the smallest p-value in the likelihood analysis, even if none of the results are at a level where any strong claims can be made.
The consistent and significant correlation of UHECRs with the high-energy cascades observed in both the cross-correlation and stacking analyses forms a potentially interesting result, which we will continue to monitor in the future.

Method and discovery potential
This analysis looks for excesses in the neutrino point-source data set, consisting of through-going tracks and described in Section 2.1.3, around the directions of the highestenergy UHECRs. We use a likelihood method which is a generalization of the one described in Section 4.2 [24,49], but now the measured positions of the UHECRs are stacked as sources.
The log of the likelihood function is defined as: where n s is the number of signal events in the sample and γ is the spectral index of the neutrino source candidates, assumed to collectively follow an unbroken power-law spectrum ∝ E −γ . N ν is the total number of astrophysical neutrino candidate events in the sample. S tot i is the signal PDF for the stacked sources and B i is the background PDF.
Signal events from an astrophysical source j are expected to be clustered around the direction of the source r j according to the Gaussian distribution: where σ ij = σ 2 i + σ 2 j + σ 2 MD (E j ) accounts for the spread due to the angular resolution of the IceCube event i, that of the CR event j, as well as for the spread due to the assumed magnetic deflection of the CR with energy E j .
The signal PDF S j i also accounts for the energy-(and declination-) dependent response of the IceCube detector, P (δ i , E i |γ), which can be obtained from Monte Carlo simulations of events produced by a source with a power-law spectrum of index γ. The single-source signal PDF is then The total PDF obtained after stacking the N CR UHECR directions is given by where R IC (δ j , γ) is the detector acceptance of IceCube at declination δ j for a source of spectral index γ. Events which are atmospheric muons and/or neutrinos are expected to be distributed uniformly within each declination band. The distribution of the values of the energy estimator E of these events can be indicated by P (δ i , E|φ atm ). The background probability distribution can then be expressed as (5.5) where B(δ i ) is the declination dependence of the sample. The number of signal events n s and the spectral index of the sources γ are unknown, but the best estimates of these two parameters,n s andγ, can be obtained as those maximizing the likelihood. The test statistic is then obtained as The significance of an observation can then be estimated by repeating the process on datasets randomized in right ascension and computing the fraction of randomized samples that produce a value of the T S bigger than that observed in the data. The performance of this kind of search can be quantified by the median source flux required for 5σ discovery, defined similarly to that in Section 4.2 but in terms of the astrophysical neutrino flux. This is calculated by repeating the test on datasets with n s additional simulated signal events injected in the direction of the sources.
Since using all of the 318 events described in Section 2 as sources with their associated angular extensions would essentially cover the whole sky, the possible gain associated to the stacking of more sources is reduced. Restricting the UHECR sample to those with energies above a given threshold E th could give the advantage of having smaller associated magnetic spreads, potentially enhancing the discovery potential. We look for an optimum value of E th for which the per-source flux required for discovery can be minimized as a consequence of these two competing effects. igure 9. Pre-trial discovery potential at the 5σ level. The normalization of the flux per source required is shown as a function of the UHECR threshold energy E th . The results from Northern and Southern skies, as well as the whole sky, are shown for D = 3 • and 6 • .

10
A simulation was then performed to determine the discovery potential of the UHECR stacking analysis obtained for different values of E th (see Figure 9) for both the Southern and Northern hemispheres separately as well as for the whole sky sample. To calculate the discovery potentials, events were injected from point sources in the direction of the UHECRs, adopting source spectra with γ = 2 and considering two deflection hypotheses, D = 3 • and 6 • . Each point source was positioned with a deflection ∆ j from the UHECR, where ∆ j is sampled randomly from a Gaussian distribution of width σ 2 j + σ 2 MD (E j ). As seen in Figure 9, the discovery potentials are characterized by a broad flat region, based on which a value of E th = 85 EeV was chosen. This keeps the 12 highest-energy UHECRs in the Southern sky and 15 in the Northern sky.
Due to the different energy ranges between the neutrino candidate events in the Southern hemisphere (∼100 TeV -100 PeV) and in the Northern hemisphere (∼1 TeV -1 PeV), the flux normalization required for Northern sources turns out to be smaller than for Southern ones. However, thanks to the gain provided by having more stacked sources, the whole sky sensitivity turns out to be similar to the one in the Northern sky alone (Figure 9). We have thus considered just the whole sky sample with E th = 85 EeV in this analysis.

Results
Applying this analysis to the data, all observations were found to be compatible with the background hypothesis. For the deflection hypothesis of D = 3 • , the p-value was found to be 17.3% with a fitted excess of n s ∼ 123 events and a fitted source spectrum of E −3.24 , while the p-value was found to be > 50% for the D = 6 • hypothesis. Accounting for the trial factor due to having tested both the D = 3 • and the D = 6 • hypotheses, the post-trial p-value was found to be 25.4%.

Conclusions
Three analyses have been performed to investigate correlations between the 318 UHECRs measured by the Auger Observatory and the Telescope Array with various samples of IceCube neutrino events, and the results obtained are all below 3.3σ.
The smallest post-trial p-values were obtained when considering IceCube highenergy cascades. In the case of the cross-correlation analysis, the minimum was found for an angular scale of 22 • , with a post-trial p-value of 5.0×10 −4 under the assumption of an isotropic flux of UHECRs and 8.5 × 10 −3 in the a posteriori study under the assumption of an isotropic flux of neutrinos. For the stacking done with the high-energy cascades and a deflection hypothesis of D = 6 • , a post-trial p-value of 8.0 × 10 −4 was obtained under the assumption of an isotropic flux of UHECRs and 1.3 × 10 −3 under the a posteriori assumption of an isotropic flux of neutrinos. The two results correspond to comparable magnetic deflections of CRs. These correlation excesses arise mostly from pairs of events in the region of the TA hot spot and also in regions close to the Super-Galactic plane where some excess of events is seen in the Auger sample, but at present the results are not significant enough to draw any strong conclusions and are also compatible with being fluctuations of an isotropic distribution. It will be interesting to further study these findings with increased statistics in the future. On the other hand, the p-values obtained in the studies involving neutrino tracks are large, being totally compatible with the absence of correlations between the arrival directions of the neutrino tracks and the UHECRs. Similarly, the p-values obtained from the stacking analyses with the point-source sample and the 27 highest-energy UHECRs show no significant excess of neutrinos in the directions of these UHECRs.
We note that the absence of strong correlations between neutrinos and UHECRs may have different causes. In particular, one has to keep in mind that if the strong suppression in the cosmic-ray flux observed above 40-50 EeV is due to the attenuation of protons in the CMB through photopion production (or of nuclei due to photodisintegrations), the majority of the observed UHECRs above 50 EeV should have been produced in sources relatively nearby, i.e., within ∼ 200 Mpc. On the other hand, the astrophysical neutrinos may come from sources at any distance since neutrinos are not limited by absorption. Hence, only a small fraction of the extragalactic neutrinos, those from sources relatively nearby, may have the same sources as the observed UHECRs. Considering for instance that the neutrino sources extend to at least z=1, the overall extragalactic contribution arising from sources closer than 200 Mpc is not expected to exceed a few percent (in the case of no source evolution and uniform source distribution) and could be particularly suppressed if the evolution of the sources is strong. Another difficulty in observing common UHECR and neutrino sources is that while the neutrinos arrive straight from their sources, the deflections of charged cosmic rays in Galactic and extragalactic magnetic fields induce sizeable delays between the arrival times of CRs and those of the neutrinos, which are typically much larger than the duration of the experiments. Hence, if the sources are burst-like rather than steady, the burst of neutrinos would arrive long before the associated UHECRs and their simultaneous observation could not be feasible. Finally, the sources capable of producing PeV neutrinos, which are those able to accelerate protons up to energies of a few tens of PeV, may only rarely be able to also accelerate cosmic rays to ultrahigh-energies. Hence, it may be natural to expect that only a small fraction of the observed energetic neutrinos come from sources that also produce UHECRs.
In conclusion, further insight will potentially arise from increased statistics and eventually with the inclusion of CR composition information that may become available so as to better model possible effects of magnetic field deflections. This may help to understand if there is a contribution in the astrophysical neutrino signal observed by IceCube correlated to the sources of the observed UHECR.
The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment and effort from the technical and administrative staff in Malargüe. We are very grateful to the following agencies and organizations for financial support: Air Force. We also wish to thank the people and the officials of Millard County, Utah for their steadfast and warm support. We gratefully acknowledge the contributions from the technical staffs of our home institutions. An allocation of computer time from the Center for High Performance Computing at the University of Utah is gratefully acknowledged. Table 4. List of the neutrino cascade-like events used in the analyses. The ID number corresponds to the ID in [21,22] 170.5 6.0 11.6 Table 5. List of the neutrino track-like events used in the analyses. The first seven ID numbers correspond to the IDs in [21,22]. The rest of the ID numbers correspond to the high-energy tracks from [23] described in Section 2.1.2. For these events, the 'most probable muon energy' is an estimate of the muon energy assuming the best fit spectrum resulting from the flux fit in that analysis [23], which is different from the deposited energy specified for HESE tracks and cascades.