(Aartsen et al., 2017d) and data from 3.5 years of operation between 2015 and
2018 that was selected for the real-time gamma-ray follow-up (GFU) program of
IceCube (Aartsen et al., 2018b, 2017c). The combined data set consists of
about $1.4$ million track-like events above
$\sim$100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$. It is dominated by
atmospheric neutrinos in the Northern hemisphere and by atmospheric muons in
the Southern hemisphere. The median of the angular resolution ($\Psi$) is
better than $0.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ above energies of a few
TeV. Figure 2 shows the median of the angular resolution for the different
detector configurations of IceCube, and for data from the ANTARES detector
(see subsection 2.2). Note that analysis 1 is not based on the median of the
angular resolution, but on the estimator for the angular resolution on an
event-by-event basis ($\sigma$).

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\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ and
$2\text{\,}\mathrm{P}\mathrm{e}\mathrm{V}$. Integrated above
$60\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, which corresponds to 60 events in
total, the percentage of events of astrophysical origin, i.e. the
astrophysical purity, is larger than $75\text{\,}\mathrm{\char 37\relax}$
(Wandkowsky, 2018), while the percentage is lower below
$60\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$. The angular resolution is about
$1\text{\,}\mathrm{\SIUnitSymbolDegree}$ for track-like events and
$15\text{\,}\mathrm{\SIUnitSymbolDegree}$ for cascade-like events above
$100\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$. The resolution of the deposited
energy for tracks and cascades is around $10\text{\,}\mathrm{\char 37\relax}$
(Aartsen et al., 2014a) without accounting for systematic uncertainties, but
the cascades have a better correlation with the primary neutrino energy since
they deposit most of their energy inside the detector, while tracks do not.

The HENU data set consists of the 35 highest energy track-like events with a
reconstructed declination $\geq-5{}^{\circ}$, 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 $\sim$100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$, only events of
high probability of non-atmospheric origin have been selected by applying an
energy threshold of $\geq$200\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$ on the
reconstructed muon energy. This corresponds to an astrophysical purity of more
than $50\text{\,}\mathrm{\char 37\relax}$. 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) ${{\rm d}\phi/{\rm
d}E}=$1.01\text{\times}{10}^{-18}\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}^{-1}\mathrm{c}\mathrm{m}^{-2}\mathrm{s}^{-1}\mathrm{s}\mathrm{r}^{-1}$(E/$100\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$)^{-2.19}$
is used as well as the best fit atmospheric flux from the same reference.

The EHE data set consists of 20 events that have been collected between 2008
and 2017 (Aartsen et al., 2018c, 2017c). The selection is targeting high-
energy track-like events of good angular resolution
$\leq$1\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The selection has been
optimized to be sensitive to events in the energy range of
$0.5\text{\,}\mathrm{P}\mathrm{e}\mathrm{V}10.0\text{\,}\mathrm{P}\mathrm{e}\mathrm{V}$.
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\text{\,}\mathrm{\char 37\relax}$ purity (see Table 2 in
Aartsen et al. (2017c)).

###  2.2 The ANTARES Neutrino Telescope

Figure 2: Median of the angular resolution, $\Psi$, which is the angle between the true neutrino direction and reconstructed muon direction for IceCube and ANTARES point-source data as a function of the true neutrino energy. This is calculated based on simulation data sets. The data from partial detector configurations of IceCube are denoted by the number of operating strings, (IC40, IC59, IC79; Aartsen et al., 2014c) and are shown in blue. The data of the full detector configuration with 86 strings including the GFU data set (Aartsen et al., 2017c, d, 2018b) are identified by the years of data taking and are shown in orange. The angular resolution of ANTARES data for the 11-year data set (Albert et al., 2018; Illuminati et al., 2019) is shown in black. Detector | Analysis 1 | Analysis 2 | Analysis 3
---|---|---|---
ANTARES | PS | HENU | HENU
IceCube | PS | HESE + HENU + EHE | HESE + HENU + EHE
Data set | Description | Topology
PS | Optimized for point-source searches, $\nu_{\mu}$ candidates. | Tracks
HESE | High-energy starting events, all flavors. | Tracks and cascades
HENU | High-energy selection of $\nu_{\mu}$ candidates. | Tracks
EHE | Extremely-high energy $\nu_{\mu}$ candidates. | Tracks
Table 1: Overview over different neutrino data sets used in the different
analyses.

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\text{\,}\mathrm{m}$. The strings are spaced at distances of
about $70\text{\,}\mathrm{m}$ from each other, instrumenting a total volume of
$\sim$$0.01\text{\,}\mathrm{k}\mathrm{m}^{3}$. Each string is equipped with 25
storeys of 3 optical modules (ANTARES Collaboration et al., 2002), vertically
spaced by $14.5\text{\,}\mathrm{m}$, for a total of 885 optical modules. Each
optical module houses a 10-inch photomultiplier tube facing 45∘ downward to
optimize the detection of light from upward going charged particles. The
detector was completed in 2008.

The ANTARES and IceCube detection principles are very similar (see subsection
2.1). Particles above the Cherenkov threshold induce coherent radiation
emitted in a cone with a characteristic angle of 42∘ in water. 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 40K
decays and bioluminescence (Aguilar et al., 2007). ANTARES is thus also
capable of detecting neutrino charged- and neutral-current interactions of all
flavors.

For analysis 1, we use all track-like events from the 11-year data set used
for point-source searches (PS) recorded between 2007 and 2017 (Albert et al.,
2018; Illuminati et al., 2019). The high-energy events (HENU) for analyses 2
and 3 are selected from an earlier data set of track-like and cascade-like
events collected between 2007 and 2015 (Albert et al., 2017). In order to
ensure a high probability of astrophysical origin for analyses 2 and 3, we
require an astrophysical purity $\geq$40\text{\,}\%$$ 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
$0.4\text{\,}\mathrm{\SIUnitSymbolDegree}$ above
$10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ (Albert et al., 2017). Overall,
the angular resolution of ANTARES tracks is better than for IceCube tracks for
energies below $100\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$, and comparable
around $100\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ and above. The median
angular resolution, $\Psi$, of the 11-year data set compared to the IceCube PS
data set is shown in Figure 2.

Despite the smaller detection volume, the inclusion of ANTARES data
significantly improves the all-sky coverage of the neutrino data set used in
analysis 1. This is shown in subsection 4.1 and Figure 5 through the relative
contribution to the expected number of signal events for the individual PS
data sets of ANTARES and IceCube. Particularly in the Southern hemisphere,
where the background from atmospheric muons results in a higher energy
threshold in IceCube, the ANTARES data contributes substantially to the signal
acceptance for soft source spectra.

###  2.3 The Pierre Auger Observatory

The Pierre Auger Observatory is located in Argentina
($32\text{\,}\mathrm{\SIUnitSymbolDegree}\mathrm{S}$,
$69\text{\,}\mathrm{\SIUnitSymbolDegree}\mathrm{W}$) at a mean altitude of
$1.4\text{\,}\mathrm{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\text{\,}\mathrm{k}\mathrm{m}^{2}$. The SD area is
overlooked by the FD, which consists of 27 wide-angle optical telescopes
located at four sites. The reconstruction of SD events is described in detail
in Aab et al. (2020). The energy estimate is based on the signal at
$1000\text{\,}\mathrm{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\text{\,}\mathrm{\char 37\relax}$ (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\text{\,}\mathrm{\char 37\relax}$ (Abreu et al., 2011; Aab et al., 2020),
and the angular uncertainty is less than
$0.9\text{\,}\mathrm{\SIUnitSymbolDegree}$ (Bonifazi & The Pierre Auger
Collaboration, 2009).

The UHECR data set used here consists of 324 events observed with the SD
between January 2004 and April 2017 (Aab et al., 2018) with reconstructed
energies $\geq$52\text{\,}\mathrm{E}\mathrm{e}\mathrm{V}$$ and zenith angles
$\leq$80\text{\,}\mathrm{\SIUnitSymbolDegree}$$. 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\text{\,}\mathrm{\SIUnitSymbolDegree}45\text{\,}\mathrm{\SIUnitSymbolDegree}$
in declination. All energies are scaled up by a constant factor of
$14\text{\,}\mathrm{\char 37\relax}$ to match both Auger and TA energies on a
common energy scale (see subsection 2.4), following the recommendation of the
Auger-TA joint working group (Biteau et al., 2019). This scaling factor was
determined by cross-calibrating the measured UHECR fluxes in the declination
band around the celestial equator covered by both observatories. It is chosen
such that the fluxes in the common declination band match at
$40\text{\,}\mathrm{EeV}$ of the Auger data set and at
$53\text{\,}\mathrm{EeV}$ of the TA data set, respectively. The celestial
distribution of the UHECR arrival directions is shown in Figure 1. The
directional exposure as function of declination can be found in Figure 3,
which amounts to an integrated geometric exposure of
$91\,300\text{\,}\mathrm{k}\mathrm{m}^{2}\,\mathrm{y}\mathrm{r}\,\mathrm{s}\mathrm{r}$.

###  2.4 The Telescope Array

Figure 3: Directional exposures of TA, Auger, and their sum. The underlying
geometric exposure functions are derived from Sommers (2001). See also Fig. 1
of Biteau et al. (2019).

The Telescope Array (TA) is located in Millard County, Utah, USA
($39.3\text{\,}\mathrm{\SIUnitSymbolDegree}\mathrm{N}$,
$112.9\text{\,}\mathrm{\SIUnitSymbolDegree}\mathrm{W}$) at an altitude of
about $1.4\text{\,}\mathrm{km}$ (Kawai et al., 2008). It consists of a surface
detector (SD) array, composed of $507$ plastic scintillation detectors of
$3\text{\,}\mathrm{m}^{2}$ each. The SD stations are located on a square grid
with $1.2\text{\,}\mathrm{km}$ separation, which extends over an area of
$700\text{\,}\mathrm{k}\mathrm{m}^{2}$ (Abu-Zayyad et al., 2012). The TA SD
array observes UHECRs with a duty cycle near $100\text{\,}\mathrm{\char
37\relax}$. With its wide field of view, the SD array covers a range from
$-15\text{\,}\mathrm{\SIUnitSymbolDegree}90\text{\,}\mathrm{\SIUnitSymbolDegree}$
in declination. In addition to the SD, there are three fluorescence telescope
stations, instrumented with $1214$ telescopes each (Tokuno et al., 2012). The
telescope stations observe the sky above the SD array and measure the
longitudinal development of the air showers as they traverse the atmosphere.

The data set for this analysis follows the selection in Abbasi et al. (2014),
which has been updated to 9 years of data in Abbasi et al. (2018). The data
set is identical to the data set used for anisotropy analyses presented in
Troitsky et al. (2017). It consists of 143 events observed with the SD between
May 2008 and May 2017 with reconstructed energies
$\geq$57\text{\,}\mathrm{E}\mathrm{e}\mathrm{V}$$ and zenith angles
$\leq$55\text{\,}\mathrm{\SIUnitSymbolDegree}$$. At these energies, the
angular uncertainty is about $1.5\text{\,}\mathrm{\SIUnitSymbolDegree}$. The
statistical uncertainty on the reconstructed energy is about
$15\text{\,}\mathrm{\char 37\relax}20\text{\,}\mathrm{\char 37\relax}$, with
an additional systematic uncertainty on the energy scale of
$21\text{\,}\mathrm{\char 37\relax}$ (Abbasi et al., 2018). All energies are
scaled down by $-14\text{\,}\mathrm{\char 37\relax}$ to match both Auger and
TA energies on a common energy scale (see subsection 2.3), following the
recommendation of the Auger–TA joint working group (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 declination of all UHECR
events, showing that TA data substantially contributes to the full-sky
exposure of the combined UHECR data set. Figure 3 shows the directional
exposure as a function of declination, with an integrated total exposure of
$11\,600\text{\,}\mathrm{k}\mathrm{m}^{2}\,\mathrm{y}\mathrm{r}\,\mathrm{s}\mathrm{r}$
(Biteau et al., 2019).

##  3 Magnetic deflections

Figure 4: Deflection simulation for protons of $100$ EeV and two different
galactic magnetic models: PT2011 (Pshirkov et al., 2011) and JF2012 (Jansson &
Farrar, 2012). Shaded areas in blue and orange show the histogram split into
Galactic Northern and Southern hemispheres, with their sum shown as black
line. Additionally, Rayleigh PDFs are shown in the same colors; a Rayleigh
function with mode $\sigma$ is the 1D projection of a symmetric 2D-Gaussian
function with width $\sigma$, with the projection being
$(x,y)\rightarrow\sqrt{x^{2}+y^{2}}$. The position of the markers indicate the
mode of the Rayleigh PDF and thus the value of the default deflection, $D$.
Note that the integrals of the split PDFs are normalized to 0.5, while the
integrals of the full PDFs are normalized to 1.

Despite their extremely high rigidity, UHECRs are deflected in galactic and in
extragalactic magnetic fields (GMF/EGMF) by a non-negligible amount. Neither
the strength and correlation length of the extragalactic magnetic fields
(Durrer & Neronov, 2013; Kronberg, 1994) nor the distance of the UHECR sources
are known well. Measurements of the Faraday rotation of extragalactic sources
indicate that the extragalactic magnetic fields are weaker than
$1\text{\,}\mathrm{n}\mathrm{G}$ (Pshirkov et al., 2016). Assuming a
correlation length of $\sim$$1\text{\,}\mathrm{Mpc}$, this results in a
deflection less than $2\text{\,}\mathrm{\SIUnitSymbolDegree}$ for protons of
$100\text{\,}\mathrm{EeV}$, even at source distances of
$50\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$. In line with the previously
reported results (The IceCube Collaboration et al., 2016), the deflection
outside of our galaxy is assumed to be generally weaker than within our
galaxy, and it is not modeled explicitly but benchmarked within the
uncertainty of the deflection by GMF and the uncertainty of the rigidity due
to the unknown composition.

The measurement and modeling of the GMF is a complex task and subject of on-
going discussion111See https://icrc2021-venue.desy.de for a recent review by
Tess Jaffe on the GMF at ICRC2021:
”Constraining Magnetic Fields at Galactic Scales”.. Among different proposed
models, we use the JF2012 (Jansson & Farrar, 2012) model and the PT2011
(Pshirkov et al., 2011) model to estimate the deflection of the UHECRs in the
GMF. Both models consist of a disc and a halo component, while the JF2012
model has an additional x-shaped field component perpendicular to the disc. A
propagation simulation has been conducted with a Monte-Carlo approach to
estimate the deflection of protons with an energy of
$100\text{\,}\mathrm{EeV}$ that are distributed isotropically outside of the
GMF. The resulting deflections are shown in Figure 4, split into the Galactic
Northern and Southern hemisphere based on the arrival direction of the proton.
The estimated deflection shows a complex structure which differs considerably
for the two models. However, the mean deflection is about
$3\text{\,}\mathrm{\SIUnitSymbolDegree}$ in both cases. Due to the heavy tails
of the distributions, the mean is consistently larger than the median, thus we
choose the mean as a conservative deflection estimate. The split into the two
hemispheres shows a considerably smaller deflection in the Galactic North than
in the Galactic South due to North-South asymmetries present in both models of
the GMF. For this work we have chosen a robust benchmark modeling of the
effect of the GMF, thus avoiding biases of detailed model uncertainties. The
deflection process is assumed to be random, resulting in a symmetric 2D
Gaussian distribution of UHECR arrival directions around a given source
direction. In reverse, the source position is assumed to be located within the
2D Gaussian distribution around the arrival direction of a UHECR event. The
standard deviation of the Gaussian deflection, $\sigma_{\rm MD}$, depends on a
scaling factor, $C$, and inversely on the energy of the UHECR event, $E_{\rm
CR}$

$\sigma_{\rm MD}=\frac{D_{0}\cdot C}{E_{\rm
CR}/$100\text{\,}\mathrm{E}\mathrm{e}\mathrm{V}$}.$ (1)

The default deflection, $D_{0}$, for $C=1$ and $E_{\rm
CR}=$100\text{\,}\mathrm{E}\mathrm{e}\mathrm{V}$$ is derived from the mean of
the deflection values obtained with the simulation, which are shown in Figure
4. Note that the deflection is usually a 2D quantity with $(x,y)$ coordinates,
which in Figure 4 is projected to the absolute value,
i.e.$\sqrt{x^{2}+y^{2}}$. Furthermore, larger and thus more conservative
values of the deflection will be tested to include an uncertainty of the GMF
model. The scaling factor $C$ thus accounts for uncertainties of both GMF
model and UHECR charge, which is not known on an event-by-event basis at
highest energies.

Analysis 1 uses a default deflection of
$D_{0}=$3\text{\,}\mathrm{\SIUnitSymbolDegree}$$ over the whole sky, while
analysis 2 uses $D_{\rm
North/South}=($2.4\text{\,}\mathrm{\SIUnitSymbolDegree}$,$3.7\text{\,}\mathrm{\SIUnitSymbolDegree}$)$
for UHECRs with arrival directions from the Galactic Northern and Southern
hemisphere, respectively. The different choices have been made based on the
respectively better sensitivity for the two analyses. Analysis 1 and 2
implement the deflection into their methods as described in the respective
sections 4.1 and 4.2. Analysis 3 employs a model-independent approach with
respect to the UHECR deflection.

##  4 Analysis methods

###  4.1 Unbinned neutrino point-source search with UHECR information

The goal of analysis 1 is to find point-like neutrino sources that are
spatially correlated with UHECR arrival directions within a region derived
from their magnetic deflection estimate. The search for neutrino sources
utilizes the unbinned maximum-likelihood analysis commonly used in IceCube
(Aartsen et al., 2017d, 2019, 2020b). This enables an easy combination of the
high-statistics, full-sky neutrino data sets, which are the PS data sets of
IceCube and ANTARES described in subsection 2.1 and subsection 2.2,
respectively. In addition to this standard method, the magnetic deflection
regions, defined by the UHECR arrival directions, energy, and scaling factor,
are used for constraining the possible source regions.

The unbinned neutrino likelihood in the source direction
$\vec{\Omega}=(\alpha,\delta)$ in right ascension and declination consists of
the sum of a signal PDF, $S$, and a background PDF, $B$

$\displaystyle\mathcal{L}\left(n_{s},\gamma\right)$
$\displaystyle=\prod_{j=1}^{7}\prod_{i=1}^{N_{j}}\left[f_{j}(\sin\delta,\gamma)\frac{n_{s}}{N_{j}}\cdot
S(\vec{\Omega}_{i},E_{i},\sigma_{i}|\gamma,\vec{\Omega})\right.$
$\displaystyle+\left.\left(1-f_{j}(\sin\delta,\gamma)\frac{n_{s}}{N_{j}}\right)\cdot
B(\delta_{i},E_{i})\right].$ (2)

The likelihood product with index $i$ runs over all neutrino events within
each data set with index $j$ (see subsection 2.1 and subsection 2.2). The
likelihood product with index $j$ runs over all data sets to yield the final
likelihood function. The likelihood combination of the seven data sets by
ANTARES and IceCube is thus handled in the same way as described in Aartsen et
al. (2017d). The signal and background PDFs, $S$ and $B$, are evaluated for
four observables per neutrino event, weighted with the number of signal events
over total number of events per data set, $n_{s}/N_{j}$. The observables are
the reconstructed right ascension and declination, summarized as
$\vec{\Omega}_{i}=(\alpha_{i},\,\delta_{i})$, reconstructed energy, $E_{i}$,
and angular error estimator, $\sigma_{i}$. The signal PDF consists of two
terms: the first term is a declination-dependent reconstructed energy
distribution, where the underlying neutrino flux is modeled as a power law

$\frac{{\rm d}\Phi^{\nu}}{{\rm d}E}=\frac{{\rm d}N^{\nu}}{{\rm d}E\,{\rm
d}A\,{\rm d}t\,{\rm
d}\Omega}=\Phi^{\nu}_{0}\cdot\left(\frac{E}{$1\text{\,}\mathrm{GeV}$}\right)^{-\gamma}.$
(3)

Here, $\Phi^{\nu}_{0}$ is the flux normalization, $1\text{\,}\mathrm{GeV}$ is
the corresponding pivot energy, and $\gamma$ 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, $\sigma_{i}$, given by the angular error estimator of
each neutrino candidate on an event-by-event basis. The background PDF is
determined as function of reconstructed energy, $E_{i}$, and declination,
$\delta_{i}$, from randomized experimental data. A full description of the
signal and background PDFs, $S$ and $B$, can be found in Aartsen et al.
(2017d). The proportionality factor, $f_{j}(\sin\delta,\gamma)$, is the
relative signal acceptance per neutrino data set calculated from the expected
number of signal events via

$\displaystyle f_{j}(\sin\delta,\gamma)$ $\displaystyle=\frac{n_{j}}{n_{\rm
tot}}$ $\displaystyle=\frac{T^{j}\int({\rm d}\Phi^{\nu}/{\rm d}E)A_{\rm
eff}^{j}(E,\sin\delta)~{}{\rm d}E}{\sum_{j=1}^{7}T^{j}\int({\rm
d}\Phi^{\nu}/{\rm d}E)A_{\rm eff}^{j}(E,\sin\delta)~{}{\rm d}E}.$ (4)

This factor is used to correctly weight the signal contribution of each data
set $j$ for a given source declination and spectral index. Per data set $j$,
the livetime is denoted with $T^{j}$ and the effective area is denoted with
$A_{\rm eff}^{j}$. Figure 5 shows the proportionality factor of each data set
as a function of declination for two spectral indices, $\gamma=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).

Figure 5: The relative signal contribution, $n_{j}/n_{\rm tot}$, for analysis
1 of the different configurations of IceCube and the ANTARES point-source data
set as a function of the declination $\delta$. The left plot shows a hard
underlying power-law signal spectrum (Equation 3) with index $\gamma=2.0$,
while the right plot shows a softer spectrum with $\gamma=2.5$. The relative
signal contributions from partial detector configurations of IceCube are shown
in blue, while the contributions of the full detector configuration of IceCube
with 86 strings are shown in orange. The ANTARES contribution is shown in
gray. The lines as listed in the rightmost legend indicate the contributions
from each data set individually.

The best-fit signal parameters, $\hat{n}_{s}$ and $\hat{\gamma}$, at a given
source position, $\vec{\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 livetime and effective area of each data
set via

$n_{s}=\sum_{j=1}^{7}T^{j}\int\frac{{\rm d}\Phi^{\nu}}{{\rm d}E}A_{\rm
eff}^{j}(E,\sin\delta)~{}{\rm d}E.$ (5)

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

$\mathrm{TS}_{\nu}(\vec{\Omega})=2\left.\log\left(\frac{\mathcal{L}\left(\hat{n}_{s},\hat{\gamma}\right)}{\mathcal{L}(n_{s}=0)}\right)\right|_{\vec{\Omega}}.$
(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
$0.2\text{\,}\mathrm{\SIUnitSymbolDegree}$. This results in a skymap of the TS
values at each grid center, $\vec{\Omega}_{\rm 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{\Omega}_{k}^{\rm CR}$, which possibly originated in the direction of
$\vec{\Omega}_{\rm grid}$, is defined as a 2D Gaussian. Its width is the
quadratic sum of the magnetic deflection estimate, $\sigma_{\rm MD}$ of
Equation 1, and the UHECR angular reconstruction error, $\sigma_{{\rm
CR}}^{2}$, which is $0.9\text{\,}\mathrm{\SIUnitSymbolDegree}$ for Auger
events and $1.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ for TA events

$\mathcal{P}_{k}\propto\exp\left(-\frac{(\vec{\Omega}_{\rm
grid}-\vec{\Omega}_{k}^{\rm CR})^{2}}{2(\sigma_{\rm MD}^{2}+\sigma_{{\rm
CR}}^{2})}\right).$ (7)

The term $\mathcal{P}_{k}$ acts as a spatial constraint, which is multiplied
to the neutrino likelihood function defined in Equation 2 via
$\mathcal{L}\rightarrow\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}\rightarrow\mathrm{TS}+2\cdot\log(\mathcal{P}_{k})$. Finally, the
maximum of the combined UHECR-neutrino TS is found at the best-fit grid point,
$\hat{\vec{\Omega}}_{\rm grid}$, via

$\displaystyle\mathrm{TS}^{\rm comb}_{k}(\hat{\vec{\Omega}}_{\rm grid})$
$\displaystyle=\max\left(\mathrm{TS}_{\nu}(\vec{\Omega}_{\rm
grid})-\frac{(\vec{\Omega}_{\rm grid}-\vec{\Omega}_{k}^{\rm
CR})^{2}}{\sigma_{\rm MD}^{2}+\sigma_{{\rm CR}}^{2}}\right)$ (8)

The maximum marks the neutrino-source candidate which is the most-likely
counterpart to the respective UHECR event used to calculate the spatial
constraint. The normalization in Equation 7 adds only a constant term to the
TS, which can be omitted when calculating significances and $p$-values based
on pseudo-experiments (see below). As a third, final step, this procedure is
repeated for all UHECRs and all obtained TS values are added up yielding the
final test statistic. This search strategy was first developed in the context
of this point-source correlation analysis and first outlined in sec. 4 of
Schumacher (2019). It has already been applied also to other neutrino
correlation searches with spatially extended source localization, namely the
correlation with ANITA events (Aartsen et al., 2020c) and with gravitational
wave events (Aartsen et al., 2020d). Note that the analysis described here is
improved with respect to the previous search (The IceCube Collaboration et
al., 2016, sec. 5), as it models the displacement of a point-like neutrino
source with respect to the UHECR arrival direction based on the assumed
magnetic deflection. Here, the position and flux of a point-like source are
fit in the vicinity of the UHECR direction, while in the previous search, a
spatially-extended neutrino emission around the UHECR direction was fit.

Six different signal hypotheses are tested, which are characterized by a lower
cut on the UHECR energy, $E_{\rm CR}>E_{\rm cut}$ with $E_{\rm
cut}\in\\{70,\,85,\,100\\}\,\mathrm{EeV}$, and the scaling factor of the
deflection estimate, $C\in\\{1,2\\}$. The lower energy cut is a threshold for
selecting only the highest-energy UHECRs with the lowest deflection for this
analysis. No analogous energy cut is applied to the neutrino data. The scaling
factor $C$ is a model parameter for scaling the baseline magnetic deflection,
and it is not derived from the UHECR data. The baseline magnetic deflection is
$D_{0}=$3\text{\,}\mathrm{\SIUnitSymbolDegree}$$ for all UHECRs, which is then
converted into the spatial constraint of an individual UHECR event using
Equation 1 and Equation 7. The corresponding sensitivity and
$3\sigma$-discovery potential are evaluated based on the normalization of the
neutrino flux per source, $\Phi^{\nu}_{0}$ (see Equation 3). The $3\sigma$
threshold is chosen, since the background expectation needs to be calculated
based on pseudo-experiments, as described in the next paragraph. A
$5\sigma$-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 $\Phi^{\nu}_{0}$ in case of a
null measurement, while the discovery potential is defined as the median of
$\Phi^{\nu}_{0}$ required for a rejection of the background hypothesis with
$3\sigma$ significance.

Both sensitivity and discovery potential are calculated based on data
challenges, for which signal-like and background-like pseudo-experiments (PEs)
are generated. Since the calculation is based on PEs, the constant term in
Equation 7 and Equation 8, i.e. the normalization of the constraint term, can
be omitted. For all signal and background PEs, the UHECR arrival directions
are kept fixed. The background PEs are obtained by randomizing the
experimental neutrino data in right ascension. The signal PEs reflect the six
different signal hypotheses; one signal PE is based on experimental neutrino
data randomized in right ascension, to which Monte-Carlo neutrino events,
representing a neutrino source, are added. The location of this neutrino
source is chosen randomly within the spatial deflection constraint of one
UHECR event. This way, we mimic a neutrino source that is displaced with
respect to the UHECR direction due to the UHECR deflection. In the baseline
model of the signal PEs, all UHECR events have such an artificial neutrino
source in their vicinity. For one PE, all neutrino sources have an $E^{-2}$
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_{\rm 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}}=\left\lceil{f_{\rm corr}N_{\rm
CR}}\right\rceil$, where $N_{\rm CR}$ is the number of UHECRs passing the
energy cut. We choose the correlation fraction from three discrete values,
$f_{\rm corr}\in\\{1,\,0.5,\,0.1\\}$, where $f_{\rm corr}=1$ represents the
baseline model. Note that the correlation fraction is only used to choose the
number of neutrino sources in the signal PEs, while in a real experimental
data, it is not known a priori which UHECRs have a correlated neutrino source
and which do not. For $f_{\rm corr}<1$, only a number of
$N_{\nu}^{\mathrm{src}}<N_{\rm 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\sigma$-discovery potential in terms of
the flux normalization per neutrino source for all signal models and
hypotheses are reported in Table 2, which are also shown in Figure 6. We see
that a smaller correlation fraction increases the flux normalization per
neutrino source for both sensitivity and $3\sigma$-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.

It is noticeable that for a small correlation fraction of $f_{\rm corr}=0.1$,
the signal model with a larger scaling factor of $C=2$ yields a better
$3\sigma$-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\sigma$-discovery
potential due to the small number of sources present222There are 22, 9 and 4
sources for $f_{\rm corr}=0.1$ and the energy cuts of $E_{\rm
cut}\in\\{70,\,85,\,100\\}\,\mathrm{EeV}$, respectively.. 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.

Figure 6: Sensitivity (left) and $3\sigma$-discovery potential (right) of the
neutrino flux per source as function of the UHECR energy cut in EeV, for the
three correlation fractions and for the two scaling factors we used.

###  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

$\displaystyle\log\,\mathcal{L}(n_{s})=$
$\displaystyle\sum_{i=1}^{N_{\text{Auger}}}\log{\bigg{(}}\frac{n_{s}}{N_{\text{CR}}}S^{i}_{\text{Auger}}+\frac{N_{\text{CR}}-n_{s}}{N_{\text{CR}}}B^{i}_{\text{Auger}}{\bigg{)}}$
$\displaystyle+\sum_{i=1}^{N_{\text{TA}}}\log{\bigg{(}}\frac{n_{s}}{N_{\text{CR}}}S^{i}_{\text{TA}}+\frac{N_{\text{CR}}-n_{s}}{N_{\text{CR}}}B^{i}_{\text{TA}}{\bigg{)}},$
(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_{\rm CR}=N_{\text{Auger}}+N_{\text{TA}}$. In contrast to the likelihood
function of analysis 1, the signal PDFs, $S^{i}_{\text{Auger}}$ and
$S^{i}_{\text{TA}}$, and the background PDFs, $B^{i}_{\text{Auger}}$ and
$B^{i}_{\text{TA}}$, are stacked for all high-energy neutrinos such that
$n_{s}$ is a global parameter. The background PDFs per UHECR detector,
$B^{i}_{\rm det}$, are the normalized expectations for an isotropic UHECR flux
as function of declination, which correspond to the normalized detector
exposures (see Figure 3). The signal PDF per UHECR detector is defined as

$\displaystyle S^{i}_{\rm det}(\vec{\Omega_{i}},E_{i})=R_{\rm
det}(\delta_{i})\cdot\sum_{j=1}^{N_{\nu}}S_{j}(\vec{\Omega_{i}},\sigma(E_{i}))$
$\displaystyle{\rm with}\ \sigma^{2}=\sigma_{\rm MD}^{2}+\sigma_{{\rm
CR}}^{2},$ (10)

as a function of the UHECR arrival direction, $\vec{\Omega_{i}}$, and energy,
$E_{i}$. The term $R_{\rm det}$ is the relative exposure of each detector as a
function of the declination 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_{\rm det}$. The sum runs over all neutrino events, $N_{\nu}$,
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 Northern or Southern Galactic Hemisphere (see Figure 4), the baseline
magnetic deflection is thus

$\displaystyle D=C\cdot\left(D_{\rm North},D_{\rm
South}\right)=C\cdot\left($2.4\text{\,}\mathrm{\SIUnitSymbolDegree}$,$3.7\text{\,}\mathrm{\SIUnitSymbolDegree}$\right)$
$\displaystyle\text{with}\ C\in\\{1,2,3\\}.$ (11)

The final value of the deflection, $\sigma_{\rm 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$

${\rm TS}=2\log\frac{\mathcal{L}(\hat{n}_{s})}{\mathcal{L}(n_{s}=0)}.$ (12)

The significance is then determined based on the assumption that the
background expectation is distributed according to a $\chi^{2}$ function with
one degree of freedom, which has been verified using background-like PEs. The
test statistic is calculated separately for all track-like and all cascade-
like neutrino events, and separately for the three different deflections. The
analysis approach is essentially the same as published in The IceCube
Collaboration et al. (2016), except for the updated magnetic deflection
values, which are split into the Galactic Northern and Southern Hemisphere as
described in section 3.

The sensitivity and the $3\sigma$-discovery potential in terms of $n_{s}$ are
obtained with data challenges based on PEs. Here, the neutrino arrival
directions are kept fixed per PE, while the UHECR arrival directions and
energies are generated resembling the signal and background hypothesis. For
the background PEs, all UHECR arrival directions are derived from an isotropic
flux and thus according to the exposure of the respective UHECR experiments.
The energies of the UHECR events are sampled from a power-law proportional to
$E^{-4.2}$ for the Auger events and to $E^{-4.5}$ for the TA events, same as
in The IceCube Collaboration et al. (2016). For the signal PEs, a number of
$n_{s}$ UHECR arrival directions are distributed randomly based on their
respective spatial signal terms in Equation 10,
$S_{j}(\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\sigma$-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\sigma$-discovery
potential using neutrino tracks requires much fewer UHECRs than when using
neutrino cascades, which is expected due to the large differences in angular
reconstruction uncertainty of the two event topologies.

###  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_{\rm obs}(\alpha)$, is compared to the
mean number of pairs expected from pure background, i.e. uncorrelated arrival
directions, $\langle n_{\rm bckg}(\alpha)\rangle$. The relative excess of
pairs is defined as

$\frac{n_{\rm obs}(\alpha)}{\langle n_{\rm bckg}(\alpha)\rangle}-1.$ (13)

As the separation angle is not known a-priori, angles between
$1\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$30\text{\,}\mathrm{\SIUnitSymbolDegree}$ in steps of
$1\text{\,}\mathrm{\SIUnitSymbolDegree}$ are tested. The angle with the
largest deviation from the background expectation is chosen after the scan.

The significance of the experimental result is calculated with respect to
background PEs. Similar to analysis 2, the background PEs are generated with a
fixed set of neutrino arrival directions, and uncorrelated UHECR arrival
directions are generated according to the exposure of the respective UHECR
experiments. The energies are sampled from the same spectra as described in
subsection 4.2. As a cross-check, additional background PEs are generated with
the fixed set of UHECR arrival directions, and neutrino arrival directions are
randomized in right ascension. The different types of background PEs
approximate an isotropic flux of UHECR and high-energy neutrinos,
respectively.

Note that this analysis does not include an assumption of the magnetic
deflection of UHECRs, which makes it a robust and model-independent approach.
The analysis approach is the same as published in The IceCube Collaboration et
al. (2016).

##  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 as described in subsection 4.1, assuming the six
different combinations of signal parameters, i.e. combinations of scaling
factor, $C$, and lower energy cut, $E_{\rm 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_{\rm
cut}=$85\text{\,}\mathrm{E}\mathrm{e}\mathrm{V}$$. 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_{\rm cut}=$100\text{\,}\mathrm{EeV}$)$, the limits are
set equal to the sensitivity in order to not over-estimate the limits based on
an under-fluctuation in data. In addition to the baseline correlation fraction
of 100%, the limits are calculated for two smaller correlation fractions of
50% and 10%. In these cases, the limits are relaxed by about
$40\text{\,}\mathrm{\char 37\relax}60\text{\,}\mathrm{\char 37\relax}$ for
$f_{\rm corr}=0.5$ compared to $f_{\rm corr}=1$, and by about a factor of 3 to
almost 5 for $f_{\rm corr}=0.1$ compared to $f_{\rm corr}=1$.

Analysis parameters |  |
---|---|---
$D_{0}\cdot C$ | $3{}^{\circ}$ | $3{}^{\circ}$ | $3{}^{\circ}$ | $6{}^{\circ}$ | $6{}^{\circ}$ | $6{}^{\circ}$
$E_{\rm cut}$ | $70\text{\,}\mathrm{EeV}$ | $85\text{\,}\mathrm{EeV}$ | $100\text{\,}\mathrm{EeV}$ | $70\text{\,}\mathrm{EeV}$ | $85\text{\,}\mathrm{EeV}$ | $100\text{\,}\mathrm{EeV}$
Number of neutrino sources, $\mathbf{N_{\nu}^{\mathrm{src}}}$ |  |
$f_{\rm corr}=0.1$ | 22 | 9 | 4 | 22 | 9 | 4
$f_{\rm corr}=0.5$ | 106 | 41 | 20 | 106 | 41 | 20
$f_{\rm corr}=1$ | 211 | 82 | 40 | 211 | 82 | 40
Sensitivity |  |
$f_{\rm corr}$=0.1 | 5.6 | 7.1 | 9.8 | 5.4 | 8.5 | 9.5
$f_{\rm corr}$=0.5 | 2.4 | 2.6 | 3.1 | 2.6 | 3.4 | 4.2
$f_{\rm corr}$=1 | 1.4 | 1.7 | 2.1 | 1.8 | 2.2 | 2.9
$\mathbf{3\sigma}$ disc. potential |  |
$f_{\rm corr}$=0.1 | 7.2 | 8.5 | 11.6 | 6.7 | 8.0 | 10.8
$f_{\rm corr}$=0.5 | 3.7 | 3.8 | 4.9 | 3.7 | 4.1 | 5.1
$f_{\rm corr}$=1 | 2.6 | 2.6 | 3.2 | 2.7 | 3.1 | 3.7
90% C.L. upper limit |  |
$f_{\rm corr}$=0.1 | 6.4 | 9.2 | 9.8 | 6.7 | 10.9 | 10.6
$f_{\rm corr}$=0.5 | 2.8 | 3.6 | 3.1 | 3.3 | 4.7 | 4.5
$f_{\rm corr}$=1 | 1.7 | 2.3 | 2.1 | 2.3 | 3.2 | 3.1
Pre-trial $\mathbf{p}$-value | 0.33 | 0.23 |  $>$0.5 | 0.19 | 0.097 | 0.43
Table 2: Sensitivity, discovery potential and 90% C.L. upper limit for the
different analysis parameters, which are the UHECR energy cut, $E_{\rm cut}$,
magnetic deflection, $D_{0}\cdot C$, and correlation fraction, $f_{\rm corr}$,
for the point-source correlation analysis. The values are given as
normalization of the neutrino flux per source in the unit of
$10^{-10}\,\mathrm{GeV^{-1}\,cm^{-2}\,s^{-1}}$. The neutrino flux per source
is modeled with a power-law of the form ${\rm d}\Phi/{\rm
d}E=\Phi_{0}\cdot\left(E/$1\text{\,}\mathrm{GeV}$\right)^{-2}$. The last row
states all six experimentally obtained pre-trial $p$-values with respect to
the null hypothesis of an isotropic neutrino flux.

###  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 Northern and Southern Galactic hemisphere of (7.2∘, 11.1∘). Based on
the observed TS values, 90% C.L. limits on the number of UHECR events
correlated to the high-energy neutrinos are calculated using the PEs of the
corresponding signal hypothesis for each of the six tests. Since all
$p$-values using the track-like neutrinos are larger than 0.5, the limits are
set equal to the sensitivity in order to not over-estimate the limits based on
an under-fluctuation in data. All $p$-values and the corresponding limits are
reported in Table 3 together with the sensitivity and discovery potential.

Analysis parameters |  |  |
---|---|---|---
$D_{N/S}\cdot C$ | (2.4∘, 3.7∘) | $\times 2$ | $\times 3$
Sensitivity |  |  |
tracks | 9.8 | 9.8 | 10.3
cascades | 57 | 61 | 64
$\mathbf{3\sigma}$ disc. potential |  |  |
tracks | 20 | 21 | 21
cascades | 113 | 129 | 135
90% C.L. UL |  |  |
tracks | 9.8 | 9.8 | 10.3
cascades | 55 | 76 | 101
$\mathbf{p}$-values |  |  |
tracks |  $>$0.5 |  $>$0.5 |  $>$0.5
cascades |  $>$0.5 | 0.38 | 0.26
Table 3: Sensitivity, $3\sigma$ discovery potential, and 90% C.L. upper limits
in terms of number of UHECRs ($n_{s}$, c.f. subsection 4.2) as well as the
pre-trial $p$-values for the UHECR stacking analysis with the samples of high-
energy tracks and cascades assuming an isotropic flux of UHECRs.

###  5.3 Two-point angular correlation of UHECRs and high-energy neutrinos

Event type | bg. hypothesis | separation |  p-value
---|---|---|---
tracks | isotropic neutrinos | 14∘ | 0.23
cascades | isotropic neutrinos | 16∘ | 0.15
tracks | isotropic UHECRs | 10∘ |  $>$0.5
cascades | isotropic UHECRs | 16∘ | 0.18
Table 4: Post-trial p-value and best-fit angular separation for the 2pt-
correlation analysis obtained with the neutrino data sets of high-energy
tracks and cascades, as stated in the first column. The $p$-values are
corrected for choosing the largest deviation out of all maximum angular
separations. The respective background hypotheses are stated in the second
column. Figure 7: Relative excess of pairs, $n_{\rm obs}/\langle n_{\rm
exp}\rangle-1$, as a function of the maximum angular separation of the
neutrino and UHECR pairs. The experimental result is shown as black dots,
while the blue color bands show the regions containing the 1, 2 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.

The 2pt-correlation analysis is applied to the sets of neutrino tracks and
cascades separately, as well as for the background hypothesis of an isotropic
UHECRs flux and an isotropic high-energy neutrino flux. The significance of
the result with respect to the background hypothesis is corrected for the scan
over the separation angles. This results in four $p$-values and four
respective best-fit angular separation, as listed in Table 4. The largest
deviation from an isotropic flux of high-energy neutrinos ($p$-value=15%) is
found using cascades and a maximum angular separation of
$16\text{\,}\mathrm{\SIUnitSymbolDegree}$. 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., 2018a, b), is located at a redshift of $z=0.34$
(Paiano et al., 2018), corresponding to a comoving distance of about
$1.3\text{\,}\mathrm{G}\mathrm{p}\mathrm{c}$. This is beyond the assumed
horizon for UHECR sources in the local universe of up to about
$200\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$. 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\text{\,}\mathrm{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\text{\,}\mathrm{TeV}$ for the combined high-energy data sets and as low as
$100\text{\,}\mathrm{GeV}$ for the PS data set. The non-observation of UHECR-
neutrino correlations thus extends to even lower energies than considered by
Palladino et al. (2020).

Even if there were local sources of both UHECRs and neutrinos, they could be
transient phenomena, such that the UHECRs emitted over a short period of time
arrive much later due to their deflection, whereas the neutrinos travel on a
straight path. It has been estimated by Davoudifar (2011) that the deflection
in the EGMF causes a typical time delay on the order of $10^{5}$ years
considering a source distance of $50\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$
and an EGMF strength of $2\text{\,}\mathrm{n}\mathrm{G}$. A delay of a couple
of decades as expected from propagation in the GMF is already sufficient to
de-correlate the observed UHECRs and neutrinos, as the livetime of the
neutrino and UHECR data sets ranges between 7 and 13 years.

Another source of uncertainty is the mass composition and thus the charge of
UHECRs at the highest energies. The measurements of the UHECR composition
above $\sim$$50\text{\,}\mathrm{EeV}$ are not yet conclusive. However, several
composition constraints at the highest energies (Aab et al., 2014a, b) and the
lack of significant, magnetically-induced structures in the UHECR arrival
directions (Aab et al., 2020) suggest that less than 10% of UHECRs above
$\sim$$50\text{\,}\mathrm{EeV}$ might be light elements. Only the light UHECRs
are expected to have their source in the vicinity of their arrival direction
as derived from the benchmark deflection model. The sources of heavier UHECRs,
e.g. iron nuclei, are spatially almost uncorrelated to the UHECR arrival
directions due to the 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
analysis 1 and 2 is a simplification with respect to the PT11 and JF12 GMF
models (Pshirkov et al., 2011; Jansson & Farrar, 2012), which themselves are
subject to large uncertainties. Especially in the region of the Galactic
Plane, we expect deflections that are larger than the assumed mean value of
around $3\text{\,}\mathrm{\SIUnitSymbolDegree}$. In addition, coherent
deflection of UHECRs in the GMF or the deflection of UHECRs in the IGMF are
not explicitly accounted for. Overall, a coherent shift of UHECRs depending on
their arrival direction or an overall significantly larger deflection in the
GMF and IGMF can dilute the spatial correlation of UHECRs and neutrinos.
Nevertheless, a scaling of the overall deflection is covered partially with
the scaling factor applied in analyses 1 and 2, since the overall strength of
the deflection and the UHECR charge are largely degenerate parameters.

From a theoretical perspective, the connection of UHECRs and high-energy
neutrinos is plausible based on their similar energy budget (Waxman & Bahcall,
1999; Bahcall & Waxman, 2001; Murase et al., 2013; Ahlers & Halzen, 2018).
However, this connection could not be proven with the current data and
analysis approaches. It is to note that neutrinos with energies in the TeV–PeV
range originate most likely from cosmic rays with energies below
$\sim$$100\text{\,}\mathrm{P}\mathrm{e}\mathrm{V}$ (Murase & Waxman, 2016),
which is below the energy threshold of
$>$50\text{\,}\mathrm{E}\mathrm{e}\mathrm{V}$$ 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
$\sim$$100\text{\,}\mathrm{P}\mathrm{e}\mathrm{V}$ are confined in a
calorimetric environment and subsequently produce TeV–PeV neutrinos, while a
fraction of the cosmic rays are accelerated to the highest energies and escape
the source before interacting (Murase & Waxman, 2016; Ahlers & Halzen, 2018).
Although no indication of such a scenario has been found in this analysis, the
first indication of such a connection could be the observation of the Seyfert
Galaxy NGC 1068 as a potential neutrino source candidate (Aartsen et al.,
2020b; Inoue et al., 2020) and UHECR source candidate (Aab et al., 2018). As
such sources are numerous also within the UHECR horizon, dedicated future
searches correlating UHECR, photon, and neutrino observations might be able to
set constraints on specific source candidates, instead of relying on UHECR and
neutrino data alone. NGC 1068 as a potential neutrino source candidate
(Aartsen et al., 2020b; Inoue et al., 2020) and UHECR source candidate (Aab et
al., 2018) could serve as blueprint for a catalog of potential common sources
to be constrained with dedicated analyses.

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

## Acknowledgements

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
gratitude 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,
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, & 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 Belgian Science
Policy under IUAP VII/37 (ULB). This work was partially supported by the
grants ofThe 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 collaboration’s 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 which 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.

ANTARES, IceCube Neutrino Observatory, Pierre Auger Observatory, Telescope
Array

## References

  * Aab et al. (2014a) Aab, A., Abreu, P., Aglietta, M., et al. 2014a, Phys. Rev. D, 90, 122005, doi: 10.1103/PhysRevD.90.122005
  * Aab et al. (2014b) —. 2014b, Phys. Rev. D, 90, 122006, doi: 10.1103/PhysRevD.90.122006
  * Aab et al. (2015) Aab, A., Abreu, P., Aglietta, M., et al. 2015, ApJ, 804, 15, doi: 10.1088/0004-637X/804/1/15
  * Aab et al. (2016) —. 2016, arXiv e-prints, arXiv:1604.03637.  https://arxiv.org/abs/1604.03637
  * Aab et al. (2017) —. 2017, Science, 357, 1266, doi: 10.1126/science.aan4338
  * Aab et al. (2018) Aab, A., Abreu, P., Aglietta, M., et al. 2018, The Astrophysical Journal, 853, L29, doi: 10.3847/2041-8213/aaa66d
  * Aab et al. (2020) Aab, A., Abreu, P., Aglietta, M., et al. 2020, Journal of Instrumentation, 15, P10021, doi: 10.1088/1748-0221/15/10/P10021
  * Aab et al. (2020) Aab, A., Abreu, P., Aglietta, M., et al. 2020, Journal of Cosmology and Astroparticle Physics, 2020, 017, doi: 10.1088/1475-7516/2020/06/017
  * Aartsen et al. (2014a) Aartsen, M. G., Abbasi, R., Ackermann, M., et al. 2014a, Journal of Instrumentation, 9, P03009, doi: 10.1088/1748-0221/9/03/P03009
  * Aartsen et al. (2014b) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2014b, Phys. Rev. Lett., 113, 101101, doi: 10.1103/PhysRevLett.113.101101
  * Aartsen et al. (2014c) —. 2014c, ApJ, 796, 109, doi: 10.1088/0004-637X/796/2/109
  * Aartsen et al. (2016) Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2016, ApJ, 833, 3, doi: 10.3847/0004-637X/833/1/3
  * Aartsen et al. (2017a) —. 2017a, ApJ, 835, 45, doi: 10.3847/1538-4357/835/1/45
  * Aartsen et al. (2017b) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017b, Journal of Instrumentation, 12, P03012, doi: 10.1088/1748-0221/12/03/P03012
  * Aartsen et al. (2017c) —. 2017c, Astroparticle Physics, 92, 30, doi: 10.1016/j.astropartphys.2017.05.002
  * Aartsen et al. (2017d) Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017d, ApJ, 835, 151, doi: 10.3847/1538-4357/835/2/151
  * Aartsen et al. (2018a) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018a, Science, 361, eaat1378, doi: 10.1126/science.aat1378
  * Aartsen et al. (2018b) —. 2018b, Science, 361, 147, doi: 10.1126/science.aat2890
  * Aartsen et al. (2018c) —. 2018c, Phys. Rev. D, 98, 062003, doi: 10.1103/PhysRevD.98.062003
  * Aartsen et al. (2019) —. 2019, European Physical Journal C, 79, 234, doi: 10.1140/epjc/s10052-019-6680-0
  * Aartsen et al. (2020a) —. 2020a, Phys. Rev. Lett., 125, 121104, doi: 10.1103/PhysRevLett.125.121104
  * Aartsen et al. (2020b) —. 2020b, Phys. Rev. Lett., 124, 051103, doi: 10.1103/PhysRevLett.124.051103
  * Aartsen et al. (2020c) —. 2020c, ApJ, 892, 53, doi: 10.3847/1538-4357/ab791d
  * Aartsen et al. (2020d) —. 2020d, ApJ, 898, L10, doi: 10.3847/2041-8213/ab9d24
  * Abbasi et al. (2013) Abbasi, R., Abdou, Y., Ackermann, M., et al. 2013, Nuclear Instruments and Methods in Physics Research A, 700, 188, doi: 10.1016/j.nima.2012.10.067
  * Abbasi et al. (2021) Abbasi, R., Ackermann, M., Adams, J., et al. 2021, Phys. Rev. D, 104, 022002, doi: 10.1103/PhysRevD.104.022002
  * Abbasi et al. (2014) Abbasi, R. U., Abe, M., Abu-Zayyad, T., et al. 2014, ApJ, 790, L21, doi: 10.1088/2041-8205/790/2/L21
  * Abbasi et al. (2018) —. 2018, ApJ, 867, L27, doi: 10.3847/2041-8213/aaebf9
  * Abreu et al. (2011) Abreu, P., Aglietta, M., Ahn, E. J., et al. 2011, arXiv e-prints, arXiv:1107.4809.  https://arxiv.org/abs/1107.4809
  * Abreu et al. (2021) Abreu, P., Aglietta, M., Albury, J. M., et al. 2021, in Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021), Vol. 395, 307, doi: 10.22323/1.395.0307
  * Abu-Zayyad et al. (2012) Abu-Zayyad, T., Aida, R., Allen, M., et al. 2012, Nuclear Instruments and Methods in Physics Research A, 689, 87, doi: 10.1016/j.nima.2012.05.079
  * Ageron et al. (2011) Ageron, M., Aguilar, J. A., Al Samarai, I., et al. 2011, Nuclear Instruments and Methods in Physics Research A, 656, 11, doi: 10.1016/j.nima.2011.06.103
  * Aguilar et al. (2007) Aguilar, J. A., Albert, A., Ameli, F., et al. 2007, Nuclear Instruments and Methods in Physics Research A, 570, 107, doi: 10.1016/j.nima.2006.09.098
  * Ahlers & Halzen (2018) Ahlers, M., & Halzen, F. 2018, Progress in Particle and Nuclear Physics, 102, 73, doi: 10.1016/j.ppnp.2018.05.001
  * Albert et al. (2017) Albert, A., André, M., Anghinolfi, M., et al. 2017, Phys. Rev. D, 96, 082001, doi: 10.1103/PhysRevD.96.082001
  * Albert et al. (2018) —. 2018, ApJ, 863, L30, doi: 10.3847/2041-8213/aad8c0
  * Albert et al. (2020) Albert, A., André, M., Anghinolfi, M., et al. 2020, The Astrophysical Journal, 892, 92, doi: 10.3847/1538-4357/ab7afb
  * Alves Batista et al. (2019) Alves Batista, R., Biteau, J., Bustamante, M., et al. 2019, Frontiers in Astronomy and Space Sciences, 6, 23, doi: 10.3389/fspas.2019.00023
  * Anchordoqui et al. (2008) Anchordoqui, L. A., Hooper, D., Sarkar, S., & Taylor, A. M. 2008, Astroparticle Physics, 29, 1, doi: 10.1016/j.astropartphys.2007.10.006
  * ANTARES Collaboration et al. (2002) ANTARES Collaboration, Amram, P., Anghinolfi, M., et al. 2002, Nuclear Instruments and Methods in Physics Research A, 484, 369, doi: 10.1016/S0168-9002(01)02026-5
  * Aublin et al. (2019) Aublin, J., Coleiro, A., Kouchner, A., et al. 2019, EPJ Web Conf., 210, 03003, doi: 10.1051/epjconf/201921003003
  * Bahcall & Waxman (2001) Bahcall, J., & Waxman, E. 2001, Phys. Rev. D, 64, 023002, doi: 10.1103/PhysRevD.64.023002
  * Barbano (2019) Barbano, A. M. 2019, in Proceedings of 36th International Cosmic Ray Conference — PoS(ICRC2019), Vol. 358, 842, doi: 10.22323/1.358.0842
  * Biteau et al. (2019) Biteau, J., Bister, T., Caccianiga, L., et al. 2019, in European Physical Journal Web of Conferences, Vol. 210, European Physical Journal Web of Conferences, 01005, doi: 10.1051/epjconf/201921001005
  * Bonifazi & The Pierre Auger Collaboration (2009) Bonifazi, C., & The Pierre Auger Collaboration. 2009, Nuclear Physics B Proceedings Supplements, 190, 20, doi: 10.1016/j.nuclphysbps.2009.03.063
  * Davoudifar (2011) Davoudifar, P. 2011, in 32nd International Cosmic Ray Conference, Vol. 2, 230, doi: 10.7529/ICRC2011/V02/1183
  * Dawson (2019) Dawson, B. 2019, in Proceedings of 36th International Cosmic Ray Conference — PoS(ICRC2019), Vol. 358, 231, doi: 10.22323/1.358.0231
  * Durrer & Neronov (2013) Durrer, R., & Neronov, A. 2013, Astron. Astrophys. Rev., 21, 62, doi: 10.1007/s00159-013-0062-7
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Haack & Wiebusch (2017) Haack, C., & Wiebusch, C. 2017, in Proceedings of 35th International Cosmic Ray Conference — PoS(ICRC2017), Vol. 301, 1005, doi: 10.22323/1.301.1005
  * Illuminati et al. (2019) Illuminati, G., Aublin, J., & Navas, S. 2019, in Proceedings of 36th International Cosmic Ray Conference — PoS(ICRC2019), Vol. 358, 920, doi: 10.22323/1.358.0920
  * Inoue et al. (2020) Inoue, Y., Khangulyan, D., & Doi, A. 2020, The Astrophysical Journal, 891, L33, doi: 10.3847/2041-8213/ab7661
  * Jansson & Farrar (2012) Jansson, R., & Farrar, G. R. 2012, The Astrophysical Journal, 757, 14, doi: 10.1088/0004-637X/757/1/14
  * Kawai et al. (2008) Kawai, H., et al. 2008, Nuclear Physics B - Proceedings Supplements, 175-176, 221, doi: 10.1016/j.nuclphysbps.2007.11.002
  * Kopper (2017) Kopper, C. 2017, in Proceedings of 35th International Cosmic Ray Conference — PoS(ICRC2017), Vol. 301, 981, doi: 10.22323/1.301.0981
  * Kronberg (1994) Kronberg, P. P. 1994, Rept. Prog. Phys., 57, 325, doi: 10.1088/0034-4885/57/4/001
  * Meszaros (2018) Meszaros, P. 2018, PhT, 71, 36
  * Murase (2015) Murase, K. 2015, AIP Conference Proceedings, 1666, 040006, doi: 10.1063/1.4915555
  * Murase et al. (2013) Murase, K., Ahlers, M., & Lacki, B. C. 2013, Physical Review D, 88, 121301, doi: 10.1103/PhysRevD.88.121301
  * Murase et al. (2014) Murase, K., Inoue, Y., & Dermer, C. D. 2014, Phys. Rev. D, 90, 023007, doi: 10.1103/PhysRevD.90.023007
  * Murase & Waxman (2016) Murase, K., & Waxman, E. 2016, Physical Review D, 94, 103006, doi: 10.1103/PhysRevD.94.103006
  * Nagano & Watson (2000) Nagano, M., & Watson, A. A. 2000, Reviews of Modern Physics, 72, 689, doi: 10.1103/RevModPhys.72.689
  * Paiano et al. (2018) Paiano, S., Falomo, R., Treves, A., & Scarpa, R. 2018, ApJ, 854, L32, doi: 10.3847/2041-8213/aaad5e
  * Palladino et al. (2020) Palladino, A., van Vliet, A., Winter, W., & Franckowiak, A. 2020, Monthly Notices of the Royal Astronomical Society, 494, 4255, doi: 10.1093/mnras/staa1003
  * Particle Data Group et al. (2020) Particle Data Group, Zyla, P. A., Barnett, R. M., et al. 2020, Progress of Theoretical and Experimental Physics, 2020, 083C01, doi: 10.1093/ptep/ptaa104
  * Pshirkov et al. (2011) Pshirkov, M. S., Tinyakov, P. G., Kronberg, P. P., & Newton-McGee, K. J. 2011, The Astrophysical Journal, 738, 192, doi: 10.1088/0004-637X/738/2/192
  * Pshirkov et al. (2016) Pshirkov, M. S., Tinyakov, P. G., & Urban, F. R. 2016, Phys. Rev. Lett., 116, 191302, doi: 10.1103/PhysRevLett.116.191302
  * Rodrigues et al. (2021) Rodrigues, X., Heinze, J., Palladino, A., van Vliet, A., & Winter, W. 2021, Phys. Rev. Lett., 126, 191101, doi: 10.1103/PhysRevLett.126.191101
  * Schumacher (2019) Schumacher, L. 2019, in European Physical Journal Web of Conferences, Vol. 207, European Physical Journal Web of Conferences, 02010, doi: 10.1051/epjconf/201920702010
  * Sommers (2001) Sommers, P. 2001, Astroparticle Physics, 14, 271, doi: 10.1016/S0927-6505(00)00130-4
  * Stettner (2019) Stettner, J. 2019, in Proceedings of 36th International Cosmic Ray Conference — PoS(ICRC2019), Vol. 358, 1017, doi: 10.22323/1.358.1017
  * The IceCube Collaboration (2013) The IceCube Collaboration. 2013, Science, 342, 1242856, doi: 10.1126/science.1242856
  * The IceCube Collaboration et al. (2016) The IceCube Collaboration, The Pierre Auger Collaboration, & The Telescope Array Collaboration. 2016, J. Cosmology Astropart. Phys, 2016, 037, doi: 10.1088/1475-7516/2016/01/037
  * The Pierre Auger Collaboration (2015) The Pierre Auger Collaboration. 2015, Nuclear Instruments and Methods in Physics Research A, 798, 172, doi: 10.1016/j.nima.2015.06.058
  * Tkachev et al. (2021) Tkachev, I., Fujii, T., Ivanov, D., et al. 2021, in Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021), Vol. 395, 392, doi: 10.22323/1.395.0392
  * Tokuno et al. (2012) Tokuno, H., Tameda, Y., Takeda, M., et al. 2012, Nuclear Instruments and Methods in Physics Research A, 676, 54, doi: 10.1016/j.nima.2012.02.044
  * Troitsky et al. (2017) Troitsky, S., Fukushima, M., Ikeda, D., et al. 2017, in Proceedings of 35th International Cosmic Ray Conference — PoS(ICRC2017), Vol. 301, 548, doi: 10.22323/1.301.0548
  * Wandkowsky (2018) Wandkowsky, N. 2018, in XXVIII International Conference on Neutrino Physics and Astrophysics, 501, doi: 10.5281/zenodo.1301088
  * Waxman & Bahcall (1999) Waxman, E., & Bahcall, J. N. 1999, Physical Review D, 59, 023002, doi: 10.1103/PhysRevD.59.023002