<figure><img src="image_10.png" /><figcaption>Figure 6. Sensitivity (left) and 3σ-discovery potential (right) of the neutrino flux per source as a function of the UHECR energy cut in EeV, for the three correlation fractions and for the two scaling factors we used.</figcaption></figure>

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.

## 4.2. Unbinned Likelihood-stacking Analysis of UHECRs and High-energy Neutrinos

The second analysis is based on using the neutrinos with a high probability of astrophysical origin as markers for the location of the sources of UHECRs and neutrinos. The UHECR events are stacked using an unbinned likelihood analysis with the arrival directions of the 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

$$ \log \mathcal{L}(n_s) = \sum_{i=1}^{N_{\text{Auger}}} \log \left( \frac{n_s}{N_{\text{CR}}} S_{\text{Auger}}^i + \frac{N_{\text{CR}} - n_s}{N_{\text{CR}}} B_{\text{Auger}}^i \right) + \sum_{i=1}^{N_{\text{TA}}} \log \left( \frac{n_s}{N_{\text{CR}}} S_{\text{TA}}^i + \frac{N_{\text{CR}} - n_s}{N_{\text{CR}}} B_{\text{TA}}^i \right), \quad (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_{\text{Auger}}$ and $N_{\text{TA}}$, where the total number of UHECR events is $N_{\text{CR}} = N_{\text{Auger}} + N_{\text{TA}}$. In contrast to the likelihood function of analysis 1, the signal pdfs, $S_{\text{Auger}}^i$ and $S_{\text{TA}}^i$, and the background pdfs, $B_{\text{Auger}}^i$ and $B_{\text{TA}}^i$, are stacked for all

high-energy neutrinos such that $n_s$ is a global parameter. The background pdfs per UHECR detector, $B_{\text{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

$$ S_{\text{det}}^i(\vec{\Omega}_i, E_i) = R_{\text{det}}(\delta_i) \cdot \sum_{j=1}^{N_i} S_j(\vec{\Omega}_i, \sigma(E_i)) \\ \text{with } \sigma^2 = \sigma_{\text{MD}}^2 + \sigma_{\text{CR}}^2, \qquad (10) $$

as a function of the UHECR arrival direction, $\vec{\Omega}_i$, and energy, $E_i$. The term $R_{\text{det}}$ is the relative exposure of each detector as a function of the decl. per UHECR event, $\delta_i$. Figure 3 shows the absolute directional exposures, which are each normalized to 1 over the whole sky to obtain $R_{\text{det}}$. The sum runs over all neutrino events, $N_i$, where the terms of $S_j(\vec{\Omega}_i, \sigma(E_i))$ are the spatial likelihood maps obtained from the neutrino directional reconstruction smeared with the UHECR uncertainty with width $\sigma(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{\Omega}_i$. Similar to analysis 1, three different deflections corresponding to scaling factors of $C \in \{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{aligned} D &= C \cdot (D_{\text{North}}, D_{\text{South}}) \\ &= C \cdot (2.4, 3.7) \text{ with } C \in \{1, 2, 3\}. \end{aligned} \qquad (11) $$

The final value of the deflection, $\sigma_{\text{MD}}$, is then calculated based on the UHECR energy using Equation (1).
<caption>**Table 3**</caption>

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

<table><thead><tr><th colspan="4">Analysis parameters</th></tr><tr><th>D<sub>N/S</sub> · C</th><th>(2°,4, 3°,7)</th><th>(4°,8, 7°,4)</th><th>(7°,2, 11°,1)</th></tr></thead><tbody><tr><td><strong>Sensitivity</strong></td><td></td><td></td><td></td></tr><tr><td>tracks</td><td>9.8</td><td>9.8</td><td>10.3</td></tr><tr><td>cascades</td><td>57</td><td>61</td><td>64</td></tr><tr><td><strong>3σ disc. potential</strong></td><td></td><td></td><td></td></tr><tr><td>tracks</td><td>20</td><td>21</td><td>21</td></tr><tr><td>cascades</td><td>113</td><td>129</td><td>135</td></tr><tr><td><strong>90% C.L. upper limit</strong></td><td></td><td></td><td></td></tr><tr><td>tracks</td><td>9.8</td><td>9.8</td><td>10.3</td></tr><tr><td>cascades</td><td>55</td><td>76</td><td>101</td></tr><tr><td><strong>Pretrial p-values</strong></td><td></td><td></td><td></td></tr><tr><td>tracks</td><td>&gt;0.5</td><td>&gt;0.5</td><td>&gt;0.5</td></tr><tr><td>cascades</td><td>&gt;0.5</td><td>0.38</td><td>0.26</td></tr></tbody></table>

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$,

$$TS = 2 \log \frac{\mathcal{L}(\hat{n}_s)}{\mathcal{L}(n_s = 0)} \quad (12)$$

The significance is then determined based on the assumption that the background expectation is distributed according to a $\chi^2$ 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 resem- bling the signal and background hypothesis. For the back- ground 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{\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.

### 4.3. Two-point Angular Correlation of UHECRs and High-energy Neutrinos

The 2pt-correlation analysis relies on counting pairs of UHECRs and the high-energy neutrinos, where the angular separation between their arrival directions is within a maximum angular separation, $\alpha$. The observed number of pairs within this radius, $n_{\text{obs}}(\alpha)$, is compared to the mean number of pairs expected from pure background, i.e., uncorrelated arrival directions, $\langle n_{\text{bckg}}(\alpha) \rangle$. The relative excess of pairs is defined as

$$\frac{n_{\text{obs}}(\alpha)}{\langle n_{\text{bckg}}(\alpha) \rangle} - 1. \quad (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 rando- mized 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).

## 5. Results

### 5.1. Unbinned Neutrino Point-source Search with UHECR Information

The test statistic of experimental data is obtained with the neutrino point-source correlation analysis 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_{\text{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_{\text{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_{\text{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
<caption>**Table 3**</caption>

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

<table><thead><tr><th colspan="4">Analysis parameters</th></tr><tr><th>D<sub>N/S</sub> · C</th><th>(2°,4, 3°,7)</th><th>(4°,8, 7°,4)</th><th>(7°,2, 11°,1)</th></tr></thead><tbody><tr><td><strong>Sensitivity</strong></td><td></td><td></td><td></td></tr><tr><td>tracks</td><td>9.8</td><td>9.8</td><td>10.3</td></tr><tr><td>cascades</td><td>57</td><td>61</td><td>64</td></tr><tr><td><strong>3σ disc. potential</strong></td><td></td><td></td><td></td></tr><tr><td>tracks</td><td>20</td><td>21</td><td>21</td></tr><tr><td>cascades</td><td>113</td><td>129</td><td>135</td></tr><tr><td><strong>90% C.L. upper limit</strong></td><td></td><td></td><td></td></tr><tr><td>tracks</td><td>9.8</td><td>9.8</td><td>10.3</td></tr><tr><td>cascades</td><td>55</td><td>76</td><td>101</td></tr><tr><td><strong>Pretrial p-values</strong></td><td></td><td></td><td></td></tr><tr><td>tracks</td><td>&gt;0.5</td><td>&gt;0.5</td><td>&gt;0.5</td></tr><tr><td>cascades</td><td>&gt;0.5</td><td>0.38</td><td>0.26</td></tr></tbody></table>

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$,

$$TS = 2 \log \frac{\mathcal{L}(\hat{n}_s)}{\mathcal{L}(n_s = 0)} \quad (12)$$

The significance is then determined based on the assumption that the background expectation is distributed according to a $\chi^2$ 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 resem- bling the signal and background hypothesis. For the back- ground 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{\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.

### 4.3. Two-point Angular Correlation of UHECRs and High-energy Neutrinos

The 2pt-correlation analysis relies on counting pairs of UHECRs and the high-energy neutrinos, where the angular separation between their arrival directions is within a maximum angular separation, $\alpha$. The observed number of pairs within this radius, $n_{\text{obs}}(\alpha)$, is compared to the mean number of pairs expected from pure background, i.e., uncorrelated arrival directions, $\langle n_{\text{bckg}}(\alpha) \rangle$. The relative excess of pairs is defined as

$$\frac{n_{\text{obs}}(\alpha)}{\langle n_{\text{bckg}}(\alpha) \rangle} - 1. \quad (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 rando- mized 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).

## 5. Results

### 5.1. Unbinned Neutrino Point-source Search with UHECR Information

The test statistic of experimental data is obtained with the neutrino point-source correlation analysis 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_{\text{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_{\text{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_{\text{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
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Relative excess of pairs, $n_{obs}/\langle n_{exp} \rangle - 1$, as a function of the maximum angular separation of the neutrino and UHECR pairs. The experimental result is shown as black circles, while the blue color bands show the regions containing the $1\sigma$, $2\sigma$, and $3\sigma$ fluctuations from background PEs based on an isotropic distribution of high-energy neutrinos. The results for the track-like neutrinos are shown on the left and for the cascade-like neutrino on the right.</figcaption></figure>

<table><caption>**Table 4**
Posttrial *p*-value and Best-fit Angular Separation for the 2pt-correlation Analysis Obtained with the Neutrino Data Sets of High-energy Tracks and Cascades, as Stated in the First Column</caption><thead><tr><th>Event Type</th><th>Background Hypothesis</th><th>Separation</th><th>Posttrial p-value</th></tr></thead><tbody><tr><td>Tracks</td><td>isotropic neutrinos</td><td>14°</td><td>0.23</td></tr><tr><td>Cascades</td><td>isotropic neutrinos</td><td>16°</td><td>0.15</td></tr><tr><td>Tracks</td><td>isotropic UHECRs</td><td>10°</td><td>&gt;0.5</td></tr><tr><td>Cascades</td><td>isotropic UHECRs</td><td>16°</td><td>0.18</td></tr></tbody></table>

**Note.** The *p*-values are corrected for choosing the largest deviation out of all maximum angular separations. The respective background hypotheses are stated in the second column.

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$.

## 5.2. Unbinned Likelihood-stacking Analysis of UHECRs and High-energy Neutrinos

The UHECR stacking analysis is performed for track-like and cascade-like high-energy neutrinos separately for three different scaling factors of $C \in \{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 (7°, 2, 11°, 1). Based on the observed TS values, 90% C.L. limits on the number of UHECR events correlated to the high-energy neutrinos are calculated using the PEs of the corresponding signal hypothesis for each of the six tests. Since all *p*-values using the track-like neutrinos are larger than 0.5, the limits are set equal to the sensitivity in order to not 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.

## 5.3. Two-point Angular Correlation of UHECRs and High-energy Neutrinos

The 2pt-correlation analysis is applied to the sets of neutrino tracks and cascades separately, as well as for the background hypothesis of an isotropic UHECRs flux and an isotropic high-energy neutrino flux. The significance of the result with respect to the background hypothesis is corrected for the scan over the separation angles. This results in four *p*-values and four respective best-fit angular 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.

# 6. Conclusions

Three complementary analyses have been performed on the UHECR data sets provided by the Pierre Auger and the Telescope Array Collaborations combined with the high-energy and full-sky neutrino data sets provided by the IceCube and ANTARES Collaborations. For both the UHECR and neutrino data sets, the combination of data from the two respective observatories provides a field of view over the entire sky. None of the analyses found a result incompatible with the assumed background hypotheses of either an isotropic neutrino flux or an isotropic UHECR flux. This is an important update on the results reported in The IceCube Collaboration et al. (2016), where *p*-values close to the $3\sigma$ level were reported when applying analyses 2 and 3 to cascade-like neutrino events.

Based on the results, 90% C.L. upper limits are calculated for analysis 1 on the flux of neutrinos from point-like sources correlated with UHECRs and for analysis 2 on the number of UHECR events correlated with high-energy neutrinos. Analysis 1 reports upper limits on the correlated neutrino flux per
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Relative excess of pairs, $n_{obs}/\langle n_{exp} \rangle - 1$, as a function of the maximum angular separation of the neutrino and UHECR pairs. The experimental result is shown as black circles, while the blue color bands show the regions containing the $1\sigma$, $2\sigma$, and $3\sigma$ fluctuations from background PEs based on an isotropic distribution of high-energy neutrinos. The results for the track-like neutrinos are shown on the left and for the cascade-like neutrino on the right.</figcaption></figure>

<table><caption>**Table 4**
Posttrial *p*-value and Best-fit Angular Separation for the 2pt-correlation Analysis Obtained with the Neutrino Data Sets of High-energy Tracks and Cascades, as Stated in the First Column</caption><thead><tr><th>Event Type</th><th>Background Hypothesis</th><th>Separation</th><th>Posttrial p-value</th></tr></thead><tbody><tr><td>Tracks</td><td>isotropic neutrinos</td><td>14°</td><td>0.23</td></tr><tr><td>Cascades</td><td>isotropic neutrinos</td><td>16°</td><td>0.15</td></tr><tr><td>Tracks</td><td>isotropic UHECRs</td><td>10°</td><td>&gt;0.5</td></tr><tr><td>Cascades</td><td>isotropic UHECRs</td><td>16°</td><td>0.18</td></tr></tbody></table>

**Note.** The *p*-values are corrected for choosing the largest deviation out of all maximum angular separations. The respective background hypotheses are stated in the second column.

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$.

## 5.2. Unbinned Likelihood-stacking Analysis of UHECRs and High-energy Neutrinos

The UHECR stacking analysis is performed for track-like and cascade-like high-energy neutrinos separately for three different scaling factors of $C \in \{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 (7°, 2, 11°, 1). Based on the observed TS values, 90% C.L. limits on the number of UHECR events correlated to the high-energy neutrinos are calculated using the PEs of the corresponding signal hypothesis for each of the six tests. Since all *p*-values using the track-like neutrinos are larger than 0.5, the limits are set equal to the sensitivity in order to not 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.

## 5.3. Two-point Angular Correlation of UHECRs and High-energy Neutrinos

The 2pt-correlation analysis is applied to the sets of neutrino tracks and cascades separately, as well as for the background hypothesis of an isotropic UHECRs flux and an isotropic high-energy neutrino flux. The significance of the result with respect to the background hypothesis is corrected for the scan over the separation angles. This results in four *p*-values and four respective best-fit angular 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.

# 6. Conclusions

Three complementary analyses have been performed on the UHECR data sets provided by the Pierre Auger and the Telescope Array Collaborations combined with the high-energy and full-sky neutrino data sets provided by the IceCube and ANTARES Collaborations. For both the UHECR and neutrino data sets, the combination of data from the two respective observatories provides a field of view over the entire sky. None of the analyses found a result incompatible with the assumed background hypotheses of either an isotropic neutrino flux or an isotropic UHECR flux. This is an important update on the results reported in The IceCube Collaboration et al. (2016), where *p*-values close to the $3\sigma$ level were reported when applying analyses 2 and 3 to cascade-like neutrino events.

Based on the results, 90% C.L. upper limits are calculated for analysis 1 on the flux of neutrinos from point-like sources correlated with UHECRs and for analysis 2 on the number of UHECR events correlated with high-energy neutrinos. Analysis 1 reports upper limits on the correlated neutrino flux per
source given different sets of parameters that make assumptions about the lower cut on the UHECR energy and the fraction of UHECRs with a correlated neutrino source, as well as about a scaling factor of the magnetic deflection (see Table 2). Analysis 2 reports upper limits on the number of correlated UHECR events with either track-like or cascade-like neutrinos depending on the scaling factor of the magnetic deflection (see Table 3). This is the first time that upper limits from a direct correlation search of UHECRs and neutrinos are reported.

These limits are calculated based on the assumed signal and background hypotheses and are thus based on the parameters of the respective signal models. In the following, we discuss the limitations of these signal models in the light of underlying assumptions made for the neutrino–UHECR correlation. There are several plausible explanations why we have not observed a significant correlation of UHECRs and neutrinos.

Neutrino sources are presumably distributed over the whole universe, and the emitted neutrinos are able to reach Earth without deflection or absorption. For example, the first neutrino source identified with compelling evidence, TXS 0506+056 (Aartsen et al. 2018b, 2018c), is located at a redshift of $z = 0.34$ (Paiano et al. 2018), corresponding to a comoving distance of about 1.3 Gpc. This is beyond the assumed horizon for UHECR sources in the local universe of up to about 200 Mpc. Therefore, the fraction of correlated sources within the UHECR horizon is small compared to the total number of neutrino sources. This was quantified by Palladino et al. (2020) for muon neutrinos above 200 TeV based on the non-observation of neutrino multiplets above this energy. They concluded that the high-energy flux measured with IceCube must come from a large number of sources such that no multiplets are observed. Necessarily, most of the neutrino sources lie beyond the UHECR horizon, which is used in Palladino et al. (2020) to explain the lack of UHECR–neutrino correlations found in The IceCube Collaboration et al. (2016) and Schumacher (2019), and the same applies to the current study. In these studies, however, the energy threshold for neutrinos lies around 20 TeV for the combined high-energy data sets and as low as 100 GeV for the PS data set. The non-observation of UHECR-neutrino correlations thus extends to even lower energies than considered by Palladino et al. (2020). Even if there were local sources of both UHECRs and neutrinos, they could be transient phenomena, such that the UHECRs emitted over a short period of time arrive much later owing to their deflection, whereas the neutrinos travel on a straight path. It has been estimated by Davoudifar (2011) that the deflection in the EGMF causes a typical time delay on the order of $10^5$ yr considering a source distance of 50 Mpc and an EGMF strength of 2 nG. A delay of a couple of decades as expected from propagation in the GMF is already sufficient to decorrelate the observed UHECRs and neutrinos, as the live time of the neutrino and UHECR data sets ranges between 7 and 13 yr.

Another source of uncertainty is the mass composition and thus the charge of UHECRs at the highest energies. The measurements of the UHECR composition above ~50 EeV are not yet conclusive. However, several composition constraints at the highest energies (Aab et al. 2014a, 2014b) and the lack of significant, magnetically induced structures in the UHECR arrival directions (Aab et al. 2020a) suggest that less than 10% of UHECRs above ~50 EeV might be light elements. Only the light UHECRs are expected to have their source in the vicinity

of their arrival direction as derived from the benchmark deflection model. The sources of heavier UHECRs, e.g., iron nuclei, are spatially almost uncorrelated to the UHECR arrival directions owing to the 26 times larger magnetic deflection compared to protons. We see in analysis 1 that the limits on the neutrino flux are significantly relaxed when assuming a correlation fraction of 50% or 10%. This correlation fraction, defined as the number of UHECRs with a neutrino source in their vicinity, is a simplified model of the UHECR composition.

The modeling of the magnetic deflection as a 2D Gaussian as assumed for analyses 1 and 2 is a simplification with respect to the PT11 and JF12 GMF models (Pshirkov et al. 2011; Jansson & Farrar 2012), which themselves are subject to large uncertainties. Especially in the region of the Galactic plane, we expect deflections that are larger than the assumed mean value of around 3°. In addition, coherent deflection of UHECRs in the GMF and the deflection of UHECRs in the IGMF are not explicitly accounted for. Overall, a coherent shift of UHECRs depending on their arrival direction or an overall significantly larger deflection in the GMF and IGMF can dilute the spatial correlation of UHECRs and neutrinos. Nevertheless, a scaling of the overall deflection is covered partially with the scaling factor applied in analyses 1 and 2, since the overall strength of the deflection and the UHECR charge are largely degenerate parameters.

From a theoretical perspective, the connection of UHECRs and high-energy neutrinos is plausible based on their similar energy budget (Waxman & Bahcall 1999; Bahcall & Waxman 2001; Murase et al. 2013; Ahlers & Halzen 2018). However, this connection could not be proven with the current data and analysis approaches. It is to note that neutrinos with energies in the TeV–PeV range originate most likely from cosmic rays with energies below ~100 PeV (Murase & Waxman 2016), which is below the energy threshold of >50 EeV of the UHECR data sets used here. A direct connection to UHECRs can only be proven with ultra-high-energy neutrinos in the EeV range that have not been discovered yet. The non-observation of a UHECR-neutrino correlation rather points to the possibility that efficient UHECR sources are less efficient neutrino sources and vice versa: sources with efficient UHECR acceleration and emission require an optically thin proton and radiation environment, while sources with dense proton and radiation targets are efficient in neutrino production but not in UHECR emission (Murase et al. 2014; Rodrigues et al. 2021). This argument, however, can be relaxed when considering models, where cosmic rays below ~100 PeV are confined in a calorimetric environment and subsequently produce TeV–PeV neutrinos, while a fraction of the cosmic rays are accelerated to the highest energies and escape the source before interacting (Murase & Waxman 2016; Ahlers & Halzen 2018). Although no indication of such a scenario has been found in this analysis, the first indication of such a connection could be the observation of the Seyfert Galaxy NGC 1068 as a potential neutrino source candidate (Aartsen et al. 2020c; Inoue et al. 2020) and UHECR source candidate (Aab et al. 2018). As such sources are numerous also within the UHECR horizon, dedicated future searches correlating UHECR, photon, and neutrino observations might be able to set constraints on specific source candidates, instead of relying on UHECR and neutrino data alone. NGC 1068 could serve as blueprint for a
source given different sets of parameters that make assumptions about the lower cut on the UHECR energy and the fraction of UHECRs with a correlated neutrino source, as well as about a scaling factor of the magnetic deflection (see Table 2). Analysis 2 reports upper limits on the number of correlated UHECR events with either track-like or cascade-like neutrinos depending on the scaling factor of the magnetic deflection (see Table 3). This is the first time that upper limits from a direct correlation search of UHECRs and neutrinos are reported.

These limits are calculated based on the assumed signal and background hypotheses and are thus based on the parameters of the respective signal models. In the following, we discuss the limitations of these signal models in the light of underlying assumptions made for the neutrino–UHECR correlation. There are several plausible explanations why we have not observed a significant correlation of UHECRs and neutrinos.

Neutrino sources are presumably distributed over the whole universe, and the emitted neutrinos are able to reach Earth without deflection or absorption. For example, the first neutrino source identified with compelling evidence, TXS 0506+056 (Aartsen et al. 2018b, 2018c), is located at a redshift of $z = 0.34$ (Paiano et al. 2018), corresponding to a comoving distance of about 1.3 Gpc. This is beyond the assumed horizon for UHECR sources in the local universe of up to about 200 Mpc. Therefore, the fraction of correlated sources within the UHECR horizon is small compared to the total number of neutrino sources. This was quantified by Palladino et al. (2020) for muon neutrinos above 200 TeV based on the non-observation of neutrino multiplets above this energy. They concluded that the high-energy flux measured with IceCube must come from a large number of sources such that no multiplets are observed. Necessarily, most of the neutrino sources lie beyond the UHECR horizon, which is used in Palladino et al. (2020) to explain the lack of UHECR–neutrino correlations found in The IceCube Collaboration et al. (2016) and Schumacher (2019), and the same applies to the current study. In these studies, however, the energy threshold for neutrinos lies around 20 TeV for the combined high-energy data sets and as low as 100 GeV for the PS data set. The non-observation of UHECR-neutrino correlations thus extends to even lower energies than considered by Palladino et al. (2020). Even if there were local sources of both UHECRs and neutrinos, they could be transient phenomena, such that the UHECRs emitted over a short period of time arrive much later owing to their deflection, whereas the neutrinos travel on a straight path. It has been estimated by Davoudifar (2011) that the deflection in the EGMF causes a typical time delay on the order of $10^5$ yr considering a source distance of 50 Mpc and an EGMF strength of 2 nG. A delay of a couple of decades as expected from propagation in the GMF is already sufficient to decorrelate the observed UHECRs and neutrinos, as the live time of the neutrino and UHECR data sets ranges between 7 and 13 yr.

Another source of uncertainty is the mass composition and thus the charge of UHECRs at the highest energies. The measurements of the UHECR composition above ~50 EeV are not yet conclusive. However, several composition constraints at the highest energies (Aab et al. 2014a, 2014b) and the lack of significant, magnetically induced structures in the UHECR arrival directions (Aab et al. 2020a) suggest that less than 10% of UHECRs above ~50 EeV might be light elements. Only the light UHECRs are expected to have their source in the vicinity

of their arrival direction as derived from the benchmark deflection model. The sources of heavier UHECRs, e.g., iron nuclei, are spatially almost uncorrelated to the UHECR arrival directions owing to the 26 times larger magnetic deflection compared to protons. We see in analysis 1 that the limits on the neutrino flux are significantly relaxed when assuming a correlation fraction of 50% or 10%. This correlation fraction, defined as the number of UHECRs with a neutrino source in their vicinity, is a simplified model of the UHECR composition.

The modeling of the magnetic deflection as a 2D Gaussian as assumed for analyses 1 and 2 is a simplification with respect to the PT11 and JF12 GMF models (Pshirkov et al. 2011; Jansson & Farrar 2012), which themselves are subject to large uncertainties. Especially in the region of the Galactic plane, we expect deflections that are larger than the assumed mean value of around 3°. In addition, coherent deflection of UHECRs in the GMF and the deflection of UHECRs in the IGMF are not explicitly accounted for. Overall, a coherent shift of UHECRs depending on their arrival direction or an overall significantly larger deflection in the GMF and IGMF can dilute the spatial correlation of UHECRs and neutrinos. Nevertheless, a scaling of the overall deflection is covered partially with the scaling factor applied in analyses 1 and 2, since the overall strength of the deflection and the UHECR charge are largely degenerate parameters.

From a theoretical perspective, the connection of UHECRs and high-energy neutrinos is plausible based on their similar energy budget (Waxman & Bahcall 1999; Bahcall & Waxman 2001; Murase et al. 2013; Ahlers & Halzen 2018). However, this connection could not be proven with the current data and analysis approaches. It is to note that neutrinos with energies in the TeV–PeV range originate most likely from cosmic rays with energies below ~100 PeV (Murase & Waxman 2016), which is below the energy threshold of >50 EeV of the UHECR data sets used here. A direct connection to UHECRs can only be proven with ultra-high-energy neutrinos in the EeV range that have not been discovered yet. The non-observation of a UHECR-neutrino correlation rather points to the possibility that efficient UHECR sources are less efficient neutrino sources and vice versa: sources with efficient UHECR acceleration and emission require an optically thin proton and radiation environment, while sources with dense proton and radiation targets are efficient in neutrino production but not in UHECR emission (Murase et al. 2014; Rodrigues et al. 2021). This argument, however, can be relaxed when considering models, where cosmic rays below ~100 PeV are confined in a calorimetric environment and subsequently produce TeV–PeV neutrinos, while a fraction of the cosmic rays are accelerated to the highest energies and escape the source before interacting (Murase & Waxman 2016; Ahlers & Halzen 2018). Although no indication of such a scenario has been found in this analysis, the first indication of such a connection could be the observation of the Seyfert Galaxy NGC 1068 as a potential neutrino source candidate (Aartsen et al. 2020c; Inoue et al. 2020) and UHECR source candidate (Aab et al. 2018). As such sources are numerous also within the UHECR horizon, dedicated future searches correlating UHECR, photon, and neutrino observations might be able to set constraints on specific source candidates, instead of relying on UHECR and neutrino data alone. NGC 1068 could serve as blueprint for a
catalog of potential common sources to be constrained with
dedicated analyses.

In summary, the three analyses presented here reflect complementary approaches for tackling the question of a common origin of UHECRs and high-energy neutrinos. Despite substantially enlarged data sets and improved methods, the initially reported results in The IceCube Collaboration et al. (2016) could not be strengthened. For future analyses, we expect substantial gains in sensitivity if the charge of the UHECRs could be estimated on an event-by-event basis, as is expected for future measurements by Auger Prime (Aab et al. 2016).

The authors of the ANTARES collaboration acknowledge the financial support of the funding agencies: Centre National de la Recherche Scientifique (CNRS), Commissariat à l'énergie atomique et aux énergies alternatives (CEA), Commission Européenne (FEDER fund and Marie Curie Program), Institut Universitaire de France (IUF), LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001), Région Île-de-France (DIM-ACAV), Région Alsace (contrat CPER), Région Provence-Alpes-Côte d'Azur, Département du Var and Ville de La Seyne-sur-Mer, France; Bundesministerium für Bildung und Forschung (BMBF), Germany; Istituto Nazionale di Fisica Nucleare (INFN), Italy; Nederlandse organisatie voor Wetenschappelijk Onderzoek (NWO), the Netherlands; Council of the President of the Russian Federation for young scientists and leading scientific schools supporting grants, Russia; Executive Unit for Financing Higher Education, Research, Development and Innovation (UEFISCDI), Romania; Ministerio de Ciencia, Innovación, Investigación y Universidades (MCIU): Programa Estatal de Generación de Conocimiento (refs. PGC2018-096663-B-C41, -A-C42, -B-C43, -B-C44) (MCIU/FEDER), Generalitat Valenciana: Prometeo (PROMETEO/2020/019), Grisolía (refs. GRISOLIA/2018/119, /2021/192) and GenT (refs. CIDEGENT/2018/034, /2019/043, /2020/049, /2021/023) programs, Junta de Andalucía (ref. A-FQM-053-UGR18), La Caixa Foundation (ref. LCF/BQ/IN17/11620019), EU: MSC program (ref. 101025085), Spain; Ministry of Higher Education, Scientific Research and Innovation, Morocco, and the Arab Fund for Economic and Social Development, Kuwait. We also acknowledge the technical support of Ifremer, AIM and Foselev Marine for the sea operation and the CC-IN2P3 for the computing facilities. The ANTARES collaboration acknowledges the significant contributions to this manuscript from Julien Aublin.

The IceCube collaboration acknowledges the significant
contributions to this manuscript from Cyril Alispach, Anastasia
Barbano, and Lisa Schumacher. The authors of the IceCube
Collaboration acknowledge the support from the following
agencies: USA—U.S. National Science Foundation-Office of
Polar Programs, U.S. National Science Foundation-Physics
Division, U.S. National Science Foundation-EPSCoR, Wisconsin
Alumni Research Foundation, Center for High Throughput
Computing (CHTC) at the University of Wisconsin–Madison,
Open Science Grid (OSG), Extreme Science and Engineering
Discovery Environment (XSEDE), Frontera computing project at
the Texas Advanced Computing Center, U.S. Department of
Energy–National Energy Research Scientific Computing Center,
Particle astrophysics research computing center at the University
of Maryland, Institute for Cyber-Enabled Research at Michigan

State University, and Astroparticle physics computational facility at Marquette University; Belgium—Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany—Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden—Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; Australia—Australian Research Council; Canada—Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada; Denmark—Villum Fonden and Carlsberg Foundation; New Zealand—Marsden Fund; Japan—Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea—National Research Foundation of Korea (NRF); Switzerland—Swiss National Science Foundation (SNSF); United Kingdom—Department of Physics, University of Oxford. L.S. and C.H.W. acknowledge financial support from the “Deutsche Forschungsgemeinschaft”.

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:

Argentina—Comisión Nacional de Energía Atómica; Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT); Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET); Gobierno de la Provincia de Mendoza; Municipalidad de Malargüe; NDM Holdings and Valle Las Leñas; in gratitud for their continuing cooperation over land access; Australia—the Australian Research Council; Belgium—Fonds de la Recherche Scientifique (FNRS); Research Foundation Flanders (FWO); Brazil—Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq); Financiadora de Estudos e Projetos (FINEP); Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ); São Paulo Research Foundation (FAPESP) Grants No. 2019/10151-2, No. 2010/07359-6 and No. 1999/05404-3; Ministério da Ciência, Tecnologia, Inovações e Comunicações (MCTIC); Czech Republic—Grant No. MSMT CR LTT18004, LM2015038, LM2018102, CZ.02.1.01/0.0/0.0/16_013/0001402, CZ.02.1.01/0.0/0.0/18_046/0016010, and CZ.02.1.01/0.0/0.0/17_049/0008422; France—Centre de Calcul IN2P3/CNRS; Centre National de la Recherche Scientifique (CNRS); Conseil Régional Ile-de-France; Département Physique Nucléaire et Corpusculaire (PNC-IN2P3/CNRS); Département Sciences de l'Univers (SDU-INSU/CNRS); Institut Lagrange de Paris (ILP) Grant No. LABEX ANR-10-LABX-63 within the Investissements d'Avenir Programme Grant No. ANR-11-IDEX-0004-02; Germany—Bundesministerium für Bildung und Forschung (BMBF); Deutsche Forschungsgemeinschaft (DFG); Finanzministerium Baden-Württemberg; Helmholtz Alliance for Astroparticle Physics (HAP); Helmholtz-Gemeinschaft Deutscher Forschungszentren (HGF); Ministerium für Innovation, Wissenschaft und Forschung des Landes Nordrhein-Westfalen; Ministerium für Wissenschaft, Forschung und Kunst des Landes Baden-Württemberg; Italy—Istituto Nazionale di Fisica Nucleare
catalog of potential common sources to be constrained with
dedicated analyses.

In summary, the three analyses presented here reflect complementary approaches for tackling the question of a common origin of UHECRs and high-energy neutrinos. Despite substantially enlarged data sets and improved methods, the initially reported results in The IceCube Collaboration et al. (2016) could not be strengthened. For future analyses, we expect substantial gains in sensitivity if the charge of the UHECRs could be estimated on an event-by-event basis, as is expected for future measurements by Auger Prime (Aab et al. 2016).

The authors of the ANTARES collaboration acknowledge the financial support of the funding agencies: Centre National de la Recherche Scientifique (CNRS), Commissariat à l'énergie atomique et aux énergies alternatives (CEA), Commission Européenne (FEDER fund and Marie Curie Program), Institut Universitaire de France (IUF), LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001), Région Île-de-France (DIM-ACAV), Région Alsace (contrat CPER), Région Provence-Alpes-Côte d'Azur, Département du Var and Ville de La Seyne-sur-Mer, France; Bundesministerium für Bildung und Forschung (BMBF), Germany; Istituto Nazionale di Fisica Nucleare (INFN), Italy; Nederlandse organisatie voor Wetenschappelijk Onderzoek (NWO), the Netherlands; Council of the President of the Russian Federation for young scientists and leading scientific schools supporting grants, Russia; Executive Unit for Financing Higher Education, Research, Development and Innovation (UEFISCDI), Romania; Ministerio de Ciencia, Innovación, Investigación y Universidades (MCIU): Programa Estatal de Generación de Conocimiento (refs. PGC2018-096663-B-C41, -A-C42, -B-C43, -B-C44) (MCIU/FEDER), Generalitat Valenciana: Prometeo (PROMETEO/2020/019), Grisolía (refs. GRISOLIA/2018/119, /2021/192) and GenT (refs. CIDEGENT/2018/034, /2019/043, /2020/049, /2021/023) programs, Junta de Andalucía (ref. A-FQM-053-UGR18), La Caixa Foundation (ref. LCF/BQ/IN17/11620019), EU: MSC program (ref. 101025085), Spain; Ministry of Higher Education, Scientific Research and Innovation, Morocco, and the Arab Fund for Economic and Social Development, Kuwait. We also acknowledge the technical support of Ifremer, AIM and Foselev Marine for the sea operation and the CC-IN2P3 for the computing facilities. The ANTARES collaboration acknowledges the significant contributions to this manuscript from Julien Aublin.

The IceCube collaboration acknowledges the significant
contributions to this manuscript from Cyril Alispach, Anastasia
Barbano, and Lisa Schumacher. The authors of the IceCube
Collaboration acknowledge the support from the following
agencies: USA—U.S. National Science Foundation-Office of
Polar Programs, U.S. National Science Foundation-Physics
Division, U.S. National Science Foundation-EPSCoR, Wisconsin
Alumni Research Foundation, Center for High Throughput
Computing (CHTC) at the University of Wisconsin–Madison,
Open Science Grid (OSG), Extreme Science and Engineering
Discovery Environment (XSEDE), Frontera computing project at
the Texas Advanced Computing Center, U.S. Department of
Energy–National Energy Research Scientific Computing Center,
Particle astrophysics research computing center at the University
of Maryland, Institute for Cyber-Enabled Research at Michigan

State University, and Astroparticle physics computational facility at Marquette University; Belgium—Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany—Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden—Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; Australia—Australian Research Council; Canada—Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada; Denmark—Villum Fonden and Carlsberg Foundation; New Zealand—Marsden Fund; Japan—Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea—National Research Foundation of Korea (NRF); Switzerland—Swiss National Science Foundation (SNSF); United Kingdom—Department of Physics, University of Oxford. L.S. and C.H.W. acknowledge financial support from the “Deutsche Forschungsgemeinschaft”.

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:

Argentina—Comisión Nacional de Energía Atómica; Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT); Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET); Gobierno de la Provincia de Mendoza; Municipalidad de Malargüe; NDM Holdings and Valle Las Leñas; in gratitud for their continuing cooperation over land access; Australia—the Australian Research Council; Belgium—Fonds de la Recherche Scientifique (FNRS); Research Foundation Flanders (FWO); Brazil—Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq); Financiadora de Estudos e Projetos (FINEP); Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ); São Paulo Research Foundation (FAPESP) Grants No. 2019/10151-2, No. 2010/07359-6 and No. 1999/05404-3; Ministério da Ciência, Tecnologia, Inovações e Comunicações (MCTIC); Czech Republic—Grant No. MSMT CR LTT18004, LM2015038, LM2018102, CZ.02.1.01/0.0/0.0/16_013/0001402, CZ.02.1.01/0.0/0.0/18_046/0016010, and CZ.02.1.01/0.0/0.0/17_049/0008422; France—Centre de Calcul IN2P3/CNRS; Centre National de la Recherche Scientifique (CNRS); Conseil Régional Ile-de-France; Département Physique Nucléaire et Corpusculaire (PNC-IN2P3/CNRS); Département Sciences de l'Univers (SDU-INSU/CNRS); Institut Lagrange de Paris (ILP) Grant No. LABEX ANR-10-LABX-63 within the Investissements d'Avenir Programme Grant No. ANR-11-IDEX-0004-02; Germany—Bundesministerium für Bildung und Forschung (BMBF); Deutsche Forschungsgemeinschaft (DFG); Finanzministerium Baden-Württemberg; Helmholtz Alliance for Astroparticle Physics (HAP); Helmholtz-Gemeinschaft Deutscher Forschungszentren (HGF); Ministerium für Innovation, Wissenschaft und Forschung des Landes Nordrhein-Westfalen; Ministerium für Wissenschaft, Forschung und Kunst des Landes Baden-Württemberg; Italy—Istituto Nazionale di Fisica Nucleare
(INFN); Istituto Nazionale di Astrofisica (INAF); Ministero dell'Istruzione, dell'Università e della Ricerca (MIUR); CETEMPS Center of Excellence; Ministero degli Affari Esteri (MAE); México—Consejo Nacional de Ciencia y Tecnología (CONACYT) No. 167733; Universidad Nacional Autónoma de México (UNAM); PAPIIT DGAPA-UNAM; The Netherlands—Ministry of Education, Culture and Science; Netherlands Organisation for Scientific Research (NWO); Dutch national e-infrastructure with the support of SURF Cooperative; Poland—Ministry of Education and Science, grant No. DIR/WK/2018/11; National Science Centre, Grants No. 2016/22/M/ST9/00198, 2016/23/B/ST9/01635, and 2020/39/B/ST9/01398; Portugal—Portuguese national funds and FEDER funds within Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia (COMPETE); Romania—Ministry of Research, Innovation and Digitization, CNCS/CCCDI—UEFISCDI, projects PN19150201/16N/2019, PN1906010, TE128 and PED289, within PNCDI III; Slovenia—Slovenian Research Agency, grants P1-0031, P1-0385, IO-0033, N1-0111; Spain—Ministerio de Economía, Industria y Competitividad (FPA2017-85114-P and PID2019-104676GB-C32), Xunta de Galicia (ED431C 2017/07), Junta de Andalucía (SOMM17/6104/UGR, P18-FR-4314) Feder Funds, RENATA Red Nacional Temática de Astropartículas (FPA2015-68783-REDT) and María de Maeztu Unit of Excellence (MDM-2016-0692); USA—Department of Energy, Contracts No. DE-AC02-07CH11359, No. DE-FR02-04ER41300, No. DE-FG02-99ER41107, and No. DE-SC0011689; National Science Foundation, Grant No. 0450696; The Grainger Foundation; Marie Curie-IRSES/EPLANET; European Particle Physics Latin American Network; and UNESCO.

The Telescope Array experiment is supported by the Japan Society for the Promotion of Science (JSPS) through Grants-in-Aid for Priority Area 431, for Specially Promoted Research JP21000002, for Scientific Research (S) JP19104006, for Specially Promoted Research JP15H05693, for Scientific Research (S) JP15H05741, for Science Research (A) JP18H03705, for Young Scientists (A) JPH26707011, and for Fostering Joint International Research (B) JP19KK0074, by the joint research program of the Institute for Cosmic Ray Research (ICRR), The University of Tokyo; by the Pioneering Program of RIKEN for the Evolution of Matter in the Universe (r-EMU); by the U.S. National Science Foundation awards PHY-1404495, PHY-1404502, PHY-1607727, PHY-1712517, PHY-1806797, PHY-2012934, and PHY-2112904; by the National Research Foundation of Korea (2017K1A4A3015188, 2020R1A2C1008230, and 2020R1A2C2102800); by the Ministry of Science and Higher Education of the Russian Federation under the contract 075-15-2020-778, IISN project No. 4.4501.18; and by Belgian Science Policy under IUAP VII/37 (ULB). This work was partially supported by the grants of the joint research program of the Institute for Space-Earth Environmental Research, Nagoya University, and Inter-University Research Program of the Institute for Cosmic Ray Research of University of Tokyo. The foundations of Dr. Ezekiel R. and Edna Wattis Dumke, Willard L. Eccles, and George S. and Dolores Doré Eccles all helped with generous donations. The state of Utah supported the project through its Economic Development Board, and the University of Utah through the Office of the Vice President for Research. The experimental site became available through the cooperation of the Utah School and Institutional Trust Lands Administration (SITLA), U.S. Bureau of Land Management (BLM), and the

U.S. Air Force. We appreciate the assistance of the state of Utah and Fillmore offices of the BLM in crafting the Plan of Development for the site. Patrick A. Shea assisted the collaboration with valuable advice and supported the collaborations' efforts. The people and the officials of Millard County, Utah, have been a source of steadfast and warm support for our work, which we greatly appreciate. We are indebted to the Millard County Road Department for their efforts to maintain and clear the roads that get us to our sites. We gratefully acknowledge the contribution 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.

**Facility:** ANTARES, IceCube Neutrino Observatory, Pierre Auger Observatory, Telescope Array.

## ORCID iDs

F. Halzen https://orcid.org/0000-0001-6224-2417

## References

Aab, A., Abreu, P., Aglietta, M., et al. 2014a, PhRvD, 90, 122005

Aab, A., Abreu, P., Aglietta, M., et al. 2014b, PhRvD, 90, 122006

Aab, A., Abreu, P., Aglietta, M., et al. 2015, ApJ, 804, 15

Aab, A., Abreu, P., Aglietta, M., et al. 2016, arXiv:1604.03637

Aab, A., Abreu, P., Aglietta, M., et al. 2017, Sci, 357, 1266

Aab, A., Abreu, P., Aglietta, M., et al. 2018, ApJL, 853, L29

Aab, A., Abreu, P., Aglietta, M., et al. 2020a, JCAP, 2020, 017

Aab, A., Abreu, P., Aglietta, M., et al. 2020b, JInst, 15, P10021

Aartsen, M. G., Abbasi, R., Ackermann, M., et al. 2014a, JInst, 9, P03009

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2016, ApJ, 833, 3

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017a, ApJ, 835, 45

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017b, ApJ, 835, 151

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014b, ApJ, 796, 109

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014c, PhRvL, 113, 101101

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017c, ApH, 92, 30

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017d, JInst, 12, P03012

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018a, PhRvD, 98, 062003

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018b, Sci, 361, eaat1378

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018c, Sci, 361, 147

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, EPJC, 79, 234

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020a, ApJ, 892, 53

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020b, ApJL, 898, L10

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020c, PhRvL, 124, 051103

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020d, PhRvL, 125, 121104

Abbasi, R., Abdou, Y., Ackermann, M., et al. 2013, NIMPA, 700, 188

Abbasi, R., Ackermann, M., Adams, J., et al. 2021, PhRvD, 104, 022002

Abbasi, R. U., Abe, M., Abu-Zayyad, T., et al. 2014, ApJL, 790, L21

Abbasi, R. U., Abe, M., Abu-Zayyad, T., et al. 2018, ApJL, 867, L27

Abreu, P., Aglietta, M., Ahn, E. J., et al. 2011, arXiv:1107.4809

Abu-Zayyad, T., Aida, R., Allen, M., et al. 2012, NIMPA, 689, 87

Ageron, M., Aguilar, J. A., Al Samarai, I., et al. 2011, NIMPA, 656, 11

Aguilar, J. A., Albert, A., Ameli, F., et al. 2007, NIMPA, 570, 107

Ahlers, M., & Halzen, F. 2018, PrPN, 102, 73

Albert, A., André, M., Anghinolfi, M., et al. 2017, PhRvD, 96, 082001

Albert, A., André, M., Anghinolfi, M., et al. 2018, ApJL, 863, L30

Albert, A., André, M., Anghinolfi, M., et al. 2020, ApJ, 892, 92

Alves Batista, R., Biteau, J., Bustamante, M., et al. 2019, FrASS, 6, 23

Anchordoqui, L. A., Hooper, D., Sarkar, S., & Taylor, A. M. 2008, ApH, 29, 1

ANTARES Collaboration, Amram, P., Anghinolfi, M., et al. 2002, NIMPA, 484, 369

Aublin, J., Coleiro, A., Kouchner, A., et al. 2019, EPJ Web Conf., 210 (Les Ulis) in EDP Sciences, 03003

Bahcall, J., & Waxman, E. 2001, PhRvD, 64, 023002

Barbano, A. M. 2019, Proc. 36th Int. Cosmic Ray Conf.—PhS (ICRC2019), 358 (Trieste) in SISSA, 842

Biteau, J., Abreu, P., Aglietta, M., et al. 2022, Proc. 37th Int. Cosmic Ray Conf.—PhS (ICRC2021), 395 (Trieste: SISSA), 307

Biteau, J., Bister, T., Caccianiga, L., et al. 2019, EPJ Web Conf., 210 (Les Ulis: EDP Sciences), 01005

Bonifazi, C. & The Pierre Auger Collaboration 2009, NuPhS, 190, 20
(INFN); Istituto Nazionale di Astrofisica (INAF); Ministero dell'Istruzione, dell'Università e della Ricerca (MIUR); CETEMPS Center of Excellence; Ministero degli Affari Esteri (MAE); México—Consejo Nacional de Ciencia y Tecnología (CONACYT) No. 167733; Universidad Nacional Autónoma de México (UNAM); PAPIIT DGAPA-UNAM; The Netherlands—Ministry of Education, Culture and Science; Netherlands Organisation for Scientific Research (NWO); Dutch national e-infrastructure with the support of SURF Cooperative; Poland—Ministry of Education and Science, grant No. DIR/WK/2018/11; National Science Centre, Grants No. 2016/22/M/ST9/00198, 2016/23/B/ST9/01635, and 2020/39/B/ST9/01398; Portugal—Portuguese national funds and FEDER funds within Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia (COMPETE); Romania—Ministry of Research, Innovation and Digitization, CNCS/CCCDI—UEFISCDI, projects PN19150201/16N/2019, PN1906010, TE128 and PED289, within PNCDI III; Slovenia—Slovenian Research Agency, grants P1-0031, P1-0385, I0-0033, N1-0111; Spain—Ministerio de Economía, Industria y Competitividad (FPA2017-85114-P and PID2019-104676GB-C32), Xunta de Galicia (ED431C 2017/07), Junta de Andalucía (SOMM17/6104/UGR, P18-FR-4314) Feder Funds, RENATA Red Nacional Temática de Astropartículas (FPA2015-68783-REDT) and María de Maeztu Unit of Excellence (MDM-2016-0692); USA—Department of Energy, Contracts No. DE-AC02-07CH11359, No. DE-FR02-04ER41300, No. DE-FG02-99ER41107, and No. DE-SC0011689; National Science Foundation, Grant No. 0450696; The Grainger Foundation; Marie Curie-IRSES/EPLANET; European Particle Physics Latin American Network; and UNESCO.

The Telescope Array experiment is supported by the Japan Society for the Promotion of Science (JSPS) through Grants-in-Aid for Priority Area 431, for Specially Promoted Research JP21000002, for Scientific Research (S) JP19104006, for Specially Promoted Research JP15H05693, for Scientific Research (S) JP15H05741, for Science Research (A) JP18H03705, for Young Scientists (A) JPH26707011, and for Fostering Joint International Research (B) JP19KK0074, by the joint research program of the Institute for Cosmic Ray Research (ICRR), The University of Tokyo; by the Pioneering Program of RIKEN for the Evolution of Matter in the Universe (r-EMU); by the U.S. National Science Foundation awards PHY-1404495, PHY-1404502, PHY-1607727, PHY-1712517, PHY-1806797, PHY-2012934, and PHY-2112904; by the National Research Foundation of Korea (2017K1A4A3015188, 2020R1A2C1008230, and 2020R1A2C2102800); by the Ministry of Science and Higher Education of the Russian Federation under the contract 075-15-2020-778, IISN project No. 4.4501.18; and by Belgian Science Policy under IUAP VII/37 (ULB). This work was partially supported by the grants of the joint research program of the Institute for Space-Earth Environmental Research, Nagoya University, and Inter-University Research Program of the Institute for Cosmic Ray Research of University of Tokyo. The foundations of Dr. Ezekiel R. and Edna Wattis Dumke, Willard L. Eccles, and George S. and Dolores Doré Eccles all helped with generous donations. The state of Utah supported the project through its Economic Development Board, and the University of Utah through the Office of the Vice President for Research. The experimental site became available through the cooperation of the Utah School and Institutional Trust Lands Administration (SITLA), U.S. Bureau of Land Management (BLM), and the

U.S. Air Force. We appreciate the assistance of the state of Utah and Fillmore offices of the BLM in crafting the Plan of Development for the site. Patrick A. Shea assisted the collaboration with valuable advice and supported the collaborations' efforts. The people and the officials of Millard County, Utah, have been a source of steadfast and warm support for our work, which we greatly appreciate. We are indebted to the Millard County Road Department for their efforts to maintain and clear the roads that get us to our sites. We gratefully acknowledge the contribution 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.

*Facility:* ANTARES, IceCube Neutrino Observatory, Pierre Auger Observatory, Telescope Array.

## ORCID iDs

F. Halzen https://orcid.org/0000-0001-6224-2417

## References

Aab, A., Abreu, P., Aglietta, M., et al. 2014a, PhRvD, 90, 122005

Aab, A., Abreu, P., Aglietta, M., et al. 2014b, PhRvD, 90, 122006

Aab, A., Abreu, P., Aglietta, M., et al. 2015, ApJ, 804, 15

Aab, A., Abreu, P., Aglietta, M., et al. 2016, arXiv:1604.03637

Aab, A., Abreu, P., Aglietta, M., et al. 2017, Sci, 357, 1266

Aab, A., Abreu, P., Aglietta, M., et al. 2018, ApJL, 853, L29

Aab, A., Abreu, P., Aglietta, M., et al. 2020a, JCAP, 2020, 017

Aab, A., Abreu, P., Aglietta, M., et al. 2020b, JInst, 15, P10021

Aartsen, M. G., Abbasi, R., Ackermann, M., et al. 2014a, JInst, 9, P03009

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2016, ApJ, 833, 3

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017a, ApJ, 835, 45

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017b, ApJ, 835, 151

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014b, ApJ, 796, 109

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014c, PhRvL, 113, 101101

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017c, ApH, 92, 30

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017d, JInst, 12, P03012

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018a, PhRvD, 98, 062003

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018b, Sci, 361, eaat1378

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018c, Sci, 361, 147

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, EPJC, 79, 234

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020a, ApJ, 892, 53

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020b, ApJL, 898, L10

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020c, PhRvL, 124, 051103

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020d, PhRvL, 125, 121104

Abbasi, R., Abdou, Y., Ackermann, M., et al. 2013, NIMPA, 700, 188

Abbasi, R., Ackermann, M., Adams, J., et al. 2021, PhRvD, 104, 022002

Abbasi, R. U., Abe, M., Abu-Zayyad, T., et al. 2014, ApJL, 790, L21

Abbasi, R. U., Abe, M., Abu-Zayyad, T., et al. 2018, ApJL, 867, L27

Abreu, P., Aglietta, M., Ahn, E. J., et al. 2011, arXiv:1107.4809

Abu-Zayyad, T., Aida, R., Allen, M., et al. 2012, NIMPA, 689, 87

Ageron, M., Aguilar, J. A., Al Samarai, I., et al. 2011, NIMPA, 656, 11

Aguilar, J. A., Albert, A., Ameli, F., et al. 2007, NIMPA, 570, 107

Ahlers, M., & Halzen, F. 2018, PrPN, 102, 73

Albert, A., André, M., Anghinolfi, M., et al. 2017, PhRvD, 96, 082001

Albert, A., André, M., Anghinolfi, M., et al. 2018, ApJL, 863, L30

Albert, A., André, M., Anghinolfi, M., et al. 2020, ApJ, 892, 92

Alves Batista, R., Biteau, J., Bustamante, M., et al. 2019, FrASS, 6, 23

Anchordoqui, L. A., Hooper, D., Sarkar, S., & Taylor, A. M. 2008, ApH, 29, 1

ANTARES Collaboration, Amram, P., Anghinolfi, M., et al. 2002, NIMPA, 484, 369

Aublin, J., Coleiro, A., Kouchner, A., et al. 2019, EPJ Web Conf., 210 (Les Ulis) in EDP Sciences, 03003

Bahcall, J., & Waxman, E. 2001, PhRvD, 64, 023002

Barbano, A. M. 2019, Proc. 36th Int. Cosmic Ray Conf.—PoS (ICRC2019), 358 (Trieste) in SISSA, 842

Biteau, J., Abreu, P., Aglietta, M., et al. 2022, Proc. 37th Int. Cosmic Ray Conf.—PoS (ICRC2021), 395 (Trieste: SISSA), 307

Biteau, J., Bister, T., Caccianiga, L., et al. 2019, EPJ Web Conf., 210 (Les Ulis: EDP Sciences), 01005

Bonifazi, C. & The Pierre Auger Collaboration 2009, NuPhS, 190, 20
Davoudifar, P. 2011, 32nd Int. Cosmic Ray Conf., 2 (Trieste: SISSA), 230

Dawson, B. 2019, Proc. 36th Int. Cosmic Ray Conf.—PoS (ICRC2019), 358 (Trieste) in SISSA, 231

Durrer, R., & Neronov, A. 2013, A&ARv, 21, 62

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Haack, C., & Wiebusch, C. 2017, Proc. 35th Int. Cosmic Ray Conf.—PoS (ICRC2017), 301 (Trieste: SISSA), 1005

Illuminati, G., Aublin, J., & Navas, S. 2019, Proc. 36th Int. Cosmic Ray Conf. —PoS (ICRC2019), 358 (Trieste) in SISSA, 920

Inoue, Y., Khangulyan, D., & Doi, A. 2020, ApJ, 891, L33

Jansson, R., & Farrar, G. R. 2012, ApJ, 757, 14

Kawai, H., et al. 2008, NuPhS, 175, 221

Kopper, C. 2017, Proc. 35th Int. Cosmic Ray Conf.—PoS (ICRC2017), 301 (Trieste) in SISSA, 981

Kronberg, P. P. 1994, RPPh, 57, 325

Mészáros, P. 2018, PhT, 71, 36

Murase, K. 2015, AIP Conf. Proc. 1666 (Melville, NY) in AIP, 040006

Murase, K., Ahlers, M., & Lacki, B. C. 2013, PhRvD, 88, 121301

Murase, K., Inoue, Y., & Dermer, C. D. 2014, PhRvD, 90, 023007

Murase, K., & Waxman, E. 2016, PhRvD, 94, 103006

Nagano, M., & Watson, A. A. 2000, RvMP, 72, 689

Paiano, S., Falomo, R., Treves, A., & Scarpa, R. 2018, ApJL, 854, L32

Palladino, A., van Vliet, A., Winter, W., & Franckowiak, A. 2020, MNRAS, 494, 4255

Particle Data Group, Zyla, P. A., Barnett, R. M., et al. 2020, PTEP, 2020, 083C01

Pshirkov, M. S., Tinyakov, P. G., Kronberg, P. P., & Newton-McGee, K. J. 2011, ApJ, 738, 192

Pshirkov, M. S., Tinyakov, P. G., & Urban, F. R. 2016, PhRvL, 116, 191302

Rodrigues, X., Heinze, J., Palladino, A., van Vliet, A., & Winter, W. 2021, PhRvL, 126, 191101

Schumacher, L. 2019, EPJ Web Conf., 207 (Les Ulis: EDP Sciences), 02010

Sommers, P. 2001, APh, 14, 271

Stettner, J. 2019, Proc. 36th Int. Cosmic Ray Conf.—PoS (ICRC2019), 358 (Trieste) in SISSA, 1017

The IceCube Collaboration 2013, Sci, 342, 1242856

The IceCube Collaboration The Pierre Auger Collaboration & The Telescope Array Collaboration 2016, JCAP, 2016, 037

The Pierre Auger Collaboration 2015, NIMPA, 798, 172

Tkachev, I., Fujii, T., Ivanov, D., et al. 2021, Proc. 37th Int. Cosmic Ray Conf. —PoS (ICRC2021), 395 (Trieste) in SISSA, 392

Tokuno, H., Tameda, Y., Takeda, M., et al. 2012, NIMPA, 676, 54

Troitsky, S., Fukushima, M., Ikeda, D., et al. 2017, Proc. 35th Int. Cosmic Ray Conf.—PoS (ICRC2017), 301 (Trieste) in SISSA, 548

Wandkowsky, N. 2018, XXVIII Int. Conf. on Neutrino Physics and Astrophysics, 501

Waxman, E., & Bahcall, J. N. 1999, PhRvD, 59, 023002
Davoudifar, P. 2011, 32nd Int. Cosmic Ray Conf., 2 (Trieste: SISSA), 230

Dawson, B. 2019, Proc. 36th Int. Cosmic Ray Conf.—PoS (ICRC2019), 358 (Trieste) in SISSA, 231

Durrer, R., & Neronov, A. 2013, A&ARv, 21, 62

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Haack, C., & Wiebusch, C. 2017, Proc. 35th Int. Cosmic Ray Conf.—PoS (ICRC2017), 301 (Trieste: SISSA), 1005

Illuminati, G., Aublin, J., & Navas, S. 2019, Proc. 36th Int. Cosmic Ray Conf. —PoS (ICRC2019), 358 (Trieste) in SISSA, 920

Inoue, Y., Khangulyan, D., & Doi, A. 2020, ApJ, 891, L33

Jansson, R., & Farrar, G. R. 2012, ApJ, 757, 14

Kawai, H., et al. 2008, NuPhS, 175, 221

Kopper, C. 2017, Proc. 35th Int. Cosmic Ray Conf.—PoS (ICRC2017), 301 (Trieste) in SISSA, 981

Kronberg, P. P. 1994, RPPh, 57, 325

Mészáros, P. 2018, PhT, 71, 36

Murase, K. 2015, AIP Conf. Proc. 1666 (Melville, NY) in AIP, 040006

Murase, K., Ahlers, M., & Lacki, B. C. 2013, PhRvD, 88, 121301

Murase, K., Inoue, Y., & Dermer, C. D. 2014, PhRvD, 90, 023007

Murase, K., & Waxman, E. 2016, PhRvD, 94, 103006

Nagano, M., & Watson, A. A. 2000, RvMP, 72, 689

Paiano, S., Falomo, R., Treves, A., & Scarpa, R. 2018, ApJL, 854, L32

Palladino, A., van Vliet, A., Winter, W., & Franckowiak, A. 2020, MNRAS, 494, 4255

Particle Data Group, Zyla, P. A., Barnett, R. M., et al. 2020, PTEP, 2020, 083C01

Pshirkov, M. S., Tinyakov, P. G., Kronberg, P. P., & Newton-McGee, K. J. 2011, ApJ, 738, 192

Pshirkov, M. S., Tinyakov, P. G., & Urban, F. R. 2016, PhRvL, 116, 191302

Rodrigues, X., Heinze, J., Palladino, A., van Vliet, A., & Winter, W. 2021, PhRvL, 126, 191101

Schumacher, L. 2019, EPJ Web Conf., 207 (Les Ulis: EDP Sciences), 02010

Sommers, P. 2001, APh, 14, 271

Stettner, J. 2019, Proc. 36th Int. Cosmic Ray Conf.—PoS (ICRC2019), 358 (Trieste) in SISSA, 1017

The IceCube Collaboration 2013, Sci, 342, 1242856

The IceCube Collaboration The Pierre Auger Collaboration & The Telescope Array Collaboration 2016, JCAP, 2016, 037

The Pierre Auger Collaboration 2015, NIMPA, 798, 172

Tkachev, I., Fujii, T., Ivanov, D., et al. 2021, Proc. 37th Int. Cosmic Ray Conf. —PoS (ICRC2021), 395 (Trieste) in SISSA, 392

Tokuno, H., Tameda, Y., Takeda, M., et al. 2012, NIMPA, 676, 54

Troitsky, S., Fukushima, M., Ikeda, D., et al. 2017, Proc. 35th Int. Cosmic Ray Conf.—PoS (ICRC2017), 301 (Trieste) in SISSA, 548

Wandkowsky, N. 2018, XXVIII Int. Conf. on Neutrino Physics and Astrophysics, 501

Waxman, E., & Bahcall, J. N. 1999, PhRvD, 59, 023002