# Search for magnetically-induced signatures in the arrival directions of
ultra-high-energy cosmic rays measured at the Pierre Auger Observatory

A. Aab  P. Abreu  M. Aglietta  J.M. Albury  I. Allekotte  A. Almela  J.
Alvarez Castillo  J. Alvarez-Muñiz  R. Alves Batista  G.A. Anastasi  L.
Anchordoqui  B. Andrada  S. Andringa  C. Aramo  P.R. Araújo Ferreira  H.
Asorey  P. Assis  G. Avila  A.M. Badescu  A. Bakalova  A. Balaceanu  F.
Barbato  R.J. Barreira Luz  K.H. Becker  J.A. Bellido  C. Berat  M.E. Bertaina
X. Bertou  P.L. Biermann  T. Bister  J. Biteau  A. Blanco  J. Blazek  C. Bleve
M. Boháčová  D. Boncioli  C. Bonifazi  L. Bonneau Arbeletche  N. Borodai  A.M.
Botti  J. Brack  T. Bretz  F.L. Briechle  P. Buchholz  A. Bueno  S. Buitink
M. Buscemi  K.S. Caballero-Mora  L. Caccianiga  L. Calcagni  A. Cancio  F.
Canfora  I. Caracas  J.M. Carceller  R. Caruso  A. Castellina  F. Catalani  G.
Cataldi  L. Cazon  M. Cerda  J.A. Chinellato  K. Choi  J. Chudoba  L. Chytka
R.W. Clay  A.C. Cobos Cerutti  R. Colalillo  A. Coleman  M.R. Coluccia  R.
Conceição  A. Condorelli  G. Consolati  F. Contreras  F. Convenga  C.E.
Covault  S. Dasso  K. Daumiller  B.R. Dawson  J.A. Day  R.M. de Almeida  J. de
Jesús  S.J. de Jong  G. De Mauro  J.R.T. de Mello Neto  I. De Mitri  J. de
Oliveira  D. de Oliveira Franco  V. de Souza  J. Debatin  M. del Río  O.
Deligny  N. Dhital  A. Di Matteo  M.L. Díaz Castro  C. Dobrigkeit  J.C.
D’Olivo  Q. Dorosti  R.C. dos Anjos  M.T. Dova  J. Ebr  R. Engel  I. Epicoco
M. Erdmann  C.O. Escobar  A. Etchegoyen  H. Falcke  J. Farmer  G. Farrar  A.C.
Fauth  N. Fazzini  F. Feldbusch  F. Fenu  B. Fick  J.M. Figueira  A. Filipčič
T. Fodran  M.M. Freire  T. Fujii  A. Fuster  C. Galea  C. Galelli  B. García
A.L. Garcia Vegas  H. Gemmeke  F. Gesualdi  A. Gherghel-Lascu  P.L. Ghia  U.
Giaccari  M. Giammarchi  M. Giller  J. Glombitza  F. Gobbi  G. Golup  M. Gómez
Berisso  P.F. Gómez Vitale  J.P. Gongora  N. González  I. Goos  D. Góra  A.
Gorgi  M. Gottowik  T.D. Grubb  F. Guarino  G.P. Guedes  E. Guido  S. Hahn  R.
Halliday  M.R. Hampel  P. Hansen  D. Harari  V.M. Harvey  A. Haungs  T.
Hebbeker  D. Heck  G.C. Hill  C. Hojvat  J.R. Hörandel  P. Horvath  M.
Hrabovský  T. Huege  J. Hulsman  A. Insolia  P.G. Isar  J.A. Johnsen  J.
Jurysek  A. Kääpä  K.H. Kampert  B. Keilhauer  J. Kemp  H.O. Klages  M.
Kleifges  J. Kleinfeller  M. Köpke  G. Kukec Mezek  B.L. Lago  D. LaHurd  R.G.
Lang  M.A. Leigui de Oliveira  V. Lenok  A. Letessier-Selvon  I. Lhenry-Yvon
D. Lo Presti  L. Lopes  R. López  R. Lorek  Q. Luce  A. Lucero  A. Machado
Payeras  M. Malacari  G. Mancarella  D. Mandat  B.C. Manning  J. Manshanden
P. Mantsch  S. Marafico  A.G. Mariazzi  I.C. Mariş  G. Marsella  D. Martello
H. Martinez  O. Martínez Bravo  M. Mastrodicasa  H.J. Mathes  J. Matthews  G.
Matthiae  E. Mayotte  P.O. Mazur  G. Medina-Tanco  D. Melo  A. Menshikov
K.-D. Merenda  S. Michal  M.I. Micheletti  L. Miramonti  D. Mockler  S.
Mollerach  F. Montanet  C. Morello  M. Mostafá  A.L. Müller  M.A. Muller  K.
Mulrey  R. Mussa  M. Muzio  W.M. Namasaka  L. Nellen  M. Niculescu-Oglinzanu
M. Niechciol  D. Nitz  D. Nosek  V. Novotny  L. Nožka  A Nucita  L.A. Núñez
M. Palatka  J. Pallotta  M.P. Panetta  P. Papenbreer  G. Parente  A. Parra  M.
Pech  F. Pedreira  J. Pȩkala  R. Pelayo  J. Peña-Rodriguez  J. Perez Armand
M. Perlin  L. Perrone  C. Peters  S. Petrera  T. Pierog  M. Pimenta  V.
Pirronello  M. Platino  B. Pont  M. Pothast  P. Privitera  M. Prouza  A.
Puyleart  S. Querchfeld  J. Rautenberg  D. Ravignani  M. Reininghaus  J. Ridky
F. Riehn  M. Risse  P. Ristori  V. Rizi  W. Rodrigues de Carvalho  J.
Rodriguez Rojo  M.J. Roncoroni  M. Roth  E. Roulet  A.C. Rovero  P. Ruehl
S.J. Saffi  A. Saftoiu  F. Salamida  H. Salazar  G. Salina  J.D. Sanabria
Gomez  F. Sánchez  E.M. Santos  E. Santos  F. Sarazin  R. Sarmento  C.
Sarmiento-Cano  R. Sato  P. Savina  C. Schäfer  V. Scherini  H. Schieler  M.
Schimassek  M. Schimp  F. Schlüter  D. Schmidt  O. Scholten  P. Schovánek
F.G. Schröder  S. Schröder  S.J. Sciutto  M. Scornavacche  R.C. Shellard  G.
Sigl  G. Silli  O. Sima  R. Šmída  P. Sommers  J.F. Soriano  J. Souchard  R.
Squartini  M. Stadelmaier  D. Stanca  S. Stanič  J. Stasielak  P. Stassi  A.
Streich  M. Suárez-Durán  T. Sudholz  T. Suomijärvi  A.D. Supanitsky  J. Šupík
Z. Szadkowski  A. Taboada  A. Tapia  C. Timmermans  P. Tobiska  C.J. Todero
Peixoto  B. Tomé  G. Torralba Elipe  A. Travaini  P. Travnicek  C. Trimarelli
M. Trini  M. Tueros  R. Ulrich  M. Unger  M. Urban  L. Vaclavek  M. Vacula
J.F. Valdés Galicia  I. Valiño  L. Valore  A. van Vliet  E. Varela  B. Vargas
Cárdenas  A. Vásquez-Ramírez  D. Veberič  C. Ventura  I.D. Vergara Quispe  V.
Verzi  J. Vicha  L. Villaseñor  J. Vink  S. Vorobiov  H. Wahlberg  A.A. Watson
M. Weber  A. Weindl  L. Wiencke  H. Wilczyński  T. Winchen  M. Wirtz  D.
Wittkowski  B. Wundheiler  A. Yushkov  O. Zapparrata  E. Zas  D. Zavrtanik  M.
Zavrtanik  L. Zehrer  A. Zepeda  M. Ziolkowski  and F. Zuccarello

###### Abstract

We search for signals of magnetically-induced effects in the arrival
directions of ultra-high-energy cosmic rays detected at the Pierre Auger
Observatory. We apply two different methods. One is a search for sets of
events that show a correlation between their arrival direction and the inverse
of their energy, which would be expected if they come from the same point-like
source, they have the same electric charge and their deflection is relatively
small and coherent. We refer to these sets of events as “multiplets”. The
second method, called “thrust”, is a principal axis analysis aimed to detect
the elongated patterns in a region of interest. We study the sensitivity of
both methods using a benchmark simulation and we apply them to data in two
different searches. The first search is done assuming as source candidates a
list of nearby active galactic nuclei and starburst galaxies. The second is an
all-sky blind search. We report the results and we find no statistically
significant features. We discuss the compatibility of these results with the
indications on the mass composition inferred from data of the Pierre Auger
Observatory.

##  1 Introduction

The identification of the sources of ultra-high-energy cosmic rays (UHECRs) is
one of the main unsolved challenges in astrophysics. Since UHECRs are charged
particles, their path from their sources to Earth is modified by the
extragalactic and Galactic magnetic fields they traverse, most notably by the
latter. The knowledge of these intervening fields is still poor, despite the
considerable experimental effort in the area (see e.g. [1, 2, 3] and
references therein).

Another element that makes the magnetic deflections difficult to predict is
the uncertain composition of UHECRs. This is due to the fact that measurements
of the maximum of the shower development, which depends on the mass of the
primary particle, have low statistics at the highest energies and its
interpretation depends on the modelling of the hadronic interactions. The
results from the Pierre Auger Observatory [4, 5, 6] indicate that the
composition becomes heavier with increasing energy. However, measurements do
not rule out a light nuclei fraction at the highest energies, which might
originate in a few nearby sources, different from the average ones. In such a
case, the search for magnetically-induced signatures in the arrival directions
of UHECRs could help identify this type of sources. Several analysis
techniques have been designed to capture this kind of effect, including
indirect ways such as in [7]. In this work, we show the results of two
different methods applied to perform this search in a direct way, which are
referred to as “multiplets” and “thrust”.

Multiplets are defined as a set of events that show a correlation between
their arrival direction and the inverse of their energy, which is expected if
they come from the same source, they have the same electric charge and their
deflections are small and remain coherent. The observation of multiplets could
enable the accurate identification of the direction of the source and could
also provide a new means to probe the Galactic magnetic field by inferring the
value of its integrated component orthogonal to the trajectory of cosmic rays.

A search for multiplets was performed by the Pierre Auger Collaboration in
[8]. In that work, data up to 31st December 2010 were used, which amounted to
a total exposure of 25,800 km2 sr yr. The results obtained were not
statistically significant. The largest multiplet found above 20 EeV had 12
events and the probability that it would appear by chance in an isotropic
distribution of events was $6\%$. There were also two independent 10-plets and
the chance probability of having at least three multiplets with a multiplicity
equal to or larger than 10 was $20\%$. With the larger dataset used in this
work, the number of events added to these multiplets is not statistically
significant when comparing it to the number of events that would be added if
the arrival directions were isotropic, with a $p$-value of $60.5\%$.

In the thrust method, an observable is built from a principal axis analysis in
a localized region of the sky, measuring the elongation of a pattern with
respect to the center of the region of interest (ROI). The thrust method was
applied to data by the Pierre Auger Collaboration in [9], with events detected
up to 19th March 2013, amounting to a total exposure of 32,800 km2 sr yr. The
measured distributions of the thrust observables with the centroid
corresponding to the highest energy cosmic rays did not reveal any local
patterns in the arrival directions of UHECRs.

In this work, we update these analyses with more statistics. The data set
used, which amounts to a total exposure of 101,900 km2 sr yr, is described in
Section 2. This is an increase of a factor 4 with respect to [8] and a factor
3 with respect to [9]. The methods used are discussed in detail in Section 3.
Their expected sensitivity is shown in Section 5, using a benchmark simulation
described in Section 4. The methods are applied to a targeted search on a
selection of nearby active galactic nuclei (AGNs) and starburst galaxies
(SBGs), candidates to be sites of ultra-high-energy (UHE) acceleration, and
the results are presented in Section 6. Moreover, in Section 7, we apply the
methods to an all-sky blind search. Finally, in Section 8, we present our
conclusions.

##  2 The Pierre Auger Observatory and the data set

The Pierre Auger Observatory [10, 11], located in Malargüe, Argentina (35.2∘
S, 69.5∘ W, 1400 m a.s.l.), consists of a large surface detector (SD) and an
air-fluorescence detector (FD), which are used to observe extensive air
showers generated by UHECRs in a complementary way. The SD is composed of an
array of 1660 water-Cherenkov stations, arranged in an equilateral triangular
grid, spread over an area of 3000 km2. It measures, with a duty cycle of
nearly 100%, the energy flow in the shower carried by electrons, photons and
muons that reach ground level. The FD, comprising 27 telescopes at four sites,
overlooks the surface array and observes, with a duty cycle close to 15%, the
fluorescence light emitted by nitrogen molecules excited by the particles
produced as the air shower develops in the atmosphere.

In this study, we select data recorded with the surface detector between 1st
January 2004 and 31st August 2018 with zenith angles $\theta\leq 80^{\circ}$.
The energy thresholds we use are 20 EeV and 40 EeV, which yield data sets of
6568 and 1119 events, respectively. The “vertical” events ($\theta\leq
60^{\circ}$) are required to have at least four active stations surrounding
the one with the highest signal, while the “inclined” events
($\theta>60^{\circ}$) are required to have at least five active stations. For
the vertical events, the reconstructed core must be inside an equilateral or
isosceles triangle of active stations. This ensures an accurate reconstruction
of the shower geometry and energy.

The arrival directions of these cosmic rays are determined from the relative
arrival times of the shower front in the triggered stations. The angular
resolution, defined as the radius around the true cosmic-ray direction that
would contain 68% of the reconstructed shower directions, is better than 1∘
for the energies considered here [12, 13]. The energy estimate is obtained
from the signals recorded by the SD stations [14, 15] and is calibrated using
“hybrid” events (detected simultaneously by SD and FD) and taking advantage of
the quasi-calorimetric energy determination obtained with the FD [16]. The
statistical uncertainty in the energy determination is smaller than 9% for the
energies considered and the systematic uncertainty on the absolute energy
scale is 14% [16].

At energies above 4 EeV, the SD trigger is fully efficient, and the
determination of the exposure is purely geometrical [17]. For the time period
considered and for the applied energy and zenith angle selection, it amounts
to 101,900 km2 sr yr.

##  3 Methods

###  3.1 Multiplets

If the deflection is small, it is a good approximation [18] to consider that
the relation between the arrival direction of the cosmic ray $\vec{\Theta}$
and the actual source direction $\vec{\Theta}_{\text{s}}$,

$\vec{\Theta}=\vec{\Theta}_{\text{s}}+\frac{Ze}{E}\int^{L}_{0}{\rm
d}\vec{\ell}\times\vec{B}(\vec{\ell})\,.$ (3.1)

This equation can be simplified to yield the following linear relation with
the inverse of the energy $E$ of the cosmic ray:

$\vec{\Theta}\simeq\vec{\Theta}_{\text{s}}+\frac{\vec{D}(\vec{\Theta}_{\text{s}})}{E}\,.$
(3.2)

Here $\vec{D}$ is the integral along the line of sight of the perpendicular
component of the magnetic field $\vec{B}$ times the charge $Ze$ of the
particle, and is approximated to a constant for a fixed source direction. The
“deflection power” $D\equiv|\vec{D}|$ will be given in units $1^{\circ}\times
100$ EeV.

To determine if a set of cosmic rays forms a multiplet, we first use a
coordinate system ($x,y$) in the tangent plane to the celestial sphere
(centered on the average direction of the events). The system is then rotated
to a system ($u,w$) with an angle such that the covariance

$\text{Cov}(w,1/E)=\frac{1}{N}\sum_{i=1}^{N}(w_{i}-\langle
w\rangle)(1/E_{i}-\langle 1/E\rangle)$ (3.3)

between the cosmic ray coordinates $w_{i}$ and their inverse energy $1/E_{i}$
is zero and the covariance $\text{Cov}(u,1/E)$ is maximal. This angle is given
by

$\alpha=\arctan\left(\frac{\text{Cov}(y,1/E)}{\text{Cov}(x,1/E)}\right).$
(3.4)

The correlation between the coordinate $u$ and $1/E$ is measured through the
correlation coefficient

$C(u,1/E)=\frac{\text{Cov}(u,1/E)}{\sqrt{\text{Var}(u)\text{Var}(1/E)}}$ (3.5)

where the variances are given by $\text{Var}(x)=\left\langle(x-\langle
x\rangle)^{2}\right\rangle$.

The cuts applied to identify possible candidates among background chance
alignments are: i) a minimum value for the correlation coefficient,
$C(u,1/E)>C_{\text{min}}$, and ii) a maximum value for the spread in the
transverse direction of the multiplet, $\max(|w_{i}-\langle
w\rangle|)<W_{\text{max}}$. In [8], the two parameters chosen were
$C_{\text{min}}=0.9$ and $W_{\text{max}}=1.5^{\circ}$. They were adopted
following a search for a compromise between maximizing the signal from a true
source and minimizing the background arising from chance alignments. To do so,
numerical simulations were performed with sets of protons from randomly-
located extragalactic sources, which were propagated through a bisymmetric
magnetic field with even symmetry and a simple random deflection with root
mean square amplitude $1.5^{\circ}\times$(20 EeV/$E$) to simulate the
turbulent field. In this work, we tested variations of these parameters using
the more up-to-date simulation described in Section 4 and, as it is shown in
Section 5.1, these values for the cuts still perform well. A third cut used
here consists in searching for multiplets which extend up to $20^{\circ}$ in
the sky and with energies above 40 EeV, so that the linear correlation of the
deflection with $1/E$ is expected to be valid for “light” nuclei, such as
emitted by proton and helium sources. We also test the sensitivity of the
method making a search of events with energies between 20 EeV and 40 EeV,
where the linear approximation is still valid for proton sources.

###  3.2 Thrust ratio

Deviations from a clear correlation between the energy and the deflection
angle can arise from large turbulent field components or from a mixed nuclear
composition of the cosmic rays, as deflections scale also with the nuclear
charges $Z_{i}$ of the particles. In this case, the thrust ratio,
$T_{2}/T_{3}$, can still capture an elongation of the arising pattern in the
arrival directions. The thrust observables $T_{k}$ $(k=1,2,3)$ are constructed
by successively maximizing the values

$T_{k}=\underset{\hat{n}_{k}}{\text{max}}\left(\frac{\sum_{i}|\omega_{i}^{-1}\,\vec{p}_{i}\cdot\hat{n}_{k}|}{\sum_{i}|\omega_{i}^{-1}\,\vec{p}_{i}|}\right)$
(3.6)

with respect to the axes $\hat{n}_{k}$ starting with $k=1$. The sum iterates
over all cosmic rays in the chosen region of interest (ROI), $\vec{p}_{i}$ is
the arrival direction of particle $i$ weighted with its energy and
$\omega_{i}$ is the exposure towards this direction. The principal axes
$\hat{n}_{k}$ are perpendicular to each other
($\hat{n}_{1}\perp\hat{n}_{2}\perp\hat{n}_{3}$), thus, we obtain by
construction for the thrust observables the relation $T_{1}>T_{2}>T_{3}$.

While the thrust axis $\hat{n}_{1}$ points radially to the local barycenter of
the energy distribution within the ROI, the axes $\hat{n}_{2}$ and
$\hat{n}_{3}$ form an orthonormal coordinate system tangential to the local
spherical plane. In this system, a signal-like structure would be
characterized by a strong collimation of arrival directions around the
$\hat{n}_{2}$-axis, and therefore by a high thrust observable $T_{2}$. To be
less affected by the radial cosmic-ray distribution within the ROI, we choose
instead the thrust ratio $T_{2}/T_{3}$ as the relevant observable. Based on
the simulations performed in Section 4, we choose for the ROI radius a value
of $r_{\text{ROI}}=0.3$ rad (${\sim}17.2^{\circ}$). A more detailed
description of the procedure can be found in [9].

##  4 Target selection and benchmark simulation

In this Section, we describe the selection of astrophysical objects that have
been discussed as promising candidates for UHE acceleration and which are used
in a targeted search. We also give details of the simulation we used to
evaluate the sensitivities of the two methods presented in the previous
Section.

###  4.1 Target selection

In [19], we reported on a potentially interesting indication of anisotropy in
the arrival direction of UHECRs given by an excess flux towards nearby
galaxies. A maximum-likelihood ratio test was used with two free parameters,
namely the fraction of events correlating with astrophysical objects and an
angular scale characterizing the clustering of cosmic rays around
extragalactic sources. Sky models of cosmic-ray density were constructed using
two catalogs of source classes studied by _Fermi_ -LAT [20, 21]: AGNs and
SBGs. The relative luminosity of each source in UHECRs was assumed to be
proportional to the $\gamma$-ray luminosity for the former and to the radio
emission for the latter, for which the sample of detected sources in the
$\gamma$-ray band remains scarce. When comparing both models to isotropic
expectations, the most promising result was found with SBGs, with a
statistical significance of $4.0\sigma$. The result found with AGNs was also
interesting, with a $2.7\sigma$ deviation from isotropy. Results were updated
in [24] and the statistical significances increased to $4.5\sigma$ and
$3.1\sigma$ for the SBG and AGN models, respectively. The Telescope Array
Collaboration tested these results with their observed events, but they found
not to have enough statistics to confirm these results [25].

In this work, we include a search for magnetically-induced effects in the
arrival directions of events close to potential candidates for UHE
acceleration. For the multiplet search, we look for sets of events with a
reconstructed source position within $3^{\circ}$ of these candidates. We also
apply the thrust method centering the ROI in the position of the source
candidates. We call these two analyses a “targeted search”.

To define the astrophysical objects that are chosen as source candidates, we
used the catalogs from [20, 21, 22, 23]. We divided the sample between AGNs
and SBGs. Since we are interested in “light” UHECRs, we computed the flux at
Earth as if these candidates would accelerate helium nuclei, accounting for
the interactions from the sources to Earth. If with the methods used in this
work we were able to detect such events, we would also be able to detect
protons in the case they were also accelerated at these energies. In Figure 1,
we show the flux contribution of these sources relative to the expected one
from the brightest source in each sample, as a function of the distance to the
object, taking also into account the non-uniform exposure on the celestial
sphere of the Observatory. We simulated events with an energy above 40 EeV and
with an energy spectrum at the source being ${\rm d}N/{\rm d}E\propto E^{-2}$.
We assume that these sources, different from the average ones, would
contribute to a small fraction of the total measured flux. Hence, the spectrum
assumed for these simulations would be compatible to the measured one. We
present the results marginalized over the uncertainties on the measured energy
and on the distance to the candidate.111The distance estimates employed in
this work are based on [26] and [27], exploiting in particular the Tully-
Fisher correlation or the tip of the red-giant branch. The typical uncertainty
on high-quality distance estimators is estimated to 10% [29], which we adopt
as the minimum uncertainty in the marginalization procedure. The largest
variance in distance estimates over different methods is found for NGC 1068,
for which redshift-based estimates go up to 19 Mpc [28]. Adopting the latter
distance for NGC 1068 does not impact the selection of this source within the
sample of candidate for emitters of elongated patterns of UHECR. We selected
the candidates using a cut on the helium flux contribution of $1\%$ relative
to the object with the highest contribution (Cen A and NGC 4945 for the AGN
and SBG catalogs, respectively). We also kept NGC 4631 since its contribution,
$0.96\%$, is barely below the threshold. In Table 1, we list the ten source
candidates that remain after this cut, which include three AGNs and seven
SBGs.

Figure 1: Relative flux with respect to the brightest object in each catalog
(AGNs in blue empty circles and SBGs in red filled circles) as a function of
the distance to the source candidate, accounting for the directional exposure
of the Observatory and assuming that the accelerated and detected events are
helium nuclei with an energy greater than 40 EeV. The results are marginalized
over the uncertainties on the measured energy and on the distance to the
candidate. Table 1: Selection of astrophysical objects which could be
candidates for UHE acceleration. Assuming that these candidates would
accelerate helium nuclei, these objects are expected to contribute at least
${\sim}1\%$ with respect to the strongest source candidate in each catalog,
accounting for propagation effects and the exposure of the Observatory. AGN
source candidates
---
Target |  Gal. longitude [∘] |  Gal. latitude [∘] |  Distance [Mpc]
Cen A |  309.5 |  $\phantom{+}19.4$ |  $\phantom{0}3.7$ [26]
M 87 |  283.8 |  $\phantom{+}74.5$ |  $17$ [26]
Fornax A |  240.2 |  $-56.7$ |  $17$ [26]
SBG source candidates
Target |  Gal. longitude [∘] |  Gal. latitude [∘] |  Distance [Mpc]
NGC 253 |  $\phantom{0}97.4$ |  $-88.0$ |  $\phantom{0}3.6$ [26]
NGC 4945 |  305.3 |  $\phantom{+}13.3$ |  $\phantom{0}3.7$ [26]
Circinus |  311.3 |  $-3.8$ |  $\phantom{0}4.2$ [27]
M 83 |  314.6 |  $\phantom{+}32.0$ |  $\phantom{0}4.7$ [26]
NGC 4631 |  142.8 |  $\phantom{+}84.2$ |  $\phantom{0}7.4$ [26]
NGC 1808 |  241.2 |  $-35.9$ |  $\phantom{0}9.1$ [26]
NGC 1068 |  172.1 |  $-51.9$ |  $10$ [27]

###  4.2 Benchmark simulation

To evaluate the performance of the methods presented in Section 3, we tested
their sensitivity in a model where the UHECRs are accelerated in the source
candidates presented in Section 4.1 (listed in Table 1) and then are
propagated through a model of the intervening magnetic fields. For the
propagation through the Galactic magnetic field, we use a lensing technique
where particles of negative charge are backtracked from Earth to outside the
Galaxy [30].

Extragalactic and Galactic magnetic fields deflect the UHECRs from their
sources to Earth. The former fields are the least known, but assuming that
they have a strength smaller than ${\sim}10^{-9}$ G and a correlation length
smaller than $1$ Mpc [31, 3, 32, 33], the deflections of helium nuclei coming
from the considered candidates at a distance less than 20 Mpc are negligible
with respect to the deflections that occur during the propagation within our
Galaxy.

The Galactic magnetic field can be further divided into regular and turbulent
components. The regular part of the Galactic magnetic field is expected to
give a dominant contribution to the UHECR deflections. There are several
models of the regular Galactic field in the literature and there are large
uncertainties (see e.g. [34, 35, 36]) so that their predictions cannot yet be
used directly as an ingredient in the analysis. To test the sensitivities of
the methods in this work, we chose to use the regular model of Jansson and
Farrar [37], and two different scenarios for the turbulent component. One,
from [38], has a “striated” random field (having the same orientation as the
local coherent field but with a strength and sign which varies on a small
scale) and a random Kolmogorov field with coherence length of $60$ pc
(referred to as “GMF-A” model). The model is empirical and data-driven, and
the striated poloidal, toroidal and spiral components allow for effects such
as a Galactic wind and the compression of an isotropic turbulent field. A
second model for the random field is based on the analysis of the Planck team,
which argues that the synchrotron emission map of WMAP used by [38] to fix the
strength of the turbulent component of GMF-A is likely too large [39], and
hence rescales the random field of [38] correspondingly. Therefore, we also
simulate a weaker version of the turbulent field with the same coherence
length but strength scaled-down by a factor $3$, and without the striated
field; we call this the “GMF-B” model.

To evaluate the sensitivity of the methods, we test scenarios where one of the
source candidates, listed in Table 1, accelerates a certain number of cosmic
rays, $N_{\text{s}}$, with an energy spectrum following a simple power law
with $E^{-2}$. Two different scenarios are assumed for the chemical
composition of the accelerated particles. In Scenario 1, the source
accelerates helium nuclei with an energy between 40 EeV and 200 EeV. In
Scenario 2, we assume that helium nuclei are accelerated with energies between
40 EeV to 80 EeV, and that there is a proton fraction between 20 EeV to 40
EeV. The number of accelerated helium nuclei is chosen to be, arbitrarily in
this example, twice as high as the number of the proton events. These energy
ranges were chosen assuming that the acceleration at the source follows a
simple rigidity-dependent mechanism as expected from acceleration by
electromagnetic fields. In Table 2, we show a brief summary of the chosen
scenarios.

In both scenarios, we propagate the cosmic rays from the sources through GMF-A
and GMF-B models. Hence, there are four simulation setups to be analyzed. We
simulate a total number of events222These numbers are slightly smaller than
the actual number of events in data, but the difference is expected to have a
minor impact on the sensitivity results that are presented. $N_{\rm tot}=900$
above 40 EeV in Scenario 1, and $N_{\rm tot}=6000$ above 20 EeV in Scenario 2.
In each simulation, each source candidate is studied separately, there are
$N_{s}$ events originating in the chosen object, plus $N_{\rm
tot}-N_{\text{s}}$ isotropically distributed events, which follow the
geometrical exposure of the Observatory and the measured energy spectrum. The
arrival directions and energies of the $N_{\text{s}}$ events are simulated
accounting for the experimental resolutions. Thus, we apply an experimental
energy uncertainty333This is a worst-case scenario compared to the energy
uncertainty reported in [16], which was found to be smaller than 9% for the
energies considered, by making improvements in the calibration, atmospheric
treatment, reconstruction, and invisible energy. of $14\%$ and an angular
resolution of $1^{\circ}$.

Table 2: Summary of the details of the scenarios chosen to study the sensitivity of the methods. | Scenario 1 | Scenario 2
---|---|---
Minimum Energy | 40 EeV | 20 EeV
Maximum Energy | 200 EeV | 80 EeV
Composition | He |  p for $E<40\,\text{EeV}$, He for $E>40\,\text{EeV}$
GMF models | A and B | A and B

In Figure 2, we show two examples of a simulation in Scenario 1, with Cen A as
the source candidate, a number of injected signal cosmic rays of
$N_{\text{s}}=50$ and both Galactic magnetic field models GMF-A (left panel)
and GMF-B (right panel), in Galactic coordinates.

Figure 2: Examples of two simulations with $N_{\text{s}}=50$ events from Cen A
in Scenario 1, for the GMF-A (left) and the weaker version GMF-B (right), in
Galactic coordinates.

##  5 Sensitivity studies

###  5.1 Sensitivity of the multiplets method

In this Section, we show the results of the sensitivity studies with the
multiplet method, assuming the acceleration Scenario 1 (described in the
previous Section 4), where helium nuclei are accelerated with energies between
40 EeV and 200 EeV. We perform $1000$ simulations where a given $N_{\text{s}}$
source events are simulated from each source candidate.

For each simulation of $N_{\text{s}}$ events, the significance of the
reconstructed multiplet is calculated by computing the fraction of simulations
of isotropically distributed events (following the geometric exposure of the
Observatory), with the same total number of events $N_{\rm tot}$ and with the
same energy spectrum, in which a multiplet with the same or larger
multiplicity and passing the same cuts appears by chance.

When injecting $N_{\text{s}}$ events, the reconstructed multiplet does not
necessarily have a multiplicity equal to $N_{\text{s}}$, since there can be
simulated events that do not pass the required cuts, and there can be
isotropically distributed events which align by chance with simulated ones
from the chosen source. To account for this fact, in the sensitivity studies
we present the mean chance probabilities computed as:

$\left<P(N_{\text{s}})\right>=\sum_{n}f(n)P_{\rm ch}(n),$ (5.1)

where $f(n)$ is the fraction of $n$-plets obtained after having simulated
$N_{\text{s}}$ source events and $P_{\rm ch}(n)$ is the isotropic chance
probability for a multiplet with $n$ or more events.

Figure 3: Sensitivity results for the multiplet method in Scenario 1 ($E>40$
EeV). Mean chance probability $\left<P(N_{\text{s}})\right>$ for different
number of injected cosmic rays $N_{\text{s}}$ as a function of the chosen
source, for GMF-A (left panel) and a comparison between GMF-A and GMF-B (right
panel). The dashed line indicates the one-sided $3\sigma$ chance probability.

In the left panel of Figure 3, we present the mean chance probabilities for
different cases of $N_{\text{s}}$ source events (shown with different
markers), for the ten source candidates considered, in Scenario 1 and with the
source events propagated through the GMF-A model. Due to the non-uniform
exposure of the Observatory and to the dependence of the deflections and of
the magnification in the Galactic magnetic field on the location of the source
[40, 41], the number of events from the source required to obtained a mean
chance probability corresponding to a $3\sigma$ level depends highly on the
source candidate. For example, $N_{\text{s}}=9$ injected cosmic rays for the
sources M 87 and NGC 1068 are needed to reach a $3\sigma$ level, while
$N_{\text{s}}=18$ injected cosmic rays are needed for the source NGC 4945.
This corresponds to signal fractions of $1\%$ and $2\%$, respectively.

In the right panel of Figure 3, we show a comparison of the two Galactic
magnetic field models, GMF-A and GMF-B, for $N_{\text{s}}$ from 9 to 11. As
expected, the multiplet method performs better for almost all source
candidates with the GMF-B model, where the random deflections are weaker.

The values obtained for the number of events from the sources are reasonable.
For instance, if we use the results published in [19] for the AGN and SBG
models, the anisotropic fractions are $(7\pm 4)\%$ and $(10\pm 4)\%$,
respectively. In the former case, if the three AGNs listed in Table 1 are the
largest contributors, and if we assume that their relative flux contribution
is the one shown in Figure 1, their respective contribution would be $6.8\%$
for Cen A and ${\sim}0.1\%$ for M 87 and Fornax A. Translating these fractions
to number of source events, it would mean $N_{\text{s}}\approx 61$ for Cen A
and $N_{\text{s}}\approx 0.8$ and $N_{\text{s}}\approx 0.9$ for M 87 and
Fornax A, respectively. If we compare to the sensitivity studies done for
instance with GMF-A model (see the left panel of Figure 3), to obtain a mean
$p$-value corresponding to a $3\sigma$ level, one would need
$N_{\text{s}}\approx 16$ for Cen A, $N_{\text{s}}\approx 9$ for M 87 and
$N_{\text{s}}\approx 12$ for Fornax A. Though the results are not as favorable
for M 87 and Fornax A as they are for Cen A, their relative luminosity in
UHECRs might not be proportional to their $\gamma$-ray luminosity as we assume
here [42].

A similar comparison can be done with the SBG model. The number of events from
the sources needed for a $3\sigma$ detection are well within the expectations
from the results of [19] for the largest contributors, NGC 4945 and NGC 253
and not as much for NGC 1808 and NGC 4631, where less than one event is
expected and $N_{\text{s}}$ between 9 and 11 are needed (as it happens in the
case of M 87 and Fornax A).

We should point out that the angular scales of $7^{\circ}$ and $13^{\circ}$,
where the maximum deviation from isotropy was found in [19] for the AGN and
SBG models, could be compatible with a He composition given the uncertainties
on the GMF model. For example, using the simulations from Scenario 1, for all
sources, at least $95\%$ of the events would be within the angular window for
the SBG model. In the case of the AGN model, at least $40\%$ of the events
would be contained within such an angular window for the ten sources, with six
sources having a containment of at least $85\%$.

Using the simulation of Scenario 1 with the GMF-A model and a number of
injected source events of $N_{\text{s}}=10$, we tested variations of the cuts
on the correlation coefficient $C_{\text{min}}$ and on the spread in the
transverse direction of the multiplet $W_{\text{max}}$. The mean chance
probability $\left<P(N_{\text{s}})\right>$ as a function of the chosen source,
is shown in Figure 4. The original cuts of $C_{\text{min}}=0.9$ and
$W_{\text{max}}=1.5^{\circ}$ perform best for most sources. When making the
cuts looser, more of the injected signal events can be recovered in the
reconstructed multiplet. However, the significance of a larger multiplet
passing looser bounds is not necessarily better than the significance of the
smaller multiplet with stricter cuts, since in the former case, the
probability of chance alignments in the isotropically distributed events also
increases.

Figure 4: Mean chance probability $\left<P(N_{\text{s}})\right>$ for
$N_{\text{s}}=10$ as a function of the chosen source, for the GMF-A model and
in Scenario 1. The different markers show the results for different values of
the cuts $C_{\text{min}}$ and $W_{\text{max}}$. The dashed line indicates the
one-sided $3\sigma$ chance probability.

We studied the sensitivity of the method with Scenario 2. To perform the
search in this scenario, we first searched for multiplets with events with
energies above 40 EeV. After that, we selected the largest multiplet and we
searched for events between 20 EeV and 40 EeV which would correlate with the
events above 40 EeV. To pre-select these events, we chose those at an angular
distance to the event with the highest energy in the multiplet expected if the
deflection power was the one obtained with the higher energy events, with an
uncertainty of $20\%$, assuming that the events between 20 EeV and 40 EeV were
protons and those above 40 EeV were He. We once more applied a cut to the
maximum angular distance of $20^{\circ}$ so that the linear approximation
between $1/E$ and the arrival directions is still expected to be valid. We
simulated $N_{\text{s}}=10$ source events from each source candidate with
energies above 40 EeV and $N_{\text{s}}=5$ source events with energies between
20 EeV and 40 EeV and we compared the results to those obtained with Scenario
1 with the same number of source events $N_{\text{s}}=10$ above 40 EeV. We
present the results in Figure 5. Given the current statistics collected at the
Pierre Auger Observatory, the sensitivity of the method is better in Scenario
1 when comparing it to Scenario 2, for the same number of source events above
40 EeV. Hence, in this work, the method is only applied to events with
energies above 40 EeV.

Figure 5: Mean chance probability $\left<P(N_{\text{s}})\right>$ as a function
of the chosen source. The black empty squares and asterisks are the results
for Scenario 2, with $N_{\text{s}}=10$ above 40 EeV and $N_{s}=5$ between 20
EeV and 40 EeV, for the two GMF models considered, respectively. The blue
filled squares and empty circles are the results for Scenario 1, with
$N_{\text{s}}=10$ above 40 EeV, for the two GMF models considered,
respectively. The dashed line indicates the one-sided $3\sigma$ chance
probability.

###  5.2 Sensitivity of the thrust method

In the following the results of the sensitivity studies from the thrust ratio
are presented for both introduced scenarios from Section 4, with energy
thresholds of 20 EeV and 40 EeV. Again $1000$ simulations are performed where
$N_{\text{s}}$ events are assigned to the respective source candidate. For
each simulation, $100{,}000$ arrival directions were drawn isotropically in
every ROI for each of the $1000$ simulations, where the identical number of
cosmic rays follow the exposure in the ROI. By comparing the thrust ratio of
the ROI with the distribution obtained with the isotropic realizations, we
derive the isotropic chance probability $p$.

First, we test the sensitivity of the observable with respect to the only free
parameter, the radius of the region of interest $r_{\text{ROI}}$. In Figure 6,
we show the number of injected cosmic-ray source events $N_{\text{s}}$ which
leads to a $3\sigma$ significance to reject isotropy as a function of the
parameter $r_{\text{ROI}}$ for Scenario 1 in the left panel and Scenario 2 in
the right panel. Different colors and markers indicate different source
candidates from the catalog in Table 1. While there is an overall agreement of
the sensitivity curves between Scenarios 1 and 2, the optimal choice for the
parameter $r_{\text{ROI}}$ depends strongly on the considered source
candidate. This effect is due to the different deflection power of the
Galactic magnetic field in the different regions of the sky resulting in
patterns that are partially not covered within the ROI for a too small radius
$r_{\text{ROI}}$. To obtain an overall good performance, we choose an
intermediate radius of the ROI of $r_{\text{ROI}}=0.3$ rad in the following.

Figure 6: Number of injected cosmic ray source events $N_{\text{s}}$, which
leads to an one-sided $3\sigma$ chance probability as a function of the ROI
radius $r_{\text{ROI}}$. Applied on the benchmark simulation with the stronger
GMF-A model, for Scenario 1 ($E>40$ EeV) on the left panel and Scenario 2
($E>20$ EeV) on the right panel.

In the next step, we study the sensitivity of the thrust ratio as a function
of the number of injected signal cosmic rays $N_{\text{s}}$. In contrast to
the multiplet method, the patterns produced by the source candidates have a
larger impact on the sensitivity as can be seen by the large spread between
the lines in Figure 7. In case of Scenario 1, with the 40 EeV energy
threshold, a relatively high number of $N_{\text{s}}\approx 30$ to
$N_{\text{s}}\approx 60$ injected signal cosmic rays is needed to pass the
$3\sigma$ confidence level, which corresponds to signal fractions between
$3\%$ and $6\%$. This performance is inferior to that of the multiplet search,
where signal fractions of $1\%$ to $2\%$ are needed to reach the same
confidence level. This difference is explained by the fact that the thrust
ratio does not make a direct use of the energy information. The weaker
turbulent magnetic field model GMF-B just yields a minor improvement on the
performance: the same sensitivity as in the case of GMF-A is reached with two
fewer signal cosmic rays, on average.

Figure 7: Expected sensitivity for the thrust ratio observable on the
benchmark simulation with the stronger GMF-A model, for Scenario 1 ($E>40$
EeV) on the left panel and Scenario 2 ($E>20$ EeV) on the right panel. The red
line indicates the $3\sigma$ confidence level.

For Scenario 2, with the 20 EeV energy threshold and the GMF-A model, most of
the targets reach the $3\sigma$ confidence level between $N_{\text{s}}\approx
40$ and $N_{\text{s}}\approx 60$ signal cosmic rays, corresponding to signal
fractions between $0.7\%$ and $1.0\%$, as can be seen in the right panel of
Figure 7. The data above this energy threshold will only be analyzed with this
method in the following.

##  6 Results of the targeted search

In this Section, we present the results of the targeted search with the source
candidates listed in Table 1, applying the multiplet and thrust methods to the
data collected at the Pierre Auger Observatory described in Section 2.

The results are shown in Table 3 for each candidate, for the multiplet search
with an energy threshold of 40 EeV, and for the thrust ratio with energy
thresholds of 20 EeV and 40 EeV. All the probabilities found with respect to
an isotropic distribution of arrival directions are above $1\%$, thus no
significant signature of an alignment pattern around the chosen targets could
be found. These probabilities are not penalized for the trial of ten different
candidates.

Table 3: Isotropic chance probabilities for the targeted search with the
multiplet and thrust ratio methods applied to data collected at the Pierre
Auger Observatory. These probabilities are not penalized for the use of ten
different candidates. Isotropic chance probabilities
---
Target | Multiplets (40 EeV) | Thrust ratio (20 EeV) | Thrust ratio (40 EeV)
Cen A |  $1.2\times 10^{-2}$ |  $0.75$ |  $0.42$
M 87 |  $0.61$ |  $0.44$ |  $0.85$
Fornax A |  $0.96$ |  $0.21$ |  $1.9\times 10^{-2}$
NGC 253 |  $0.54$ |  $0.98$ |  $0.88$
NGC 4945 |  $0.25$ |  $2.9\times 10^{-2}$ |  $3.7\times 10^{-2}$
Circinus |  $0.99$ |  $0.82$ |  $0.58$
M 83 |  $0.20$ |  $0.14$ |  $0.54$
NGC 4631 |  — |  $0.59$ |  $0.85$
NGC 1808 |  $0.61$ |  $0.63$ |  $0.77$
NGC 1068 |  $0.75$ |  $6.0\times 10^{-2}$ |  $0.29$

The most interesting but non-significant result for the multiplet search, with
$E>40$ EeV, is with Cen A as source candidate, where an 8-plet is found with a
chance probability of $1.2\%$ with respect to isotropic expectations. The
fitted deflection power is $D=(9\pm 2)^{\circ}\times 100$ EeV. In the upper
left panel of Figure 8, we visualize the arrival directions of the events that
form the 8-plet.

In the case of the thrust ratio observable, the smallest $p$-values are found
with the source candidates Fornax A (with $E>40$ EeV), NGC 4945 (with both
energy thresholds) and NGC 1068 (with $E>20$ EeV), with chance probabilities
between $1.9\%$ and $6.0\%$. There is barely any overlap between the cosmic
rays within the multiplet of the Cen A region and the region of interest
evaluated for the thrust observable. In the upper right panel of Figure 8, we
present the arrival directions of the events in the ROI centered on Fornax A
above 40 EeV. The strong response of the thrust ratio ($T_{2}/T_{3}=1.686$)
seems to be an artifact of an asymmetric event distribution in the ROI rather
than a promising deflection signature. The lower panel of Figure 8 shows the
candidates NGC 4945 and NGC 1068 for energies above 20 EeV.

|
---|---
|
Figure 8: Arrival directions of events in Galactic coordinates from chosen sky
regions which yield the smallest isotropic chance probabilities in Table 3.
The upper row shows the AGN candidates Cen A with the 8-plet from the
multiplet method (left) and Fornax A with the thrust ratio (right), both for
the energy threshold of 40 EeV. The lower row shows the SBG candidates NGC
4945 (left) and NGC 1068 (right) evaluated with the thrust ratio on the energy
threshold of 20 EeV. The red bars show the thrust vectors
$T_{2}\cdot\hat{n}_{2}$ (solid) and $T_{3}\cdot\hat{n}_{3}$ (dashed) as
defined in Section 3.2.

It is instructive to look at the lack of magnetically-induced structure in the
light of the other measurements carried out at the Observatory. Assuming that
the GMF models used in our simulations are a good approximation of the actual
one, and considering the anisotropic fractions obtained in [19] for the AGN
and SBG models, we should have detected a significant signal in the multiplets
and thrust searches if all the events coming from the sources are light
elements (p or He). For instance, in the case of Cen A, for a total number of
1119 events, ${\sim}76$ events would have been expected from this source.
However, the putative multiplet found with Cen A as source candidate has a
multiplicity of 8 events only. Thus, under these assumptions, the lack of a
significant high-multiplicity signal appears to be consistent with a scenario
where less than ${\sim}10\%$ of cosmic rays above 40 EeV are light elements.
Such a fraction is in line with our composition constraints at the highest
energies [4, 5]. Overall, the various facets of our data are thus consistent.

##  7 Results of the all-sky blind search

We also evaluate both methods on the entire sky, for the multiplet method
above energies of $40$ EeV and for the thrust method for both energy
thresholds, $40$ EeV and $20$ EeV.

The largest multiplet found has a multiplicity of $10$ which may appear by
chance from an isotropic distribution of events with a probability of
$11.4\%$. The second largest multiplet has a multiplicity of $9$ events. The
chance probability of finding at least 2 multiplets with multiplicity larger
than or equal to 9 in isotropic simulations is $19.1\%$. The respective
deflection power values are $D=(8.0\pm 1.3)^{\circ}\times 100$ EeV and
$D=(12\pm 2)^{\circ}\times 100$ EeV. The arrival directions of the events that
form the 10-plet and 9-plet are shown in the left panel of Figure 9, as well
as the position of what would be their reconstructed source.

When searching for 8-plets, seven independent ones were found. We define
“independent” multiplets those sets of events that have less than half the
number of events in the set in common with others. The chance probability of
finding at least 9 multiplets with multiplicity larger than or equal to $8$ in
simulations of isotropic arrival directions is $5.6\%$. Thus, no statistically
significant results were found. In the right panel of Figure 9, we present the
arrival directions of the events in all the 8-plets as well as the
reconstructed positions of their potential sources.

Figure 9: Arrival directions in Galactic coordinates of the events that form
the 10-plet and 9-plet (left) and 8-plets (right) found on the all-sky blind
search with the multiplets method with an energy threshold of $E>40$ EeV. The
radius of the circles is proportional to the energy of the event. Plus signs
indicate the positions of the potential sources for each multiplet. The dashed
line represents the Equatorial plane.

For the thrust ratio observable, we used the HEALPix scheme [44] to locate the
ROIs equally on the sphere. For each ROI, isotropic realizations with the same
number of cosmic rays following the geometrical exposure within the ROI were
used to determine local $p$-values. Then, the lowest local $p$-value in the
sky is compared with the lowest local $p$-values of $10{,}000$ isotropic skies
to determine the post-trial significance. To overcome computational issues,
the isotropic distributions for the thrust ratio observable have been
parametrized and pre-calculated on a grid of equatorial declinations and
number of cosmic rays within the ROI (cf. Appendix). We select only ROIs with
declination lower than $+20^{\circ}$ and a cosmic ray number of at least $15$
to guarantee a sufficiently accurate parametrization.

The lowest but non-significant isotropic chance probability for the thrust
ratio observable was found for the energy threshold of 20 EeV with a ROI
centered at Galactic longitude $\ell=-9.8^{\circ}$ and Galactic latitude of
$b=-17.0^{\circ}$. The corresponding thrust-ratio was found to be
$T_{2}/T_{3}=1.362$ which translates to a post-trial $p$-value of $25\%$. For
the higher energy threshold of 40 EeV, we obtained a post-trial $p$-value of
$46.7\%$ with a thrust-ratio of $T_{2}/T_{3}=2.703$. With the ROI centered at
$\ell=-151.9^{\circ}$ and $b=-54.3^{\circ}$, we found the same artifact of an
unevenly filled ROI close to the Fornax A region (cf. Section 6).

Figure 10: Arrival directions of events in Galactic coordinates inside the ROI
located in the direction which yields the smallest local $p$-value in the all-
sky blind search for the energy thresholds of 40 EeV (left) and 20 EeV
(right).

##  8 Conclusions

Using a set of UHECRs collected at the Pierre Auger Observatory, with an
exposure of 101,900 km2 sr yr, we searched for indications of magnetically-
induced signatures in the arrival directions of the events. We applied two
different methods: a search for multiplets and the thrust method. We performed
both a search using a list of astrophysical sources, candidates to be sites of
UHE acceleration, and an all-sky blind search. We found no statistically
significant feature with any of the methods for the energy thresholds
considered, 40 EeV for the multiplets method and 20 EeV and 40 EeV for the
thrust one. The largest deviations from isotropic expectations were found with
a $p$-value of 1.9% around the Fornax A region for the thrust search above 40
EeV and with a $p$-value of 1.2% around the Cen A region for the multiplet
search above 40 EeV. Given the state of the art of the Galactic magnetic field
models, the lack of signal is along the lines of a scenario where less than
${\sim}10\%$ of cosmic rays above 40 EeV are light elements, which is
consistent with our composition constraints at the highest energies.

Future measurements will allow for disentangling of the electromagnetic and
muonic components on an event-by-event basis with the upgraded version of the
Pierre Auger Observatory [43]. Ultimately, the analyses presented in this
study will benefit from these measurements which will enable an accurate
analysis of mass-sensitive parameters for each event.

## Acknowledgments

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; 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 Science and Higher
Education, grant No. DIR/WK/2018/11; National Science Centre, Grants No.
2013/08/M/ST9/00322, No. 2016/23/B/ST9/01635 and No. HARMONIA
5–2013/10/M/ST9/00062, UMO-2016/22/M/ST9/00198; 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 – Romanian
Ministry of Education and Research, the Program Nucleu within MCI
(PN19150201/16N/2019 and PN19060102) and project PN-
III-P1-1.2-PCCDI-2017-0839/19PCCDI/2018 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 FPA2017-85197-P),
Xunta de Galicia (ED431C 2017/07), Junta de Andalucía (SOMM17/6104/UGR), 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.

## Appendix

The distribution of the thrust-ratio observable for an isotropic arrival
distribution and an energy spectrum as measured by the Pierre Auger
Observatory can be described by the density function

$f(t;a,b)=\frac{a^{3/b}\,b}{\Gamma(3/b)}\,(t-1)^{2}\,\exp(-a(t-1)^{b})\,\,,$
(8.1)

where $t=T_{2}/T_{3}$ is the thrust ratio, $\Gamma(x)$ is the Gamma function
and $a$ and $b$ are fit parameters. For both scenarios presented in Section 4,
the density function (8.1) has been fitted to $T_{2}/T_{3}$-distributions as
resulting from isotropically distributed cosmic rays. The parameters $a$ and
$b$ of this density function depend only on two variables, the equatorial
declination $\delta_{\text{ROI}}$ of the ROI center and the number of cosmic
rays within the ROI, $N_{\text{ROI}}$. We evaluated $a$ and $b$ of the density
function for $30$ different declinations $\delta_{\text{ROI}}$ ranging from
$-90^{\circ}$ to $+45^{\circ}$ and for each of them for $25$ different cosmic-
ray numbers $N_{\text{ROI}}$ ranging from $-4$ to $+4$ standard deviations in
the Poissonian distribution, which depend on the integrated exposure within
the ROI. For sufficiently high cosmic-ray numbers $N_{\text{ROI}}$ (above $15$
events), we find that the parameter $a$ is highly correlated to the parameter
$b$ as $a=3.59\,\exp(0.831\,b^{4})+2.73.$ For any continuous value of the
equatorial declination $\delta_{\text{ROI}}$ and the number of cosmic rays
$N_{\text{ROI}}$, we use a 2-dimensional interpolation algorithm to determine
the parameter $a$ and calculate the local $p$-value with the analytical
integral of Equation (8.1),

$p_{\text{val}}=\int_{t}^{\infty}f(x;a,b)\,\mathrm{d}x=\frac{\Gamma_{i}(3/b,\,a(t-1)^{b})}{\Gamma(3/b)}\,\,,$
(8.2)

where $\Gamma_{i}(x,\,y)$ is the incomplete upper Gamma function.

## References

  * [1] J.L. Han, _Observing Interstellar and Intergalactic Magnetic Fields_ , Annu. Rev. Astron. Astrophys. 55 (2017) 111-57.
  * [2] R. Beck, R. Wielebinski, _Magnetic Fields in Galaxies_ , Chapter 13 in “Planets, Stars and Stellar Systems”, Vol. 5: “Galactic Structure and Stellar Populations”, ed. G. Gilmore, Springer, Berlin 2013, ISBN 978-90-481-8817-8 [arXiv:1302.5663].
  * [3] R. Durrer, A. Neronov, _Cosmological Magnetic Fields: Their Generation, Evolution and Observation_ , Astron. Astrophys. Rev. 21 (2013) 62 [arXiv:1303.7121].
  * [4] A. Aab et al. [The Pierre Auger Collaboration], _Depth of Maximum of Air-Shower Profiles at the Auger Observatory: Measurements at Energies above $10^{17.8}$ eV_, Phys. Rev. D 90 (2014) 122005 [arXiv:1409.4809].
  * [5] A. Aab et al. [The Pierre Auger Collaboration], _Depths of Maximum of Air-Shower Profiles at the Pierre Auger Observatory: Composition Implications_ , Phys. Rev. D 90 (2014) 122006 [arXiv:1409.5083].
  * [6] A. Yushkov for the Pierre Auger Collaboration, _Mass composition of cosmic rays with energies above $10^{17.2}$ eV from the hybrid data of the Pierre Auger Observatory_, PoS (ICRC2019) 482 [arXiv:1909.09073].
  * [7] J.P. Lundquist and P.V. Sokolsky for the Telescope Array Collaboration, _Supergalactic Structure of Multiplets with the Telescope Array Surface Detector_ , EPJ Web Conf. 210 (2019) 01006
  * [8] P. Abreu et al. [The Pierre Auger Collaboration], _Search for signatures of magnetically-induced alignment in the arrival directions measured by the Pierre Auger Observatory_ , Astropart. Phys. 35 (2012) 354-361 [arXiv:1111.2472].
  * [9] A. Aab et al. [The Pierre Auger Collaboration], _Search for patterns by combining cosmic ray energy and arrival directions at the Pierre Auger Observatory_ , Eur. Phys. J. C 75 (2015), 269 [arXiv:1410.0515].
  * [10] J. Abraham et al. [The Pierre Auger Collaboration], _Properties and performance of the prototype instrument for the Pierre Auger Observatory_ , Nucl. Instrum. Meth. A 523 (2004) 50.
  * [11] A. Aab et al. [The Pierre Auger Collaboration], _The Pierre Auger Cosmic Ray Observatory_ , Nucl. Instrum. Meth. A 798 (2015) 172 [arXiv:1502.01323].
  * [12] C. Bonifazi for the Pierre Auger Collaboration, _The angular resolution of the Pierre Auger Observatory_ , Nucl. Phys. B Proc. Suppl. 190 (2009) 20 [arXiv:0901.3138].
  * [13] D. Mockler for the Pierre Auger Collaboration, _Reconstruction of Vertical Events Recorded by the Surface Detector of the Pierre Auger Observatory_ , PoS (ICRC2019) 353 [arXiv:1909.09073].
  * [14] M. Ave for the Pierre Auger Collaboration, _Reconstruction accuracy of the surface detector array of the Pierre Auger Observatory_ , Proc. of the $30^{th}$ Int. Cosmic Ray Conf. (2007) 4 307, [arXiv:0709.2125].
  * [15] A. Aab et al. [The Pierre Auger Collaboration], _Reconstruction of inclined air showers detected with the Pierre Auger Observatory_ , J. Cosmology & Astropart. Phys. 8 (2014) 19 [arXiv:1407.3214].
  * [16] B.R. Dawson for the Pierre Auger Collaboration, _The Energy Scale of the Pierre Auger Observatory_ , PoS (ICRC2019) 231 [arXiv:1909.09073].
  * [17] J. Abraham et al. [The Pierre Auger Collaboration], _Trigger and aperture of the surface detector array of the Pierre Auger Observatory_ , Nucl. Instrum. Meth. A 613 (2010) 29 [arXiv:1111.6764].
  * [18] G. Golup, D. Harari, S. Mollerach, E. Roulet, _Source position reconstruction and constraints on the galactic magnetic field from ultra-high energy cosmic rays_ , Astropart. Phys. 32 (2009) 269-277 [arXiv:0902.1742].
  * [19] A. Aab et al. [The Pierre Auger Collaboration], _An indication of anisotropy in arrival directions of ultra-high-energy cosmic rays through comparison to the flux pattern of extragalactic gamma-ray sources_ , Astrophys. J. Lett. 853 (2018), 2 [arXiv:1801.06160].
  * [20] M. Ackermann et al., _2FHL: the second catalog of hard _Fermi_ -LAT sources_, Astrophys. J. Suppl. Ser. 222 (2016) 5 [arXiv:1508.04449].
  * [21] M. Ackermann et al., _GeV observations of star-forming galaxies with the _Fermi_ Large Area Telescope_, Astrophys. J. 755 (2012) 164 [arXiv:1206.1346].
  * [22] M. Ajello et al., _3FHL: The Third Catalog of Hard _Fermi_ -LAT Sources_, Astrophys. J. 232 (2017) 2 [arXiv:1702.00664].
  * [23] J.K. Becker, P.L. Biermann, J. Dreyer, T.M. Kneiske, _Cosmic Rays VI-Starburst galaxies at multiwavelengths_ , [arXiv:0901.1775].
  * [24] L. Caccianiga for the Pierre Auger Collaboration, _Anisotropies of the highest energy cosmic-ray events recorded by the Pierre Auger Observatory in 15 years of operation_ , PoS (ICRC2019) 206 [arXiv:1909.09073].
  * [25] R.U. Abbasi et al. [The Telescope Array Collaboration], _Testing a Reported Correlation between Arrival Directions of Ultra-high-energy Cosmic Rays and a Flux Pattern from nearby Starburst Galaxies using Telescope Array Data_ , Astrophys. J. Lett. 867 (2018), 2.
  * [26] R.B. Tully et al., _COSMICFLOWS-2: the data_ , Astronom. J. 146 (2013) 4 [arXiv:1307.7213].
  * [27] R.B. Tully et al., _Our Peculiar Motion Away from the Local Void_ , Astrophys. J. 676 1 (2008) 184-205 [arXiv:0705.4139].
  * [28] M. Cappellari et al., _The ATLAS3D project - I. A volume-limited sample of 260 nearby early-type galaxies: science goals and selection criteria_ , Mon. Not. R. Astron. Soc. 413 2 (2011) 813-836 [arXiv:1012.1551].
  * [29] R.B. Tully, L. Rizzi, E.J. Shaya and H.M. Courtois, D.I. Makarov, B.A. Jacobs, _The Extragalactic Distance Database_ , arXiv:0902.3668.
  * [30] H. P. Bretz, M. Erdmann, P. Schiffer, D. Walz, T. Winchen, _PARSEC: A Parametrized Simulation Engine for Ultra-High Energy Cosmic Ray Protons_ , Astropart. Phys. C54 (2014) 110 [arXiv:1302.3761].
  * [31] P.P. Kronberg, _Extragalactic magnetic fields_ , Rept. Prog. Phys. 57 (1994) 325.
  * [32] M.S. Pshirkov, P.G. Tinyakov and F.R. Urban, _New Limits on Extragalactic Magnetic Fields from Rotation Measures_ , Phys. Rev. Lett. 116, 191302 (2016) [arXiv:1504.06546].
  * [33] K. Subramanian, _From Primordial Seed Magnetic Fields to the Galactic Dynamo_ , Galaxies 7 (2019) 47 [arXiv:1903.03744].
  * [34] G.R. Farrar, _The Galactic Magnetic Field and Ultrahigh-Energy Cosmic Ray Deflections_ , C. R. Physique XX (2014) [arXiv:1405.3680].
  * [35] M. Erdmann, G. Müller, M. Urban, M. Wirtz, _The Nuclear Window to the Extragalactic Universe_ , Astropart. Phys. 85 (2016) 54-64 [arXiv: 1607.01645].
  * [36] M. Unger, G.R. Farrar, _Progress in the Global Modeling of the Galactic Magnetic Field_ , EPJ Web of Conferences 210 (2019) 04005 [arXiv:1901.04720].
  * [37] R. Jansson, G.R. Farrar, _A New Model of the Galactic Magnetic Field_ , Astrophys. J. 757 (2012), 14 [arXiv:1204.3662].
  * [38] R. Jansson, G.R. Farrar, _The Galactic Magnetic Field_ , Astrophys. J. 761 (2012), L11 [arXiv:1210.7820].
  * [39] R. Adam [Planck Collaboration], _Planck intermediate results. XLII. Large-scale Galactic magnetic fields_ , Astron. Astrophys. 596 (2016) A103 [arXiv:1601.00546].
  * [40] D. Harari, S. Mollerach, E. Roulet, _The toes of the ultrahigh-energy cosmic ray spectrum_ , J. High Energy Phys. 9908 (1999) 022 [arXiv:astro-ph/9906309].
  * [41] G.R. Farrar, M.S. Sutherland, _Deflections of UHECRs in the Galactic magnetic field_ , J. Cosmol. Astropart. Phys. 1905 (2019) 05 004 [arXiv:1711.02730].
  * [42] J.H. Matthews, A.R. Bell, K.M. Blundell, A.T. Araudo, _Fornax A, Centaurus A, and other radio galaxies as sources of ultrahigh energy cosmic rays_ , Mon. Not. R. Astron. Soc. 479 (2018) 1 L76-L80 [arXiv:1805.01902].
  * [43] A. Castellina for the Pierre Auger Observatory, _AugerPrime: the Pierre Auger Observatory Upgrade_ , EPJ Web Conf. 210 (2019) 06002.
  * [44] K. Gorski, E. Hivon, A. Banday, et al., _HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere_ , Astrophys. J. 622 759 (2005) 771.