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

The IceCube Neutrino Observatory (Aartsen et al. 2017d) is an ice-Cherenkov detector sensitive to neutrinos with energies ≥5 GeV. It is located at the geographic South Pole, about 1.45−2.45 km deep in the ice. Its main component consists of a volume of about 1 km³ glacial ice instrumented with 5160 photomultipliers that are connected to the surface by 86 cable strings.

Two classes of neutrino-induced events can be phenomenologically distinguished: elongated, track-like events that 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 owing to their different topologies, but they 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. 2014c; 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 and elaborate selection procedures are necessary to filter out atmospheric muons (Aartsen et al. 2017b, 2017c, 2018c). 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 the case of the very down-going region, by vetoing events where an atmospheric shower is observed in IceTop, IceCube'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 through-going 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 owing 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.

Zoom In Zoom Out Reset image size

The PS data set consists of a combination of data collected from 7 yr of operation between 2008 and 2015 that were used for point-source searches (Aartsen et al. 2017b) and data from 3.5 yr of operation between 2015 and 2018 that were selected for the real-time gamma-ray follow-up (GFU) program of IceCube (Aartsen et al. 2017c, 2018c). 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 Section 2.2). Note that analysis 1 is based not on the median of the angular resolution but on the estimator for the angular resolution on an event-by-event basis (σ).

Zoom In Zoom Out Reset image size

The HESE data set contains 76 cascade-like events and 26 track-like events that have been collected between 2010 and 2017, as presented in Wandkowsky (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 decl. ≥ −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 nonatmospheric 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¹⁸ GeV⁻¹ cm⁻² s⁻¹ sr⁻¹ (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. 2017c, 2018a). 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−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 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³. Each string is equipped with 25 storeys of three 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 Section 2.1). Particles above the Cherenkov threshold induce coherent radiation emitted in a cone with a characteristic angle of 42° in water. 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 ⁴⁰K decays and bioluminescence (Aguilar et al. 2007). ANTARES is thus also capable of detecting charged−current and neutral−current neutrino interactions of all flavors.

For analysis 1, we use all track-like events from the 11 yr 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 track-like 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 04 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 yr 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 Section 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 contribute substantially to the signal acceptance for soft source spectra.

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². 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 Aab et al. (2020b). 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; Aab et al. 2020b), and the angular uncertainty is less than 09 (Bonifazi & The Pierre Auger Collaboration 2009).

The UHECR data set used here consists of 324 events observed with the SD between 2004 January and 2017 April (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 decl. All energies are scaled up by a constant factor of 14% to match both Auger and TA energies on a common energy scale (see Section 2.4), following the recommendation of the Auger–TA joint working group (Biteau et al. 2019). This scaling factor was detemined 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 a function of decl. can be found in Figure 3, which amounts to an integrated geometric exposure of 91,300 km² yr sr.

Zoom In Zoom Out Reset image size

The 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² each. The SD stations are located on a square grid with 1.2 km separation, which extends over an area of 700 km² (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 decl. In addition to the SD, there are three fluorescence telescope stations, instrumented with 12–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 yr 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 2008 May and 2017 May with reconstructed energies ≥57 EeV and zenith angles ≤55°. At these energies, the angular uncertainty is about 15. The statistical uncertainty on the reconstructed energy is about 15%−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 Section 2.3), following the recommendation of the Auger–TA joint working group (Biteau et al. 2019).

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 decl. 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 decl., with an integrated total exposure of 11,600 km² yr sr (Biteau et al. 2019).

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. 2017b, 2019, 2020c). 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 Sections 2.1 and 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 $\vec{{\rm{\Omega }}}=(\alpha ,\delta )$ in R.A. and decl. consists of the sum of a signal probability density function (pdf), S, and a background pdf, B,

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L }\left({n}_{s},\gamma \right) & = & \displaystyle \prod _{j=1}^{7}\displaystyle \prod _{i=1}^{{N}_{j}}{f}_{j}(\sin \delta ,\gamma )\displaystyle \frac{{n}_{s}}{{N}_{j}}\cdot S({\vec{{\rm{\Omega }}}}_{i},{E}_{i},{\sigma }_{i}| \gamma ,\vec{{\rm{\Omega }}})\\ & & +\left(1-{f}_{j}(\sin \delta ,\gamma )\displaystyle \frac{{n}_{s}}{{N}_{j}}\right)\cdot B({\delta }_{i},{E}_{i}).\end{array}\end{eqnarray} \tag{ 2 }$$

The likelihood product with index i runs over all neutrino events within each data set with index j (see Sections 2.1 and 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 the same way as described in Aartsen et al. (2017b). 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 R.A. and decl., summarized as ${\vec{{\rm{\Omega }}}}_{i}=({\alpha }_{i},{\delta }_{i})$, reconstructed energy, E_i , and angular error estimator, σ_i . The signal pdf consists of two terms: the first term is a decl.-dependent reconstructed energy distribution, where the underlying neutrino flux is modeled as a power law,

$$\begin{eqnarray}&&\frac{{d}^{}{{\rm{\Phi }}}^{\nu }}{d{E}^{}}=\frac{{{dN}}^{\nu }}{{dE}\,{dA}\,{dt}\,d{\rm{\Omega }}}={{\rm{\Phi }}}_{0}^{\nu }\cdot {\left(\frac{E}{1\,\mathrm{GeV}}\right)}^{-\gamma }.\end{eqnarray} \tag{ 3 }$$

Here ${{\rm{\Phi }}}_{0}^{\nu }$ 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 a function of reconstructed energy, E_i , and decl., δ_i , from randomized experimental data. A full description of the signal and background pdfs, S and B, can be found in Aartsen et al. (2017b). The proportionality factor, ${f}_{j}(\sin \delta ,\gamma )$, is the relative signal acceptance per neutrino data set calculated from the expected number of signal events via

$$\begin{eqnarray}&&\begin{array}{l}{f}_{j}(\sin \delta ,\gamma )=\displaystyle \frac{{n}_{j}}{{n}_{\mathrm{tot}}}\\ \,\,=\,\displaystyle \frac{{T}^{j}\displaystyle \int (d{{\rm{\Phi }}}^{\nu }/{dE}){A}_{\mathrm{eff}}^{j}(E,\sin \delta )\,{dE}}{\sum _{j=1}^{7}{T}^{j}\displaystyle \int (d{{\rm{\Phi }}}^{\nu }/{dE}){A}_{\mathrm{eff}}^{j}(E,\sin \delta )\,{dE}}.\end{array}\end{eqnarray} \tag{ 4 }$$

This factor is used to correctly weight the signal contribution of each data set j for a given source decl. and spectral index. Per data set j, the live time is denoted with T^j and the effective area is denoted with ${A}_{\mathrm{eff}}^{j}$. Figure 5 shows the proportionality factor of each data set as a function of decl. 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).

Zoom In Zoom Out Reset image size

The best-fit signal parameters, ${\hat{n}}_{s}$ and $\hat{\gamma }$, at a given source position, $\vec{{\rm{\Omega }}}$, are obtained with the maximum likelihood method. The number of events and the proportionality factor are related to the neutrino flux using the respective live time and effective area of each data set via

$$\begin{eqnarray}&&{n}_{s}=\sum _{j=1}^{7}{T}^{j}\int \frac{{d}^{}{{\rm{\Phi }}}^{\nu }}{d{E}^{}}{A}_{\mathrm{eff}}^{j}(E,\sin \delta )\,{dE}.\end{eqnarray} \tag{ 5 }$$

The corresponding significance of a source at position $\vec{{\rm{\Omega }}}$ is evaluated using a likelihood ratio test, which yields the test statistic (TS)

$$\begin{eqnarray}&&{\mathrm{TS}}_{\nu }(\vec{{\rm{\Omega }}})=2{\left.\mathrm{log}\left(\displaystyle \frac{{ \mathcal L }\left({\hat{n}}_{s},\hat{\gamma }\right)}{{ \mathcal L }({n}_{s}=0)}\right)\right|}_{\vec{{\rm{\Omega }}}}.\end{eqnarray} \tag{ 6 }$$

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 02. This results in a sky map of the TS values at each grid center, ${\vec{{\rm{\Omega }}}}_{\mathrm{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 ${\vec{{\rm{\Omega }}}}_{k}^{\mathrm{CR}}$, which possibly originated in the direction of ${\vec{{\rm{\Omega }}}}_{\mathrm{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, ${\sigma }_{\mathrm{CR}}^{2}$, which is 09 for Auger events and 15 for TA events,

$$\begin{eqnarray}&&{{ \mathcal P }}_{k}\propto \exp \left(-\displaystyle \frac{{\left({\vec{{\rm{\Omega }}}}_{\mathrm{grid}}-{\vec{{\rm{\Omega }}}}_{k}^{\mathrm{CR}}\right)}^{2}}{2({\sigma }_{\mathrm{MD}}^{2}+{\sigma }_{\mathrm{CR}}^{2})}\right).\end{eqnarray} \tag{ 7 }$$

The term ${{ \mathcal P }}_{k}$ acts as a spatial constraint, which is multiplied by the neutrino likelihood function defined in Equation (2) via ${ \mathcal L }\to { \mathcal L }\cdot {{ \mathcal P }}_{k}$. Effectively, the constraint is added as a logarithmic term to the neutrino TS defined for each grid point in Equation (6) via $\mathrm{TS}\to \mathrm{TS}+2\cdot \mathrm{log}({{ \mathcal P }}_{k})$. Finally, the maximum of the combined UHECR−neutrino TS is found at the best-fit grid point, ${\hat{\vec{{\rm{\Omega }}}}}_{\mathrm{grid}}$, via

$$\begin{eqnarray}&&{\mathrm{TS}}_{k}^{\mathrm{comb}}({\hat{\vec{{\rm{\Omega }}}}}_{\mathrm{grid}})=\max \left({\mathrm{TS}}_{\nu }({\vec{{\rm{\Omega }}}}_{\mathrm{grid}})-\displaystyle \frac{{\left({\vec{{\rm{\Omega }}}}_{\mathrm{grid}}-{\vec{{\rm{\Omega }}}}_{k}^{\mathrm{CR}}\right)}^{2}}{{\sigma }_{\mathrm{MD}}^{2}+{\sigma }_{\mathrm{CR}}^{2}}\right).\end{eqnarray} \tag{ 8 }$$

The maximum marks the neutrino source candidate that 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 point-source correlation analysis and first outlined in Section 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. 2020a) and with gravitational wave events (Aartsen et al. 2020b). Note that the analysis described here is improved with respect to the previous search (The IceCube Collaboration et al. 2016, Section 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₀ = 3° for all UHECRs, which is then converted into the spatial constraint of an individual UHECR event using Equations (1) and (7). The corresponding sensitivity and 3σ-discovery potential are evaluated based on the normalization of the neutrino flux per source, ${{\rm{\Phi }}}_{0}^{\nu }$ (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 ${{\rm{\Phi }}}_{0}^{\nu }$ in the case of a null measurement, while the discovery potential is defined as the median of ${{\rm{\Phi }}}_{0}^{\nu }$ 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 Equations (7) and (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 R.A. The signal PEs reflect the six different signal hypotheses; one signal PE is based on experimental neutrino data randomized in R.A., 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 owing 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⁻² power-law 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 rounded-up product of the correlation fraction with the numbers of UHECRs, ${N}_{\nu }^{\mathrm{src}}=\unicode{x02308}{f}_{\mathrm{corr}}{N}_{\mathrm{CR}}\unicode{x02309}$, 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 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 ${N}_{\nu }^{\mathrm{src}}\lt {N}_{\mathrm{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 and 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 neutrino sources, which in turn need to have a higher flux normalization to be detected by the analysis.

Zoom In Zoom Out Reset image size

Table 2. Sensitivity, Discovery Potential, and 90% C.L. Upper Limit for the Different Analysis Parameters, Which Are The UHECR Energy Cut, E_cut, Magnetic Deflection, D₀ · C, and Correlation Fraction, f_corr, for the Point-source Correlation Analysis

Magnetic deflection, D0 · C	3°	3°	3°	6°	6°	6°
Energy cut, Ecut	70 EeV	85 EeV	100 EeV	70 EeV	85 EeV	100 EeV
Number of neutrino sources, ${{\boldsymbol{N}}}_{{\boldsymbol{\nu }}}^{\mathrm{src}}$
fcorr = 0.1	22	9	4	22	9	4
fcorr = 0.5	106	41	20	106	41	20
fcorr = 1	211	82	40	211	82	40
Sensitivity
fcorr = 0.1	5.6	7.1	9.8	5.4	8.5	9.5
fcorr = 0.5	2.4	2.6	3.1	2.6	3.4	4.2
fcorr = 1	1.4	1.7	2.1	1.8	2.2	2.9
3σ disc. potential
fcorr = 0.1	7.2	8.5	11.6	6.7	8.0	10.8
fcorr = 0.5	3.7	3.8	4.9	3.7	4.1	5.1
fcorr = 1	2.6	2.6	3.2	2.7	3.1	3.7
90% C.L. upper limit
fcorr = 0.1	6.4	9.2	9.8	6.7	10.9	10.6
fcorr = 0.5	2.8	3.6	3.1	3.3	4.7	4.5
fcorr = 1	1.7	2.3	2.1	2.3	3.2	3.1
Pretrial p -value	0.33	0.23	>0.5	0.19	0.097	0.43

Note. The values of sensitivity, discovery potential, and limit are given as normalization of the neutrino flux per source in units of 10⁻¹⁰ GeV⁻¹ cm⁻² s⁻¹. The neutrino flux per source is modeled with a power law of the form $d{\rm{\Phi }}/{dE}={{\rm{\Phi }}}_{0}\cdot {\left(E/1\,\mathrm{GeV}\right)}^{-2}$. The last row states all six experimentally obtained pretrial p-values with respect to the null hypothesis of an isotropic neutrino flux.

Download table as:  ASCIITypeset image

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 owing to the small number of sources present. ²⁴¹ 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.

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 high-energy 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

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\,{ \mathcal L }({n}_{s}) & = & \displaystyle \sum _{i=1}^{{N}_{\mathrm{Auger}}}\mathrm{log}\left(\displaystyle \frac{{n}_{s}}{{N}_{\mathrm{CR}}}{S}_{\mathrm{Auger}}^{i}+\displaystyle \frac{{N}_{\mathrm{CR}}-{n}_{s}}{{N}_{\mathrm{CR}}}{B}_{\mathrm{Auger}}^{i}\right)\\ & & +\displaystyle \sum _{i=1}^{{N}_{\mathrm{TA}}}\mathrm{log}\left(\displaystyle \frac{{n}_{s}}{{N}_{\mathrm{CR}}}{S}_{\mathrm{TA}}^{i}+\displaystyle \frac{{N}_{\mathrm{CR}}-{n}_{s}}{{N}_{\mathrm{CR}}}{B}_{\mathrm{TA}}^{i}\right),\end{array}\end{eqnarray} \tag{ 9 }$$

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}_{\mathrm{Auger}}^{i}$ and ${S}_{\mathrm{TA}}^{i}$, and the background pdfs, ${B}_{\mathrm{Auger}}^{i}$ and ${B}_{\mathrm{TA}}^{i}$, are stacked for all high-energy neutrinos such that n_s is a global parameter. The background pdfs per UHECR detector, ${B}_{\det }^{i}$, are the normalized expectations for an isotropic UHECR flux as a function of decl., which correspond to the normalized detector exposures (see Figure 3). The signal pdf per UHECR detector is defined as

$$\begin{eqnarray}\begin{array}{rcl}{S}_{\det }^{i}(\vec{{{\rm{\Omega }}}_{i}},{E}_{i}) & = & {R}_{\det }({\delta }_{i})\cdot \displaystyle \sum _{j=1}^{{N}_{\nu }}{S}_{j}(\vec{{{\rm{\Omega }}}_{i}},\sigma ({E}_{i}))\\ \mathrm{with}\,{\sigma }^{2} & = & {\sigma }_{\mathrm{MD}}^{2}+{\sigma }_{\mathrm{CR}}^{2},\end{array}\end{eqnarray} \tag{ 10 }$$

as a function of the UHECR arrival direction, $\vec{{{\rm{\Omega }}}_{i}}$, and energy, E_i . The term R_det is the relative exposure of each detector as a function of the decl. 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}(\vec{{{\rm{\Omega }}}_{i}},\sigma ({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, $\vec{{{\rm{\Omega }}}_{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 Galactic northern hemisphere or southern hemisphere (see Figure 4), the baseline magnetic deflection is thus

$$\begin{eqnarray}\begin{array}{rcl}D & = & C\cdot \left({D}_{\mathrm{North}},{D}_{\mathrm{South}}\right)\\ & = & C\cdot \left(2\mathop{.}\limits^{^\circ }4,3\mathop{.}\limits^{^\circ }7\right)\ \mathrm{with}\ C\in \{1,2,3\}.\end{array}\end{eqnarray} \tag{ 11 }$$

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, ${\hat{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,

$$\begin{eqnarray}&&\mathrm{TS}=2\mathrm{log}\displaystyle \frac{{ \mathcal L }({\hat{n}}_{s})}{{ \mathcal L }({n}_{s}=0)}.\end{eqnarray} \tag{ 12 }$$

The significance is then determined based on the assumption that the background expectation is distributed according to a χ² function with 1 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 hemispheres 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, as in The IceCube Collaboration et al. (2016). For the signal PEs, a number n_s of UHECR arrival directions are distributed randomly based on their respective spatial signal terms in Equation (10), ${S}_{j}(\vec{{{\rm{\Omega }}}_{i}},\sigma ({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 sensitivity and 3σ-discovery potential using neutrino tracks require much fewer UHECRs than when using neutrino cascades, which is expected owing to the large differences in angular reconstruction uncertainty of the two event topologies.

Table 3. Sensitivity, 3σ Discovery Potential, and 90% C.L. Upper Limits in Terms of Number of UHECRs (n_s ; see Equation (9)), as well as the Pretrial p-values for the UHECR Stacking Analysis with the Samples of High-energy Tracks and Cascades Assuming an Isotropic Flux of UHECRs

Analysis parameters
DN/S · C	(24, 37)	(48, 74)	(72, 111)
Sensitivity
tracks	9.8	9.8	10.3
cascades	57	61	64
3 σ disc. potential
tracks	20	21	21
cascades	113	129	135
90% C.L. upper limit
tracks	9.8	9.8	10.3
cascades	55	76	101
Pretrial p -values
tracks	>0.5	>0.5	>0.5
cascades	>0.5	0.38	0.26

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{\mathrm{obs}}(\alpha )}{\langle {n}_{\mathrm{bckg}}(\alpha )\rangle }-1.\end{eqnarray} \tag{ 13 }$$

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 Section 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 R.A. The different types of background PEs approximate an isotropic flux of UHECRs 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).

The test statistic of experimental data is obtained with the neutrino point-source correlation analysisas described in Section 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 overestimate the limits based on an underfluctuation 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%–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.

The UHECR stacking analysis is performed for track-like 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 Galactic northern and southern hemispheres of (72, 111). 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 overestimate the limits based on an underfluctuation in data. All p-values and the corresponding limits are reported in Table 3, together with the sensitivity and discovery potential.

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 separations, 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 p-values 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.

Zoom In Zoom Out Reset image size