# Constraining models for the origin of ultra-high-energy cosmic rays with a
novel combined analysis of arrival directions, spectrum, and composition data
measured at the Pierre Auger Observatory

A. Abdul Halim  P. Abreu  M. Aglietta  I. Allekotte  K. Almeida Cheminant  A.
Almela  R. Aloisio  J. Alvarez-Muñiz  J. Ammerman Yebra  G.A. Anastasi  L.
Anchordoqui  B. Andrada  S. Andringa  C. Aramo  P.R. Araújo Ferreira  E.
Arnone  J. C. Arteaga Velázquez  H. Asorey  P. Assis  G. Avila  E. Avocone
A.M. Badescu  A. Bakalova  A. Balaceanu  F. Barbato  A. Bartz Mocellin  J.A.
Bellido  C. Berat  M.E. Bertaina  G. Bhatta  M. Bianciotto  P.L. Biermann  V.
Binet  K. Bismark  T. Bister  J. Biteau  J. Blazek  C. Bleve  J. Blümer  M.
Boháčová  D. Boncioli  C. Bonifazi  L. Bonneau Arbeletche  N. Borodai  J.
Brack  P.G. Brichetto Orchera  F.L. Briechle  A. Bueno  S. Buitink  M. Buscemi
M. Büsken  A. Bwembya  K.S. Caballero-Mora  L. Caccianiga  I. Caracas  R.
Caruso  A. Castellina  F. Catalani  G. Cataldi  L. Cazon  M. Cerda  J.A.
Chinellato  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  M. Conte
F. Convenga  D. Correia dos Santos  P.J. Costa  C.E. Covault  M. Cristinziani
C.S. Cruz Sanchez  S. Dasso  K. Daumiller  B.R. Dawson  R.M. de Almeida  J. de
Jesús  S.J. de Jong  J.R.T. de Mello Neto  I. De Mitri  J. de Oliveira  D. de
Oliveira Franco  F. de Palma  V. de Souza  E. De Vito  A. Del Popolo  O.
Deligny  L. Deval  A. di Matteo  M. Dobre  C. Dobrigkeit  J.C. D’Olivo  L.M.
Domingues Mendes  J.C. dos Anjos  R.C. dos Anjos  J. Ebr  F. Ellwanger  M.
Emam  R. Engel  I. Epicoco  M. Erdmann  A. Etchegoyen  C. Evoli  H. Falcke  J.
Farmer  G. Farrar  A.C. Fauth  N. Fazzini  F. Feldbusch  F. Fenu  A. Fernandes
B. Fick  J.M. Figueira  A. Filipčič  T. Fitoussi  B. Flaggs  T. Fodran  T.
Fujii  A. Fuster  C. Galea  C. Galelli  B. García  C. Gaudu  H. Gemmeke  F.
Gesualdi  A. Gherghel-Lascu  P.L. Ghia  U. Giaccari  M. Giammarchi  J.
Glombitza  F. Gobbi  F. Gollan  G. Golup  M. Gómez Berisso  P.F. Gómez Vitale
J.P. Gongora  J.M. González  N. González  I. Goos  D. Góra  A. Gorgi  M.
Gottowik  T.D. Grubb  F. Guarino  G.P. Guedes  E. Guido  S. Hahn  P. Hamal
M.R. Hampel  P. Hansen  D. Harari  V.M. Harvey  A. Haungs  T. Hebbeker  C.
Hojvat  J.R. Hörandel  P. Horvath  M. Hrabovský  T. Huege  A. Insolia  P.G.
Isar  P. Janecek  J.A. Johnsen  J. Jurysek  A. Kääpä  K.H. Kampert  B.
Keilhauer  A. Khakurdikar  V.V. Kizakke Covilakam  H.O. Klages  M. Kleifges
F. Knapp  N. Kunka  B.L. Lago  N. Langner  M.A. Leigui de Oliveira  Y Lema-
Capeans  V. Lenok  A. Letessier-Selvon  I. Lhenry-Yvon  D. Lo Presti  L. Lopes
L. Lu  Q. Luce  J.P. Lundquist  A. Machado Payeras  M. Majercakova  D. Mandat
B.C. Manning  P. Mantsch  S. Marafico  F.M. Mariani  A.G. Mariazzi  I.C. Mariş
G. Marsella  D. Martello  S. Martinelli  O. Martínez Bravo  M.A. Martins  M.
Mastrodicasa  H.J. Mathes  J. Matthews  G. Matthiae  E. Mayotte  S. Mayotte
P.O. Mazur  G. Medina-Tanco  J. Meinert  D. Melo  A. Menshikov  C. Merx  S.
Michal  M.I. Micheletti  L. Miramonti  S. Mollerach  F. Montanet  L. Morejon
C. Morello  A.L. Müller  K. Mulrey  R. Mussa  M. Muzio  W.M. Namasaka  A.
Nasr-Esfahani  L. Nellen  G. Nicora  M. Niculescu-Oglinzanu  M. Niechciol  D.
Nitz  D. Nosek  V. Novotny  L. Nožka  A Nucita  L.A. Núñez  C. Oliveira  M.
Palatka  J. Pallotta  G. Parente  J. Pawlowsky  M. Pech  J. Pȩkala  R. Pelayo
L.A.S. Pereira  E.E. Pereira Martins  J. Perez Armand  C. Pérez Bertolli  L.
Perrone  S. Petrera  C. Petrucci  T. Pierog  M. Pimenta  M. Platino  B. Pont
M. Pothast  M. Pourmohammad Shahvar  P. Privitera  M. Prouza  A. Puyleart  S.
Querchfeld  J. Rautenberg  D. Ravignani  M. Reininghaus  J. Ridky  F. Riehn
M. Risse  V. Rizi  W. Rodrigues de Carvalho  E. Rodriguez  J. Rodriguez Rojo
M.J. Roncoroni  S. Rossoni  M. Roth  A.C. Rovero  P. Ruehl  A. Saftoiu  M.
Saharan  F. Salamida  H. Salazar  G. Salina  J.D. Sanabria Gomez  F. Sánchez
E.M. Santos  E. Santos  F. Sarazin  R. Sarmento  R. Sato  P. Savina  C.M.
Schäfer  V. Scherini  H. Schieler  M. Schimassek  M. Schimp  F. Schlüter  D.
Schmidt  O. Scholten  H. Schoorlemmer  P. Schovánek  F.G. Schröder  J. Schulte
T. Schulz  S.J. Sciutto  M. Scornavacche  A. Segreto  S. Sehgal  S.U.
Shivashankara  G. Sigl  G. Silli  O. Sima  F. Simon  R. Smau  R. Šmída  P.
Sommers  J.F. Soriano  R. Squartini  M. Stadelmaier  D. Stanca  S. Stanič  J.
Stasielak  P. Stassi  M. Straub  A. Streich  M. Suárez-Durán  T. Suomijärvi
A.D. Supanitsky  Z. Svozilikova  Z. Szadkowski  A. Tapia  C. Taricco  C.
Timmermans  O. Tkachenko  P. Tobiska  C.J. Todero Peixoto  B. Tomé  Z. Torrès
A. Travaini  P. Travnicek  C. Trimarelli  M. Tueros  M. Unger  L. Vaclavek  M.
Vacula  J.F. Valdés Galicia  L. Valore  E. Varela  A. Vásquez-Ramírez  D.
Veberič  C. Ventura  I.D. Vergara Quispe  V. Verzi  J. Vicha  J. Vink  J.
Vlastimil  S. Vorobiov  C. Watanabe  A.A. Watson  A. Weindl  L. Wiencke  H.
Wilczyński  D. Wittkowski  B. Wundheiler  B. Yue  A. Yushkov  O. Zapparrata
E. Zas  D. Zavrtanik  and M. Zavrtanik

###### Abstract

The combined fit of the measured energy spectrum and shower maximum depth
distributions of ultra-high-energy cosmic rays is known to constrain the
parameters of astrophysical models with homogeneous source distributions.
Studies of the distribution of the cosmic-ray arrival directions show a better
agreement with models in which a fraction of the flux is non-isotropic and
associated with the nearby radio galaxy Centaurus A or with catalogs such as
that of starburst galaxies. Here, we present a novel combination of both
analyses by a simultaneous fit of arrival directions, energy spectrum, and
composition data measured at the Pierre Auger Observatory.

We find that a model containing a flux contribution from the starburst galaxy
catalog of around 20% at 40 EeV with a magnetic field blurring of around
$20^{\circ}$ for a rigidity of 10 EV provides a fair simultaneous description
of all three observables. The starburst galaxy model is favored with a
significance of $4.5\sigma$ (considering experimental systematic effects)
compared to a reference model with only homogeneously distributed background
sources. By investigating a scenario with Centaurus A as a single source in
combination with the homogeneous background, we confirm that this region of
the sky provides the dominant contribution to the observed anisotropy signal.
Models containing a catalog of jetted active galactic nuclei whose flux scales
with the $\gamma$-ray emission are, however, disfavored as they cannot
adequately describe the measured arrival directions.

##  1 Introduction

The origin and acceleration mechanisms of ultra-high-energy cosmic rays
(UHECRs) still remain a mystery today. Several challenges have to be faced in
the search for UHECR sources, among others the stochastic nature of the
interactions during propagation, the not directly measurable cosmic-ray
composition at the highest energies, and deflections by largely uncertain
cosmic magnetic fields. However, significant advances have been made both in
theoretical modeling and in the amount and quality of measured data in recent
years. With the Pierre Auger Observatory [1] an exposure of more than
$100{,}000$ km2 sr yr [2] has been accumulated, allowing us to determine the
characteristics of arriving UHECRs with unprecedented precision. Three
properties that could provide insights into the cosmic-ray sources, and that
are investigated in this work, are the cosmic-ray energy spectrum, the
distribution of shower maxima tracing the mass composition, and the arrival
directions.

For instance, the arrival directions have already revealed a dipolar
distribution above $8~{}\mathrm{EeV}$ with a significance exceeding $5\sigma$.
It is directed away from the Galactic center, which serves as proof of an
extragalactic origin of UHECRs [3, 4]. Recently, an indication of a
correlation of the UHECR arrival directions above $\sim 40~{}\mathrm{EeV}$
with the directions of powerful extragalactic source candidates was found [5,
6, 7]. The catalogs that have been tested by the Pierre Auger Collaboration
are two different selections of jetted- and non-jetted active galactic nuclei,
a selection of starburst galaxies, a broader catalog tracing the large-scale
mass distribution, and the nearby radio galaxy Centaurus A. The largest
statistical significance of $4.0\sigma$ (1-sided) is achieved for the
selection of starburst galaxies [2]. But, mainly because all catalogs contain
strong sources located in the same direction around Centaurus A, currently no
catalog can truly be favored above the others [2].

The combination of the arrival directions with the other two observables, the
energy spectrum and the distribution of shower maxima, may help to solve this
puzzle, as they contain additional information about the nature and
propagation of UHECRs. The energy spectrum and shower maximum depth
distributions have been used in previous analyses of the Pierre Auger
Collaboration [8, 9, 10] in order to constrain generic source models of
UHECRs. In those works, the model comprises only homogeneous distributions of
identical extragalactic sources, possibly adding a Galactic component when
lower energies are included, and the arrival directions are not considered as
an observable. Through the comparison of the modeled energy spectrum and
shower maximum depth distributions to the measured ones, parameters describing
the injected spectrum and composition at the sources have been determined. The
previous works will be used as a basis for this study.

Different kinds of combined analyses of all three observables have been
conducted before by other authors, for example Refs. [11, 12, 13, 14, 15]. Due
to the large parameter space involved for sophisticated astrophysical models
which also describe the arrival directions, these works focus on specific
aspects and choose to simplify other ingredients of the model or fit
procedure. Often, a two-step approach is chosen [11, 14], in which the model
parameters are adapted first, and then a comparison of the predicted arrival
directions to the measured ones is conducted. In other studies,
simplifications regarding the modeling of the source emission or propagation
have been made, in order to focus, for example, on the determination of
parameters of the extragalactic magnetic field or of the emission parameters
of individual sources [15, 13, 12].

In this work, all three observables are combined in a simultaneous fit of the
energy spectrum, shower maximum distributions, and arrival directions measured
at the Pierre Auger Observatory, taking into account propagation effects and
determining the parameters describing the source setup and emission.
Additionally, the uncertainties and correlations associated with these
parameters are also determined during the fit, using a Bayesian approach. A
preliminary proof of concept of the method using realistic simulations was
presented before in Refs. [16, 17].

The astrophysical model is built in a straightforward way, based on the
medium-scale anisotropies seen at energies well above the ankle and
correlating with extragalactic objects [2, 6]. For that reason, this work
focuses on the highest energies, taking into account arrival direction data
above 16 EeV, and energy spectrum and shower depth distributions above 10 EeV.
The astrophysical model, following the findings of [2, 6], not only contains
homogeneously distributed sources, but also a variable contribution from the
catalogs of either starburst galaxies or jetted active galactic nuclei that
correlate with the UHECR arrival directions. Here, the contribution of each
source to the total flux is modeled depending on its distance and the injected
spectrum and composition common to all sources. Additionally, a rigidity and
possibly distance-dependent magnetic field blurring is included. As a
comparison, a model containing only Centaurus A as a single source in addition
to the homogeneously distributed background sources, is also investigated.

This work is structured as follows: in section 2, the astrophysical model is
described, going from the emission at the sources in section 2.1 over the
propagation in section 2.2 to the calculation of the energy-dependent
contributions of the source candidates in section 2.3. Then, the modeling of
the three observables on Earth is described in section 2.4. In section 3, the
used inference methods, the likelihood function, and the measured data sets
are described. The fit results are discussed in section 4, and the influence
of the most important experimental systematic uncertainties is evaluated in
section 4.1.3. Finally, a discussion of the results and a conclusion is
provided in section 5.

##  2 Astrophysical models for the origin of UHECRs

The astrophysical model in this work is based on the one already successfully
employed before [8] in which only homogeneous background sources were used.
That setup, without a contribution of catalog sources, will also be
investigated in this work as a baseline and will hence be called the reference
model.

###  2.1 Source candidates and acceleration

In the following, the setup of the full astrophysical model will be explained
in detail, starting with the source distribution including catalog and
background sources.

##### Catalogs of starburst galaxies and active galactic nuclei:

The arrival directions of UHECRs measured with the Pierre Auger Observatory
exhibit intermediate-scale correlations with candidate sources from different
catalogs [6]. As described above, the astrophysical model in this work is
built up as a sum of homogeneous background sources and prominent foreground
catalog sources whose contributions are studied individually. As the costly
calculation of individual fluxes for all source candidates prohibits the use
of catalogs with a large number of sources, the shortest two updated ones from
the latest analysis [2, 18] are selected for this analysis. Catalogs with a
large number of source candidates, e.g. tracing the extragalactic matter
distribution, can rather be studied by calculating the contributions of
different distance shells [12] than that of individual sources (which is
however outside the scope of this work). The first catalog, which reaches the
highest significance of $4.0\sigma$ (1-sided) for a correlation with the UHECR
arrival directions, is a selection of 44 starburst galaxies (SBGs). The second
catalog is a selection of 26 radio-loud, jetted active galactic nuclei
selected by their $\gamma$-ray flux ($\gamma$-AGNs) with a significance of
$3.2\sigma$.

For the SBGs, the $\gamma$-ray flux measured with Fermi-LAT [19] is used as a
proxy for the expected UHECR flux. In practice, the flux weights are derived
from the radio flux in the 1.4 GHz band as that scales linearly with the
$\gamma$-ray emission [20]. Only SBGs with a flux above 0.3 Jy, located
between 1 Mpc and 130 Mpc, are selected. The selection is based on the catalog
from [21] with a few changes: the Large and Small Magellanic Clouds are
excluded and the Circinus galaxy is added as in [18]. The distances to the
catalog sources are derived by crossmatching the selection with the HyperLEDA
database111http://leda.univ-lyon1.fr/ [18]. The sources in the jetted
$\gamma$-AGN catalog are selected from the 3FHL catalog [19]. The selection
contains only sources with a $\gamma$-ray flux measured with Fermi-LAT
$>3.3\times 10^{-11}\,\text{cm}^{-2}\,\text{s}^{-1}$ between 10 GeV and 1 TeV,
which is also used as a UHECR flux proxy. This leads to a selection of 26
jetted $\gamma$-AGNs between 1 Mpc and 150 Mpc, where the distances are again
obtained from HyperLEDA.

The most important source candidates from the SBG catalog are NGC 4945, NGC
253, M83, and M82 (outside of the Observatory’s exposure), which are all
nearby at around $3-5$ Mpc. For the $\gamma$-AGN catalog, the dominating
source candidate is the powerful blazar Markarian 421 at a large distance of
$\sim 134$ Mpc. Also, the closest radio galaxy Centaurus A with a distance of
$\sim 3.5$ Mpc is part of the $\gamma$-AGN selection. Centaurus A has been
discussed as a UHECR source in previous publications [22, 23], and an
overdensity of events around its position has been observed and studied
throughout the years [24, 2, 5, 25]. For that reason, Centaurus A will also be
tested as a single source contributing on top of the homogeneous background.
Note that the direction and distance of the strongest source in the SBG
selection, NGC 4945, are almost identical to those of Centaurus A, so that a
differentiation of the two sources is not possible within the size of the
expected magnetic field blurring [2].

##### Background sources and source evolution:

In addition to the catalog sources, the model contains a distribution of
homogeneous background sources. These should account for unidentified, mostly
distant sources of the same kind as the catalog sources. A minimum distance of
3 Mpc is chosen for the background distribution, and a maximum distance
corresponding to redshift $z\approx 2.2$ (see above). The evolution of the
emissivity of the background sources as a function of redshift $z$ is modeled
as

$\psi(z)\propto(1+z)^{m}.$ (2.1)

For starforming regions such as SBGs the evolution follows the star formation
rate, which corresponds to $m=3.4$ up to $z=1$ [26, 27]. Under the assumption
that the flux of UHECRs from redshifts beyond $z\sim 1$ is negligibly small
[28, 8], we choose to simplify the evolution by keeping $m=3.4$ constant over
the whole range of redshifts of the background up to $z=2.2$. For AGNs, the
evolution depends on the luminosity of the considered class. These assumptions
for the source evolution are in line with other studies such as Refs. [10, 29,
30, 31].

##### Emission characteristics:

We assume that both catalog and background sources accelerate CRs according to
a rigidity-dependent process, so that the high-energy cutoff gives rise to a
Peters cycle [32] leading to an increase in the average mass with energy.
Analogously to Ref. [8], five representative stable elements are studied:
hydrogen (1H), helium (4He), nitrogen (14N), silicon (28Si), and iron (56Fe).
The generation rate at the sources $Q_{\mathrm{inj}}$ (defined as the number
of nuclei ejected per unit of energy, volume and time) as a function of the
injected energy $E_{\mathrm{inj}}$ and the mass number $A_{\mathrm{inj}}$ of
each representative element is modeled as a power-law with a broken
exponential cutoff function which sets in when the maximum energy
$E_{\mathrm{cut}}=Z_{A}R_{\mathrm{cut}}$ is reached

$Q_{\mathrm{inj}}(E_{\mathrm{inj}},A_{\mathrm{inj}})=Q_{0}\ a_{A}\
\Big{(}\frac{E_{\mathrm{inj}}}{10^{18}~{}\mathrm{eV}}\Big{)}^{-\gamma}\
~{}\begin{cases}1&;\ Z_{A}R_{\mathrm{cut}}>E_{\mathrm{inj}}\\\
\exp\big{(}1-\frac{E_{\mathrm{inj}}}{Z_{A}R_{\mathrm{cut}}}\big{)}&;\,Z_{A}R_{\mathrm{cut}}\leq
E_{\mathrm{inj}}\end{cases}.$ (2.2)

Here, $\gamma$ is the spectral index, $R_{\mathrm{cut}}$ is the maximum
rigidity of the source and $Q_{0}$ is a normalization for the UHECR generation
rate. $Z_{A}$ is the atomic number of the injected species with mass number
$A_{\mathrm{inj}}$ and $a_{A}$ is the fraction of particles of that species.
Due to the way eq. 2.2 is defined, the charge fractions $a_{A}$ denote the
fractions defined at energies below where the cutoff function sets in. It is
often more intuitive (and leads to a more stable fit) to utilize integral
fractions $I_{A}$ which denote the total emitted energy fractions of the
respective element above a given threshold. They can be calculated as222Note
that this definition denotes integral fractions calculated above a common
rigidity threshold and not a common energy threshold as used in [10]. Only for
hard spectral indices $\gamma<0$ do both definitions result in approximately
the same results.

$I_{A}=\frac{a_{A}\ Z_{A}^{2-\gamma}}{\sum_{A}a_{A}\ Z_{A}^{2-\gamma}}.$ (2.3)

It is important to note that the emission spectrum as described in eq. 2.2
does not necessarily represent the accelerated one since propagation and
interaction effects inside the source environment may alter it (see e.g. Refs.
[33, 34, 35, 36, 37, 38, 39]). So, if one determines the free parameters of
the source emission ($Q_{0}$, $\gamma$, $R_{\mathrm{cut}}$, $I_{\mathrm{H}}$,
$I_{\mathrm{He}}$, $I_{\mathrm{N}}$, $I_{\mathrm{Si}}$, $I_{\mathrm{Fe}}$) by
a fit to the observed data, these describe the state of the spectrum and
composition after leaving the source environment.

###  2.2 Propagation through intergalactic space

For the propagation of UHECRs from the sources to Earth, the software CRPropa3
[40] is used. With CRPropa3, in principle, propagation through 4-dimensional
time-space can be performed. But, because of the stochastic nature of
interactions and possible deflections in the extragalactic magnetic field
(EGMF), 4-dimensional simulations can take a long time, prohibiting the
determination of multiple free model parameters. The origin, strength, and
structure of the EGMF are largely unknown [41, 42] and models can differ
substantially [43]. Following the propagation theorem [44], a suppression of
the energy spectrum at lower energies is expected due to diffusion effects
[45, 46, 47], depending on the field strength and coherence length of the
EGMF. In this work, however, we only use the energy spectrum at the highest
energies $>10^{19}$ eV in the likelihood function (see section 3.1.2). For
that reason, we neglect the suppression effect333For an investigation of the
effect of a structured EGMF on the energy spectrum and shower depth
distributions $>10^{18.7}$ eV see [48]. Further studies of the EGMF effect are
planned for the future.. Regarding the arrival directions, the modeling of the
effect of possible deflections by cosmic magnetic fields will be described in
section 2.4.3.

Because we neglect the effect of the EGMF on the propagation, we can use
1-dimensional simulations instead of full 4-dimensional ones. Thus, a
propagation database of 1-dimensional CRPropa3 simulations is generated. The
database is then reweighted to account for the three-dimensional spatial setup
as will be described in section 2.3. Interactions with the CMB and EBL are
included in the form of photo-pion production, electron pair production, and
photodisintegration. To describe the photodisintegration, the TALYS model [49,
50] is used. For describing the EBL, the Gilmore model [51] is taken. The
influence of the interaction model on the results has been examined already in
e.g. Refs. [8, 9, 10]. Additionally, the propagation considers nuclear decay,
cosmological evolution of sources and background radiation fields, and
adiabatic energy losses of the cosmic rays during propagation. For the latter
two, the standard cosmological parameters of CRPropa3 assuming a flat universe
($H_{0}=67.3$ km/s/Mpc, $\Omega_{m}=0.315$, $\Omega_{\Lambda}=0.685$) [40] are
used.

Simulations are performed for multiple distances, injected energies, and mass
numbers, using $10^{4}$ injected particles for each bin. The injected energy
is binned in 150 logarithmic bins of width
$\log_{10}(E_{\mathrm{inj}}/\mathrm{eV})=0.02$ between $10^{18.0}$ eV and
$10^{21.0}$ eV. The light-travel distances $d$ are binned in 118 logarithmic
bins between 3 Mpc and 3342 Mpc (which approximately corresponds to redshift
$z\approx 2.2$). Simulations are done for the five representative elements
described in section 2.1, so that in total this amounts to $5\times 118\times
150=88{,}500$ simulations of $10^{4}$ injected particles each. Upon detection,
the arriving UHECRs are grouped into detected mass number bins
$A_{\mathrm{det}}\in(1,\ 2\text{--}4,\ 5\text{--}22,\ 23\text{--}38,\ $>39)
and detected energy bins
$E_{\mathrm{det}}\in(10^{18.0},10^{18.02}...,10^{20.98},10^{21.0})$ eV. In the
following, the propagation database will be denoted as
$P(E^{h}_{\mathrm{inj}},A^{i}_{\mathrm{inj}},E^{j}_{\mathrm{det}},A^{k}_{\mathrm{det}},d^{l})$,
where $h,i,j,k,l$ run over the bins described above. The detected energies are
rebinned to larger bins of width $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=0.1$
to match the measurements and to be in line with previous analyses [8, 10].
These larger energy bins will then be denoted as
$E^{e}_{\mathrm{det}}\in(10^{18.0},10^{18.1},...,10^{20.4})$ eV with the
running index $e$. The initial finer energy binning is necessary to include
systematic uncertainties on the energy spectrum, which will be discussed in
section 2.5.

###  2.3 Modeling the contributions from individual sources

The propagation database can now be used to map the injected flux from eq. 2.2
to Earth. For that, we calculate the histogram
$\mu(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},d^{l})$ describing the number
of arriving UHECRs in each energy $E^{e}_{\mathrm{det}}$ and mass bin
$A^{k}_{\mathrm{det}}$ from different distances $d^{l}$. It can be calculated
via matrix multiplication as

$\frac{\Delta\mu}{\Delta
d^{l}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},d^{l})=\frac{1}{4\pi}\sum_{h}\sum_{i}Q_{\mathrm{inj}}(E^{h}_{\mathrm{inj}},A^{i}_{\mathrm{inj}})\
P(E^{h}_{\mathrm{inj}},A^{i}_{\mathrm{inj}},E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},d^{l})\Delta
E_{\mathrm{inj}}^{h}.$ (2.4)

##### Calculation of the homogeneous background flux:

From this, the background arrival histogram $\mu_{\mathrm{back}}$ for the
homogeneous background sources can be calculated as

$\mu_{\mathrm{back}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})=\sum_{l}\frac{\Delta\mu}{\Delta
d^{l}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},d^{l})\
(1+z(d^{l}))^{m-1}\Delta d^{l}.$ (2.5)

The function $z(d^{l})$ describes the redshift as a function of the distance
of the respective bin. The factor $(1+z(d^{l}))^{-1}$ accounts for the
conversion from light-travel distance $d$ to comoving distance. The additional
factor $(1+z(d^{l}))^{m}$ takes care of the source evolution as described in
section 2.1.

##### Calculation of the flux from individual catalog sources:

For the catalog sources, the distances $d(C^{s}):=d^{s}$ and relative flux
weights $w_{\mathrm{flux}}(C^{s})$ of the respective candidates $C^{s}$ from
the catalogs described in section 2.1 are used. Since the propagation database
$P$ contains explicit distance bins $d^{l}$ and not specifically the source
distances $d^{s}$, a linear interpolation is taken into account to calculate
the fluxes for the distances $d^{s}$ (here denoted by the Kronecker
$\delta(d^{s}-d^{l})$):

$\mu_{\mathrm{sig}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})=\sum_{l}\frac{\Delta\mu}{\Delta
d^{l}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},d^{l})\ \delta(d^{s}-d^{l})\
w_{\mathrm{flux}}(C^{s})\Delta d^{l}$ (2.6)

The histogram
$\mu_{\mathrm{sig}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})$
describes the contribution of each source to the flux arriving at Earth for
each mass number bin in each detected energy bin. Overall, this contribution
depends on the flux weight $w_{\mathrm{flux}}(C^{s})$. The energy dependency
of the contribution, however, is determined by the source distance in
combination with propagation effects.

If the exposure on Earth in the direction of each source $C^{s}$ would be the
same, a simple summation over all candidate sources would lead to the detected
energies and masses for CRs from all catalog sources combined. For a ground-
based observatory, however, the differing exposure for the source candidates
has to be taken into account when computing the arrival histogram. Here, the
3-dimensional setup becomes important, even if one is only interested in
calculating energies and charges and not the arrival directions of the
incoming particles. The flux contribution of each catalog source due to the
non-uniform exposure of the Pierre Auger Observatory depends not only on the
source direction, but also on the modeling of the arrival directions and the
consideration of magnetic field deflections.

For a fast calculation of the exposure influence, a weight matrix
$w_{\mathrm{exp}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})$ is defined
and then multiplied with the arrival histogram $\mu_{\mathrm{sig}}$ to get the
exposure-weighted histogram $\hat{\mu}_{\mathrm{sig}}$. The weights
$w_{\mathrm{exp}}$ for each catalog source, mass number, and energy bin are
calculated considering the rigidity-dependent arrival directions modeling
which will be described below in section 2.4.3. The sum over all catalog
sources now gives the arrival histogram for the signal part:

$\hat{\mu}_{\mathrm{sig}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})=\sum_{s}\hat{\mu}_{\mathrm{sig}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})=\sum_{s}\mu_{\mathrm{sig}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})\
w_{\mathrm{exp}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s}).$ (2.7)

##### Definition of a signal fraction:

The two histograms for the catalog sources (signal) and the background have to
be combined to get the total arrival histogram. For that, another free model
parameter $f_{0}$ is introduced, the signal fraction that weights the two
parts. To have a signal fraction that is easily comparable to the arrival
directions correlation analysis [6], $f_{0}$ is defined as the signal fraction
in the detected energy bin $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.5-19.6$
($\approx 40$ EeV), corresponding to the running index $e=15$. Accordingly,
the two arrival histograms are weighted in the following way:

$\mu_{\mathrm{tot}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})=\underbrace{f_{0}\
\frac{\hat{\mu}_{\mathrm{sig}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})}{\sum_{k}\hat{\mu}_{\mathrm{sig}}(E^{15}_{\mathrm{det}},A^{k}_{\mathrm{det}})}}_{s_{\mathrm{sig}}}+\underbrace{(1-f_{0})\
\frac{\mu_{\mathrm{back}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})}{\sum_{k}\mu_{\mathrm{back}}(E^{15}_{\mathrm{det}},A^{k}_{\mathrm{det}})}}_{s_{\mathrm{back}}}:=\mu^{e,k}.$
(2.8)

The catalog contribution as a function of the detected energy bin can be
calculated from the two summands of eq. 2.8 as follows:

$f_{S}(f_{0},E^{e}_{\mathrm{det}})=\sum_{k}\frac{s_{\mathrm{sig}}}{s_{\mathrm{sig}}+s_{\mathrm{back}}}.$
(2.9)

As the catalog sources are on average much closer than the homogeneous
background, especially for a strong evolution,
$f_{S}(f_{0},E^{e}_{\mathrm{det}})$ usually rises with energy.

###  2.4 Simulated observables

The histogram $\mu_{\mathrm{tot}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})$
describes the arriving energies and masses on Earth according to the
astrophysical model. In order to find the best-fit free model parameters, the
observables predicted from
$\mu_{\mathrm{tot}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})$ have to be
compared to the measured ones on Earth. In our case, these observables are the
energy spectrum, shower maxima ($X_{\mathrm{max}}$) distributions, and arrival
directions as measured at the Pierre Auger Observatory.

In this section, it will be described how the observables are calculated from
the arrival histogram. For that, the detector effects that the measuring
devices induce on the observables have to be considered. A forward-folding
approach will be used to include the detector effects, meaning that the
detector response is part of the model so that the comparison between
simulations and data is performed at the detector level. Additionally, the
systematic uncertainties on the detector effects are taken care of, as will be
described in section 2.5.

####  2.4.1 Energy spectrum

The effects of the measurement with the surface detector (SD) of the Pierre
Auger Observatory on the event counts per energy bin can be expressed in terms
of an energy resolution, a bias, and an acceptance. Above $10^{18.4}$ eV, the
SD is fully efficient, so the acceptance is 100% in the energy range used for
this analysis. We take the response matrix $R^{e,e^{\prime}}$ from [52] which
transforms the energy spectrum predicted by the model to the reconstructed one
including resolution and bias effects. The spectrum $J(E^{e}_{\mathrm{det}})$
is then calculated as

$J(E^{e}_{\mathrm{det}})=\frac{\sum_{e^{\prime}}R^{e,e^{\prime}}\
\sum_{k}\mu^{e^{\prime},k}}{\mathcal{E}_{\mathrm{\theta\leq 60^{\circ}}}\
\Delta E^{e}_{\mathrm{det}}}$ (2.10)

with the total vertical (zenith angle $\theta\leq 60^{\circ}$) exposure
$\mathcal{E}_{\mathrm{\theta\leq 60^{\circ}}}$, the arrival histogram
$\mu^{e,k}$ from eq. 2.8, and the width of the energy bins $\Delta
E^{e}_{\mathrm{det}}$. Here, the vertical exposure is taken as we use only
vertical events for the energy spectrum data set while for the arrival
directions inclined events are also included (see section 3.2).

####  2.4.2 Depth of the shower maximum

The Gumbel distributions
$g(\tilde{X}_{\mathrm{max}}^{x}|E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})$
[53] are used to model the expected $X_{\mathrm{max}}$ distribution for each
energy and mass bin as in Ref. [8]. Here, the parameters are taken from the
updated generalized version of the Gumbel distributions, as in Ref. [10] for
the EPOS-LHC [54] hadronic interaction model. The $X_{\mathrm{max}}$ bins
denoted by $\tilde{X}_{\mathrm{max}}^{x}$ are given as
$\tilde{X}_{\mathrm{max}}^{x}\in(550,570,...,1050)\,\text{g/cm}^{-2}$.

The shower maximum depth distributions are influenced by the acceptance
$\mathcal{A}$ and the resolution $\mathcal{R}$ of the fluorescence detector
(FD) of the Pierre Auger Observatory, with which the longitudinal shower
development is observed. The parameterizations of the acceptance and
resolution can be found in Ref. [55]. For both effects, we use the updated
parameters as in Ref. [10]. The influence of the detector effects is folded
into the Gumbel distributions in the following way:

$~{}g(X_{\mathrm{max}}^{x}|E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})=\Big{(}g(\tilde{X}_{\mathrm{max}}^{x}|E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})\
\mathcal{A}(\tilde{X}_{\mathrm{max}}^{x}|E^{e}_{\mathrm{det}})\Big{)}\otimes\mathcal{R}(X_{\mathrm{max}}^{x}|\tilde{X}_{\mathrm{max}}^{x},E^{e}_{\mathrm{det}}).$
(2.11)

For the comparison of the measured $X_{\mathrm{max}}$ with the predicted model
values, a histogram is produced:

$~{}\mu(E^{e}_{\mathrm{det}},X_{\mathrm{max}}^{x})=c\sum_{k}\mu(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})\
g(\tilde{X}_{\mathrm{max}}|E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}}):=\mu^{e,x}.$
(2.12)

Here, $c$ normalizes the right side of the equation over all
$X_{\mathrm{max}}^{x}$ bins to one.

####  2.4.3 Arrival directions

In analogy with the energy spectrum and the shower maximum depth
distributions, the arrival directions are included as an observable to compare
them to the measured data. For that, the arrival directions need to be binned
in order to calculate the probability of cosmic rays arriving in different
directions for different energies. For the binning, the healpy package [56] is
used which is based on HEALPix444http://healpix.sf.net [57]. Considering
$\texttt{nside}=64$ it divides the sky into $n_{\mathrm{pix}}=49{,}152$ pixels
of the same angular size. The binning approximately corresponds to the
detector resolution of $~{}0.9^{\circ}$ [25] and can hence be considered to
account for this effect.

The modeling of the arrival directions is done separately for the background
and the catalog sources. For the background distribution, the arrival
directions are expected to be isotropic, so they are modeled following the SD
exposure. For that, the exposure $\mathcal{E}_{\mathrm{\theta\leq
80^{\circ}}}$ up to zenith angle $\theta=80^{\circ}$ in each healpy pixel
$\mathrm{pix}^{p}$ is calculated and a normalized isotropic arrival map is
defined which contains the exposure values of the $49{,}152$ pixels indexed
with the running index $p$:

$B(\mathrm{pix}^{p})=\mathcal{E}_{\mathrm{\theta\leq
80^{\circ}}}(\mathrm{pix}^{p})/\sum_{p}\mathcal{E}_{\mathrm{\theta\leq
80^{\circ}}}(\mathrm{pix}^{p}).$ (2.13)

For the modeling of the arrival directions of the catalog sources, the effect
of magnetic fields has to be considered. As described above, we neglect the
flux suppression induced by diffusion effects in the EGMF. Similarly, the
effect of the EGMF on the arrival directions is expected to be small [58] as
we are only interested in the highest energies in this work. The effect can be
parameterized by a beam-widening, based on the assumption that UHECRs at the
highest energies travel mostly in the non-resonant scattering regime [46].

Regarding the Galactic magnetic field (GMF), several sophisticated models
exist, e.g. Refs. [59, 60, 61, 62]. These do, however, differ substantially
regarding the predictions of the arrival directions [63, 64]. Most catalog
sources are not located in the direction of the Galactic disk, but rather the
Galactic halo, whose coherent component is even less known. Additionally, the
effect of the magnetic field of the Local Sheet on the UHECR arrival
directions is fairly unknown. As this work aims at investigating scenarios in
which source candidates correlating with the UHECR arrival directions are
actually the sources of UHECRs, we at this point refrain from including strong
coherent deflections which could break the source association555For
alternative studies, taking into account the possibility that the observed
overdensities around the catalog sources could actually be generated by
different sources whose flux is deflected into the candidate source directions
by a stronger coherent Galactic magnetic field model, compare to Refs. [12,
15].

Hence, similar to the arrival-directions correlation analysis which originally
revealed the overdensities around the catalog sources [6], the arrival
directions are described using a circular blurring around the source direction
(which could be caused by an extragalactic or a mostly turbulent Galactic
magnetic field) which is inserted into a von Mises–Fisher distribution [65]
$F(x|\mu,\kappa)$, the equivalent of a 2-dimensional Gaussian on the sphere.
The concentration parameter $\kappa$ describes the width of the distribution.
It is linked with the blurring angle $\delta$ via $\kappa=1/\delta^{2}$. In
Ref. [6] this blurring angle is a fixed value for all sources. In this study,
however, we include the rigidity dependency that is expected for magnetic
field deflections. This only becomes possible in our astrophysical model
because the model contains the whole simulation from sources to Earth, and
hence a prediction for the rigidities of the arriving particles for each
source candidate.

We parameterize the blurring angle
$\delta(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})$ as

$\delta(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})=\frac{\delta_{0}}{R_{\mathrm{det}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})/10\
\mathrm{EV}}$ (2.14)

with the free model parameter $\delta_{0}$ corresponding to the spread of
protons with energies of 10 EeV (or that of nitrogen with 70 EeV). Here, the
detected rigidity $R_{\mathrm{det}}=E_{\mathrm{det}}/(e\
Z_{\mathrm{det}}(A_{\mathrm{det}}))$ is used, and one should keep in mind that
when nuclei photodisintegrate during propagation the rigidity of the leading
fragment is essentially unchanged. It is determined from $E_{\mathrm{det}}$
and $A_{\mathrm{det}}$ by using the representative atomic numbers
$Z_{\text{det}}=1$, $2$, $7$, $14$ and $26$ for the detected mass number bins
$A_{\text{det}}\in\\{1\\}$, $\\{2,3,4\\}$, $\\{5,...,22\\}$, $\\{23,...,38\\}$
and $\\{39,...,56\\}$, respectively.

The von Mises–Fisher distributions are then stored in a histogram
$F(\mathrm{pix}^{p}|E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})$ using
again the same healpy pixels $\mathrm{pix}^{p}$ as for the background. The
final arrival probability map for the catalog sources is then calculated by
matrix multiplication:

$S(E^{e}_{\mathrm{det}},\mathrm{pix}^{p})=c^{e}\sum_{k,s}\mu_{\mathrm{sig}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})\
F(\mathrm{pix}^{p}|E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})\
\mathcal{E}_{\mathrm{\theta\leq 80^{\circ}}}(\mathrm{pix}^{p})$ (2.15)

with $c^{e}$ normalizing $\sum_{p}S(E^{e}_{\mathrm{det}},\mathrm{pix}^{p})=1$
in each energy bin as expected for a probability density. The multiplication
in eq. 2.15 weights each von Mises–Fisher distribution with the probability of
the respective element in the respective energy bin and also takes into
account the exposure weights of the different sources. The same von
Mises–Fisher distributions have also been used before to calculate the
exposure weighted flux (note that for the spectrum only zenith angles
$\theta\leq 60^{\circ}$ are used, see section 3.2):

$w_{\mathrm{exp}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})=\sum_{p}F(\mathrm{pix}^{p}|E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})\
\mathcal{E}_{\mathrm{\theta\leq 60^{\circ}}}(\mathrm{pix}^{p}).$ (2.16)

Finally, the probability densities for background and catalog sources are
added using the energy-dependent signal fraction function
$f_{S}(f_{0},E^{e}_{\mathrm{det}})$ from eq. 2.9:

$~{}\mathrm{pdf}(E^{e}_{\mathrm{det}},\mathrm{pix}^{p})=f_{S}(f_{0},E^{e}_{\mathrm{det}})\
S(E^{e}_{\mathrm{det}},\mathrm{pix}^{p})+(1-f_{S}(f_{0},E^{e}_{\mathrm{det}}))\
B(\mathrm{pix}^{p}).$ (2.17)

In [6], the signal fraction and the magnetic field blurring have a fixed value
for all energies above the considered energy threshold, and no energy
dependency above the threshold is taken into account. In our astrophysical
model, however, the signal fraction function
$f_{S}(f_{0},E^{e}_{\mathrm{det}})$ is still determined only by one free model
parameter $f_{0}$, but the energy dependency of the signal contribution is
generated by propagation effects. The same argument holds for the magnetic
field blurring, which is described by one parameter only, $\delta_{0}$. But,
because the rigidity-dependency of the blurring is taken into account, the
mean blurring decreases consequently when the rigidity increases with
increasing energy. This allows us to accumulate the signal with one likelihood
function summing over all energy bins (see section 3.1.2), instead of relying
on an energy-threshold scan that would have to be penalized for.

###  2.5 Consideration of experimental systematic uncertainties

Uncertainties on the detector response can induce systematic uncertainties on
the observables. The impact of these experimental systematic uncertainties on
the analysis results and fitted model parameters can be investigated by
introducing nuisance parameters $\nu_{E/X_{\mathrm{max}}}$ in the
astrophysical model, which are here given in units of standard scores.

The energy scale uncertainty can have a significant impact on the fitted model
parameters. For the FD, and hence also for the cross-calibrated SD, the
resolution is almost energy-independent and around $14\%$ [52]. A systematic
shift of the energy scale, parameterized by the nuisance parameter $\nu_{E}$,
corresponds to a shift of the bin content in the detected energy bins.

The shower maximum depth distributions are influenced by the detector
acceptance, resolution, and scale as discussed in section 2.4.2. There are
systematic uncertainties on the parameterizations of the acceptance and
resolution, but the effect on the observables is very minor as described in
Refs. [10, 8], so they are neglected here. This is mostly due to a cut on the
field of view of the FD applied on the measured data, keeping it within a
high-acceptance region. The $X_{\mathrm{max}}$ scale uncertainty is influenced
by effects like the calibration and the atmospheric conditions, leading to an
energy-dependent uncertainty between $\pm 6$ g/cm2 and $\pm 10$ g/cm2 [55]. It
was shown in [8, 10] that the hadronic interaction model used in the Gumbel
distributions can have a major impact on the fit results, and this can
partially be parameterized by a shift of the $X_{\mathrm{max}}$ scale. The
$X_{\mathrm{max}}$ scale uncertainty is included in the astrophysical model by
shifting the $X_{\mathrm{max}}$ values used for the evaluation of the Gumbel
distributions in eq. 2.11 by $\nu_{X_{\mathrm{max}}}$, taking into account the
energy dependence of the scale uncertainty.

Note that, in contrast with [10], in this work we shift the modeled spectrum
and $X_{\mathrm{max}}$ instead of the data. This is because here our model
also contains the forward-folding of the detector effects (instead of the
unfolding approach used in [10]). However, for better comparability, the same
convention as in [10] is applied to define the shift direction (a “positive”
versus a “negative” shift). Hence, for example for $X_{\mathrm{max}}$, a
negative shift of the $X_{\mathrm{max}}$ scale (implying that the true
composition is actually heavier than measured) is represented here by shifting
the mean modeled $X_{\mathrm{max}}$ to larger values. As the underlying
uncertainties concern the data and not the model, in all figures visualizing
the effects of experimental systematic uncertainties (fig. 6 and fig. 7) we
shift the data points instead of the model.

Experimental systematic uncertainties on the arrival directions are
significantly smaller than the size of the magnetic field blurring, and hence
they are neglected for this work.

##  3 Fitting procedure and data sets

The astrophysical model described in section 2 can be used to infer the values
of the free model parameters by comparing the predicted observables to the
data measured with the Pierre Auger Observatory. The inference methods
including the likelihood function will be discussed in section 3.1, and the
measured data sets will be presented in section 3.2.

###  3.1 Inference of model parameters

The astrophysical model contains several free parameters ($\gamma$,
$R_{\mathrm{cut}}$, $I_{A}$, $f_{0}$, $\delta_{0}$, $Q_{0}$) $\in\Theta$ as
well as the nuisance parameters $\nu_{E}$ and $\nu_{X_{\mathrm{max}}}$.

####  3.1.1 Fitting techniques

The inference of the model parameters is performed with two methods. The first
one is a gradient-based minimizer based on scipy [66]. It determines the best
fit according to the likelihood function $\mathcal{L}(\Theta)$ (maximum-
likelihood estimate or MLE). The second method is a Bayesian Markov-Chain-
Monte-Carlo (MCMC) sampler which is used to investigate uncertainties on the
fit parameters via posterior distributions $p(\Theta|d,\mathcal{M})$, where
$d$ denotes the data and $\mathcal{M}$ the model. For the estimation of the
posterior distributions two ingredients are needed, the likelihood function
and the prior distributions $p(\Theta|\mathcal{M})$, according to Bayes
theorem:

$p(\Theta|d,\mathcal{M})=\frac{\mathcal{L}(\Theta)\
p(\Theta|\mathcal{M})}{p(d|\mathcal{M})}\propto p(\Theta|\mathcal{M})\
\mathcal{L}(\Theta).$ (3.1)

Here, $p(d|\mathcal{M}):=b$ is the Bayesian evidence, the probability of the
model itself given the data. It is not needed for the estimation of the
posterior distributions, but can be used to compare the probabilities of
different models [67, 68]. For the sampling of the posterior distributions, we
use the Sequential Monte Carlo sampler form PyMC [69] which can also estimate
the Bayesian evidence.

####  3.1.2 Likelihood function

The total log-likelihood function is given as a sum of the single likelihood
functions for each independent observable:

$\log\mathcal{L}=\log\mathcal{L}_{E}+\log\mathcal{L}_{X_{\mathrm{max}}}+\log\mathcal{L}_{\mathrm{ADs}}$
(3.2)

where “ADs” is used for the arrival directions. The likelihood for the energy
spectrum $\mathcal{L}_{E}$ is given as a Poissonian [8]

$\log\mathcal{L}_{E}=\sum_{e}\Big{(}n^{e}\log(\mu^{e})-\log(n^{e}!)-\mu^{e}\Big{)}.$
(3.3)

Here, $n^{e}$ is the measured number of events per energy bin $e$, and
$\mu^{e}=\sum_{e^{\prime}}R^{e,e^{\prime}}\mu^{e^{\prime}}$ gives the number
of events per energy bin $e$ following eq. 2.10, predicted by the model after
detector effects. The minimum energy for the combined fit of energy spectrum
and shower maximum depth distributions above the ankle with a homogeneous
background model only was set to
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=18.7$ in previous works [8, 70]. In
Ref. [10] it was shown that an additional sub-ankle component can still have
an influence around $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=18.7$, and that
the extragalactic component only completely dominates the flux above $\sim 10$
EeV. For that reason, the minimum energy is set to
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.0$ in this work. Nevertheless, it
is enforced that the flux predicted by the astrophysical model between
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=18.7$ and
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.0$ does not exceed the measured
one by using a one-sided Poissonian likelihood for these energy bins.

The shower maxima $X_{\mathrm{max}}$ information is binned into the histogram
$\mu^{\tilde{e},x}$ following eq. 2.12. Due to the small statistics of the
$X_{\mathrm{max}}$ measurements with the FD at the highest energies, the bins
above $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.6$ are combined into one,
consequently denoted by $\tilde{e}$ instead of $e$. Since the energy spectrum
information is already incorporated by the energy likelihood and should not
have an influence twice, the $X_{\mathrm{max}}$ likelihood is given as a
multinomial instead of a Poissonian

$\mathcal{L}_{X_{\mathrm{max}}}=\prod_{\tilde{e}}n^{\tilde{e}}!\prod_{x}\
\frac{(\mu^{\tilde{e},x})^{n^{\tilde{e},x}}}{n^{\tilde{e},x}!}.$ (3.4)

For the $X_{\mathrm{max}}$ distributions, the minimum energy is also set to
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.0$.

For the arrival directions, a similar likelihood function as in Ref. [6] is
used. The only difference is that the events are binned into the detected
energy bins $e$, and hence that the modeled energy-dependent probability maps
$\mathrm{pdf}^{e,p}$ from eq. 2.17 are used. Then, the pdf value is read out
at the arrival directions of the measured data $n^{e,p}$, also binned in
energy and pixels. This leads to the likelihood function

$\mathcal{L_{\text{ADs}}}=\prod_{e}\prod_{p}\mathrm{pdf}^{e,p}\,n^{e,p}.$
(3.5)

For the arrival directions, a minimum energy of
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.2$ is used which is higher than
that for energy and $X_{\mathrm{max}}$ in the likelihood. The reason for this
is mainly that at $E=10^{19.0}$ eV the dipole [3] is very prominent. With our
simplified treatment of the foreground and background sources, the dipole is
not expected to be described properly by the astrophysical model.
Additionally, with a minimum energy of $E=10^{19.2}\,\text{eV}\approx 16$ EeV,
this analysis is in line with the other recent arrival directions analyses we
performed [2, 71], and can use the same data set for better comparison.

As the experimental systematic uncertainties on the energy and
$X_{\mathrm{max}}$ scale, $\nu_{E}$ and $\nu_{X_{\mathrm{max}}}$, are
expressed in units of standard deviations, the nuisance parameters $\nu$ are
expected to follow a Gaussian likelihood with mean $0$ and standard deviation
$1$

$\mathcal{L}_{\mathrm{syst}}=\frac{1}{\sqrt{2\pi}}\exp(-\nu^{2}/2).$ (3.6)

When the experimental systematic uncertainties are considered in the fit,
$\mathcal{L}_{\mathrm{syst}}$ has to be included in the total likelihood (eq.
3.4).

For the Poissonian energy and the multinomial $X_{\mathrm{max}}$ likelihood, a
deviance $D$ can be defined which characterizes the goodness of fit [72]. It
is given as twice the negative likelihood ratio between the fitted model and
the saturated model $\mathcal{L}^{\mathrm{sat}}$ which describes the data
perfectly. For the arrival directions, the definition of a deviance is not
necessarily useful. This is because a saturated model would not have a
physical meaning: setting the arrival probability histogram
$\mathrm{pdf}^{e,p}$ equal to the measured arrival directions $v^{e,p}$ would
lead to a map of only very sparsely filled pixels, which is not in agreement
with the physics of UHECR sources and the assumptions of the astrophysical
model. For that reason, in the following, the deviance will be stated only for
energy, $X_{\mathrm{max}}$, and systematics, and the likelihood value will be
given for the arrival directions.

####  3.1.3 Prior distributions

For the spectral index $\gamma\in[-4,3]$ and the logarithmic rigidity cutoff
$\log_{10}(R_{\mathrm{cut}}/\mathrm{V})$ $\in[18.0,20.5]$, flat priors are
used as in previous work [8]. The lower border for $\gamma$ is extended from
$-3$ to $-4$ in this work because harder best-fit values for the spectral
index were observed than in Ref. [8]. The signal fraction prior is flat
$f_{0}\in[0,50\%]$. Before, we also tested the range $f_{0}\in[0,100\%]$, but
found that models with very large signal fractions cannot recreate the
findings of Refs. [6, 2], so the maximum allowed signal fraction was
accordingly set to 50%. The blurring prior is also flat and can have values
$\delta_{0}\in[0^{\circ},86^{\circ}$]. The upper border value (1.5 rad)
corresponds, for reasonable UHECR compositions (e.g. nitrogen at $\approx
40~{}\mathrm{EeV}$), to an almost isotropic distribution of arrival directions
at Earth.

For the elemental fractions, the sum over the five representative fractions
has to return unity. This circumstance can be included in the prior
distributions by using the unit simplex [73]. Instead of using the elemental
fractions $a_{A}$ defined below the cutoff for the fit as in Ref. [8], in this
work the integral fractions $I_{A}\in[0,1]$ (eq. 2.3) will be utilized. In the
case of very hard spectral indices even extremely small fractions $a_{A}$ can
still lead to large integrated fractions, and hence large total emissions of
the respective element. It was tested that the sensitivity of the fit in that
case is limited due to the values of $a_{A}$ spanning over multiple orders of
magnitude. This can introduce an unwanted bias on the posterior distributions.

The generation rate normalization $Q_{0}$ is a special case as it only
contributes to the energy likelihood, so it can be deduced simply by
normalizing the number of events predicted by the model to the one in the data
(uninformative prior).

Due to the Gaussian likelihood definition for the systematic uncertainties
(eq. 3.6), their expected distribution must not be taken into account again in
the prior, so a uniform prior between $-4\sigma$ and $4\sigma$ is used to keep
the nuisance parameters in a reasonable range.

###  3.2 The data sets

The Pierre Auger Observatory [1] is located in Malargüe, Argentina. It uses a
hybrid detector design. For the energy spectrum, the measurements of the
surface detector taken between 01/2004 and 08/2018 with zenith angles below
$60^{\circ}$ from Ref. [52] are used. The exposure of the data set is
$60{,}426$ km2 sr yr. Above $10^{18.7}$ eV, the data set contains $55{,}730$
events. No unfolding of the detector effects is necessary for this study
because the astrophysical model contains a forward folding of these effects as
described in section 2.4.1.

The shower maximum depth distributions are taken from the measurements of the
fluorescence detector of the Pierre Auger Observatory [74]. Above $10^{19.0}$
eV, the data set contains 1022 events which are binned into a 2-dimensional
histogram in energy and $X_{\mathrm{max}}$ as described in section 3.1.2.

For the arrival directions, the data set is based on the one from [2],
containing the largest set of events from Phase 1 of the Pierre Auger
Observatory. It consists of SD data measured between 01/01/2004 and
31/12/2020. In total, the data set contains $12{,}606$ events above the
minimum energy of 16 EeV, so it also extends to lower energies than the one
described in [2] ($>32$ EeV). The data set combines both vertical (zenith
angle smaller than $60^{\circ}$) and inclined events (zenith angle larger than
$60^{\circ}$ and up to $80^{\circ}$). The exposure amounts to $95{,}700$ km2
sr yr for the vertical subset and $26{,}300$ km2 sr yr for the inclined
subset.

##  4 Fit results

In the following, the fit results will be presented, starting with the models
containing SBGs or Centaurus A in section 4.1. The results of the $\gamma$-AGN
model will be discussed separately in section 4.2. For comparison, the results
of the reference models with only homogeneous background are given in appendix
A.

###  4.1 Starburst galaxy and Centaurus A models

####  4.1.1 Fitted model parameters

The fit results for the SBG and Centaurus A models are given in table 1. Here,
the best fit (MLE) is stated including the respective deviance and likelihood
values. Additionally, the posterior mean and the borders of the highest
posterior density interval666Defined as the shortest interval in which $68\%$
of the posterior resides. Note that these borders can be asymmetric with
respect to the posterior mean for not fully-Gaussian posterior distribution.
So, the values stated in the tables should not be confused with $1\sigma$
uncertainties. are provided. As an example for posterior distributions, fig. 1
depicts them for the SBG model. For all models, a strong correlation between
the spectral index $\gamma$ and the maximum rigidity
$\log_{10}(R_{\mathrm{cut}}/\mathrm{V})$ is visible as in Ref. [8], as well as
a correlation between the signal fraction $f_{0}$ and the blurring
$\delta_{0}$ as expected from Ref. [6]. The integral composition fractions are
typically almost uncorrelated with the other parameters.

|  Cen A, $\boldsymbol{m=0}$ (flat) |  Cen A, $\boldsymbol{m=3.4}$ (SFR) |  SBG, $\boldsymbol{m=3.4}$ (SFR)
---|---|---|---
| posterior | MLE | posterior | MLE | posterior | MLE
$\gamma$ | $-1.67_{-0.47}^{+0.48}$ | $-2.21$ | $-3.09_{-0.24}^{+0.23}$ | $-3.05$ | $-2.77_{-0.29}^{+0.27}$ | $-2.67$
$\log_{10}(R_{\mathrm{cut}}$/V) | $18.23_{-0.06}^{+0.04}$ | 18.19 | $18.10_{-0.02}^{+0.02}$ | 18.11 | $18.13_{-0.02}^{+0.02}$ | 18.13
$f_{0}$ | $0.16_{-0.14}^{+0.06}$ | $0.028$ | $0.05_{-0.03}^{+0.01}$ | $0.028$ | $0.17_{-0.08}^{+0.06}$ | 0.19
$\delta_{0}/^{\circ}$ | $56.5_{-12.8}^{+29.4}$ | 16.5 | $27.6_{-16.3}^{+2.7}$ | 16.8 | $22.2_{-4.0}^{+5.3}$ | 24.3
$I_{\mathrm{H}}$ | $5.9_{-1.7}^{+2.5}\times 10^{-2}$ | $7.1\times 10^{-2}$ | $8.3_{-8.3}^{+2.0}\times 10^{-3}$ | $1.6\times 10^{-5}$ | $6.4_{-6.4}^{+1.3}\times 10^{-3}$ | $4.3\times 10^{-5}$
$I_{\mathrm{He}}$ | $2.3_{-0.5}^{+0.3}\times 10^{-1}$ | $1.9\times 10^{-1}$ | $1.3_{-0.2}^{+0.2}\times 10^{-1}$ | $1.4\times 10^{-1}$ | $1.7_{-0.4}^{+0.3}\times 10^{-1}$ | $1.8\times 10^{-1}$
$I_{\mathrm{N}}$ | $6.3_{-0.3}^{+0.3}\times 10^{-1}$ | $6.2\times 10^{-1}$ | $7.4_{-0.3}^{+0.3}\times 10^{-1}$ | $7.3\times 10^{-1}$ | $7.4_{-0.3}^{+0.3}\times 10^{-1}$ | $7.4\times 10^{-1}$
$I_{\mathrm{Si}}$ | $6.5_{-3.3}^{+3.6}\times 10^{-2}$ | $9.9\times 10^{-2}$ | $9.2_{-2.3}^{+3.2}\times 10^{-2}$ | $1.1\times 10^{-1}$ | $5.7_{-3.1}^{+2.5}\times 10^{-2}$ | $5.4\times 10^{-2}$
$I_{\mathrm{Fe}}$ | $1.6_{-1.0}^{+0.7}\times 10^{-2}$ | $2.0\times 10^{-2}$ | $2.5_{-0.9}^{+0.8}\times 10^{-2}$ | $2.3\times 10^{-2}$ | $2.5_{-0.9}^{+0.8}\times 10^{-2}$ | $2.3\times 10^{-2}$
$\boldsymbol{\log b}$ | $-264.0\pm 0.2$ |  | $-272.6\pm 0.2$ |  | $-266.9\pm 0.1$ |
$\boldsymbol{D_{E}}\ (N_{J}=14)$ |  | 22.3 |  | 28.5 |  | 33.3
$\boldsymbol{D_{X_{\mathrm{max}}}}\ (N_{{X}_{\mathrm{max}}}=74)$ |  | 124.9 |  | 130.6 |  | 126.2
$\boldsymbol{D}$ |  | 147.2 |  | 159.1 |  | 159.5
$\boldsymbol{\log}$ $\boldsymbol{\mathcal{L}_{\mathrm{ADs}}}$ |  | $-10.5$ |  | $-10.4$ |  | $-13.3$
$\boldsymbol{\log}$ $\boldsymbol{\mathcal{L}}$ |  | $-239.1$ |  | $-245.1$ |  | $-242.4$

Table 1: Fit results for the Centaurus A and SBG models. The best-fit (MLE)
parameters and the corresponding deviance $D$ and log-likelihood values
$\log\mathcal{L}$ are stated. Also, the posterior mean and highest posterior
density interval from the MCMC sampler are given, as well as the logarithmic
Bayesian evidence $\log b$ (see text). Note that in this and all following
tables the value for the arrival directions likelihood
$\log\mathcal{L}_{\mathrm{ADs}}^{\mathrm{ref}}$ for fully isotropic arrival
directions (as in the reference models, see table 4) has been set to 0 for
better readability of the values of $\log b$, $\log\mathcal{L}$, and
$\log\mathcal{L}_{\mathrm{ADs}}$. Figure 1: Posterior distributions for the
SBG model. The pink cross marks the MLE, the black vertical line the posterior
mean.

From the reference models without catalog sources (see appendix A and fig. 8
and compare to Ref. [10]) it is clear that a strong source evolution with
$m\gtrsim 5$ is disfavored by the data. This is due to the great amount of
low-energy secondaries from photodisintegration which overshoot the measured
spectrum at low energy. They are produced over the large propagation distances
when the background sources emit predominantly at large redshifts for a strong
source evolution. The deviance for a flat evolution is smaller than that for a
SFR one by around $\Delta D\approx 10$. The spectral index is very hard
already for flat evolution, but can become as hard as $\gamma\approx-3$ for
the SFR case. The maximum rigidity is relatively small with
$\log_{10}(R_{\mathrm{cut}}/\mathrm{V})\approx 18.2$ V.

####  4.1.2 Influence of the catalog sources on the observables at Earth

When adding Centaurus A as a single source on top of the background, the
injection spectrum softens slightly compared to the reference model for both
tested evolutions (a strong evolution with $m=5.0$ still leads to a poor
description of spectrum and composition when Centaurus A is included, so it is
not further discussed here, see however fig. 8). The best-fit signal fraction
(MLE) associated with Centaurus A is estimated to be around $f_{0}\approx 3\%$
for both evolutions. But, especially for the flat evolution case, the
posterior distributions for signal fraction and blurring are quite broad. This
leads to a posterior mean value of $f_{0}=0.16^{+0.06}_{-0.14}$ in the case of
flat evolution, and $f_{0}=0.05^{+0.01}_{-0.03}$ for SFR evolution.

The modeled spectra on Earth are shown in fig. 2. For both evolutions, the
spectra look quite similar so that conclusions can be drawn independently of
the model evolution of the background sources. It is visible how the nearby
source Centaurus A contributes mainly to the flux at higher energies. The
best-fit signal contribution to the total flux rises from $\sim 3\%$ at $40$
EeV to $\sim 10\%$ at 100 EeV. The uncertainty on the signal fraction is large
as described above, so contributions up to $\sim 50\%$ at 100 EeV are still
within the 68% highest posterior density interval. The composition at the
sources is dominated by nitrogen. This leads to a nitrogen flux at Earth up to
very high energies, which has a rigidity of around $R\simeq 10$ eV at $70$
EeV. For this, the model predicts a blurring of around
$\delta_{0}=20^{\circ}$. The contribution of lighter elements to the flux at
the highest energies (both from the nearby sources and from the background)
turns out to be very small, and is completely negligible above $\sim 30$ EeV
where the helium contribution cuts off. Most of the light elements present
below those energies are actually of secondary origin, produced in the
photodisintegration of the heavier primaries from the background sources. The
first and second moments of the modeled shower maximum depth distributions are
depicted in fig. 3. The model describes the measured data fairly well above
the minimum energy employed in the fit. It is however visible that the model
mean mass is slightly heavier than the measured one, which can be resolved by
including the experimental $X_{\mathrm{max}}$ scale uncertainty as a nuisance
parameter in the model, as will be shown in sec. 4.1.3.

Figure 2: Modeled spectra on Earth for the Centaurus A model with $m=0$
(left), with $m=3.4$ (middle) and the SBG model with $m=3.4$ (right). In this
and all following figures the thick lines indicate the best-fit and the thin
lines are drawn from the posterior distribution demonstrating the uncertainty.
The markers represent the measured data of the Pierre Auger Observatory. The
grey area symbolizes the energy bins which are not (fully) included in the fit
(see section 3.1.2). In the upper row the total spectrum is depicted with the
different element contributions in different colors. The contribution by the
source catalog is indicated by the dashed line. The lower row shows the
individual spectrum of the strongest source in the respective catalog.

Figure 3: First two moments of the $X_{\mathrm{max}}$ distributions for the
Centaurus A model with $m=0$ (left), with $m=3.4$ (middle) and the SBG model
with $m=3.4$ (right).

The best-fit arrival directions for three representative energy bins are shown
in fig. 4. The growing flux contribution of Centaurus A with the energy, as
well as the decreasing blurring, is visible.

Figure 4: Modeled arrival directions pdf (eq. 2.17), left for Centaurus A
model with $m=0$, middle for Centaurus A model with $m=3.4$ and right for SBG
model with $m=3.4$. The energy bins $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})$
=19.3 (upper row), $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})$ =19.6 (middle
row) and $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})$ =19.9 (lower row) are shown
as examples. The catalog contribution and with that the level of anisotropy
rises with the energy, while the overall blurring decreases. Additionally, the
contribution of individual sources depends on the energy through their
distances and flux weights. The stars indicate the directions of the source
candidates with the size scaling with the relative flux contribution before
the observatory exposure is applied.

For the SBG model with $m=3.4$, the spectral index has softened visibly
compared to the reference model with $m=3.4$. The softening decreases the
number of high-energy particles emitted at the background sources. This
compensates the increased number of high-energy particles from the catalog
sources which can reach Earth easily due to the short propagation distances.
The signal fraction is around $f_{0}\approx 20\%$. In fig. 2, it is again
visible how the catalog contribution becomes relevant at higher energies,
where it reaches values of up to $\approx 40\%$ at 100 EeV. Additionally, fig.
2 displays the individual spectrum of the strongest source in the SBG catalog,
NGC 4945. The modeled spectrum looks rather similar to the one of Centaurus A,
only with smaller uncertainties due to the additional constraints from the
other candidates in the SBG catalog. The two source candidates NGC 4945 and
Centaurus A are located in similar directions and distances of around $3.5$
Mpc. Hence, the contribution of that sky region is modeled consistently,
independent from the number of other subdominant candidate sources in the
catalog. This is also visible in the arrival directions in fig. 4, where both
the size of the blurring and the overall anisotropy level is similar for both
NGC 4945 and Centaurus A.

To verify explicitly that the model describes the overdensity in the region of
Centaurus A and NGC 4945, we investigate the spectrum of all events in a
circular region with variable angular size centered on the direction of the
respective source candidate. The size of the region is a trade-off between
wanting to fully contain the contribution of the source candidate, and
minimizing the contamination by background and neighboring candidates. The
modeled and measured (calculated from the energy spectrum data set described
in section 3.2) spectra for one example selection angle of $20^{\circ}$
(radius) are displayed in fig. 5. Here, it is visible how the model predicts
an increased flux for the region around the source candidate. In particular,
more high-energy events are expected from the close-by candidate. Due to this,
the model spectrum from the candidate region (dashed lines) agrees better with
the measured flux (markers) than the whole-sky model spectrum (solid lines).
For all tested angles between $5^{\circ}$ and $30^{\circ}$, the reduced
$\chi^{2}$ for the candidate region model is smaller than that of the whole-
sky model by around $10\%$ to $40\%$. With increased statistics in the
candidate region, a more detailed assessment of this effect could be performed
in the future.

Figure 5: Modeled spectra on Earth inside a circular region of angular size
$20^{\circ}$ (radius) around the direction of Centaurus A (left and middle) or
NGC 4945 (right). Displayed are the Centaurus A model with $m=0$ (left), with
$m=3.4$ (middle), and the SBG model with $m=3.4$ (right). The colors are the
same as in fig. 2. The markers represent the measured data from the Pierre
Auger Observatory inside the $20^{\circ}$ region, the dashed lines the best-
fit model in that region, and the solid lines the full best-fit model (same as
in fig. 2). Note that the spectra in the $20^{\circ}$ region have not been
fitted explicitly (see section 3.1.2), so this serves as a cross check of the
model.

Note that for a source such as Centaurus A or NGC 4945 at about 4 Mpc
distance, one can estimate from the inferred spectral shape and the derived
composition that the CR luminosity above 10 EeV should be
$L_{>10\,\mathrm{EeV}}\simeq(f_{\mathrm{source}}/0.05)\times 10^{39}$ erg/s,
with $f_{\mathrm{source}}$ being the overall fractional contribution from the
specified source to the flux at 40 EeV (directly corresponding to $f_{0}$ for
the case of the Centaurus A model, and to around
$f_{0}\times\omega_{\mathrm{flux}}^{\mathrm{NGC4945}}$ for the SBG model). For
a scenario like the one of the SBGs, one would expect that the remaining
catalog sources should have CR luminosities approximately scaling as the
respective adopted flux weights.

####  4.1.3 Impact of experimental systematic uncertainties

The effect of experimental systematic uncertainties on the fit result is
tested by including the uncertainty on the energy and $X_{\mathrm{max}}$
scales as nuisance parameters in the model (see section 2.5). The fit results
for that are given in table 2. The preferred shift on the energy scale is
estimated to be around $-1.2\sigma$ for the model with flat evolution
(implying that the true energy is lower than measured) and almost no shift is
favored for the models with SFR evolution. In general, the shift of the energy
scale has a smaller impact on the deviance than the shift of the
$X_{\mathrm{max}}$ scale. Here, a negative shift is preferred, which is around
$-1\sigma$ for $m=0$ and $-1.5\sigma$ for $m=3.4$, in accordance with results
for the reference models (see appendix A). This shifts the $X_{\mathrm{max}}$
scale up by around $10-15\,\text{g/cm}^{2}$ (see section 2.5), which is in
line with the findings from Ref. [75]. According to these results, the true
$X_{\mathrm{max}}$ values are lower than reconstructed, and the true
composition is correspondingly heavier. The main effect on the model
parameters is a softening of the spectral index to around $\gamma\approx-1$
for both evolutions. The softening is also visible in the modeled spectra on
Earth in fig. 6. The modeled signal fraction and blurring are almost
completely independent of the systematic uncertainties on the energy and
$X_{\mathrm{max}}$ scales, so the contribution from the catalog sources is
also similar to the one without systematics displayed in fig. 2. The modeled
arrival directions look similar to the ones without systematics (fig. 4), so
they are not shown again. Also, the different element contributions are almost
unchanged, so still no elements lighter than the nitrogen group are expected
at the highest energies. The fluctuations between posterior draws have
increased, as is visible in fig. 6, which is because of the added freedom
allowed by the two additional fit parameters for the experimental systematic
uncertainties.

As displayed in fig. 7, the measured mean shower maximum depth is now
substantially better described than without the shift of the
$X_{\mathrm{max}}$ scale.

|  Cen A, $\boldsymbol{m=0}$ (flat) |  Cen A, $\boldsymbol{m=3.4}$ (SFR) |  SBG, $\boldsymbol{m=3.4}$ (SFR)
---|---|---|---
| posterior | MLE | posterior | MLE | posterior | MLE
$\gamma$ | $-0.89_{-0.33}^{+0.37}$ | $-0.65$ | $-1.19_{-0.39}^{+0.45}$ | $-1.41$ | $-1.02_{-0.36}^{+0.43}$ | $-1.25$
$\log_{10}(R_{\mathrm{cut}}$/V) | $18.20_{-0.05}^{+0.04}$ | 18.23 | $18.21_{-0.05}^{+0.04}$ | 18.20 | $18.24_{-0.06}^{+0.04}$ | 18.22
$f_{0}$ | $0.07_{-0.05}^{+0.01}$ | $0.029$ | $0.07_{-0.05}^{+0.01}$ | $0.031$ | $0.19_{-0.11}^{+0.07}$ | 0.23
$\delta_{0}/^{\circ}$ | $30.5_{-20.2}^{+2.0}$ | 14.4 | $27.4_{-17.0}^{+4.2}$ | 14.3 | $18.8_{-3.6}^{+5.9}$ | 21.9
$I_{\mathrm{H}}$ | $5.8_{-2.6}^{+2.9}\times 10^{-2}$ | $4.2\times 10^{-4}$ | $1.2_{-1.2}^{+0.2}\times 10^{-2}$ | $3.0\times 10^{-4}$ | $1.2_{-1.2}^{+0.1}\times 10^{-2}$ | $1.0\times 10^{-4}$
$I_{\mathrm{He}}$ | $2.7_{-0.4}^{+0.4}\times 10^{-1}$ | $3.5\times 10^{-1}$ | $9.9_{-2.9}^{+3.8}\times 10^{-2}$ | $1.2\times 10^{-1}$ | $1.1_{-0.4}^{+0.3}\times 10^{-1}$ | $1.4\times 10^{-1}$
$I_{\mathrm{N}}$ | $5.6_{-0.4}^{+0.4}\times 10^{-1}$ | $5.0\times 10^{-1}$ | $6.7_{-0.7}^{+0.7}\times 10^{-1}$ | $6.8\times 10^{-1}$ | $7.2_{-0.6}^{+0.6}\times 10^{-1}$ | $7.3\times 10^{-1}$
$I_{\mathrm{Si}}$ | $9.0_{-3.4}^{+3.9}\times 10^{-2}$ | $1.4\times 10^{-1}$ | $1.5_{-0.6}^{+0.5}\times 10^{-1}$ | $1.6\times 10^{-1}$ | $1.2_{-0.5}^{+0.5}\times 10^{-1}$ | $9.8\times 10^{-2}$
$I_{\mathrm{Fe}}$ | $2.3_{-1.2}^{+0.9}\times 10^{-2}$ | $1.8\times 10^{-2}$ | $5.1_{-1.8}^{+1.5}\times 10^{-2}$ | $4.4\times 10^{-2}$ | $4.7_{-1.7}^{+1.3}\times 10^{-2}$ | $3.8\times 10^{-2}$
$\nu_{E}/\sigma$ | $-1.24_{-0.50}^{+0.68}$ | $-1.35$ | $0.23_{-0.60}^{+0.42}$ | $0.13$ | $0.35_{-0.65}^{+0.44}$ | $0.40$
$\nu_{X\mathrm{max}}/\sigma$ | $-0.94_{-0.24}^{+0.29}$ | $-0.97$ | $-1.60_{-0.25}^{+0.30}$ | $-1.45$ | $-1.55_{-0.25}^{+0.26}$ | $-1.33$
$\boldsymbol{\log b}$ | $-254.6\pm 0.1$ |  | $-264.5\pm 0.2$ |  | $-258.6\pm 0.2$ |
$\boldsymbol{D_{\mathrm{syst}}}$ |  | 2.8 |  | 2.1 |  | 1.9
$\boldsymbol{D_{E}}\ (N_{J}=14)$ |  | 13.6 |  | 21.9 |  | 25.3
$\boldsymbol{D_{X_{\mathrm{max}}}}\ (N_{{X}_{\mathrm{max}}}=74)$ |  | 107.4 |  | 113.6 |  | 112.7
$\boldsymbol{D}$ |  | 123.8 |  | 137.7 |  | 139.9
$\boldsymbol{\log}$ $\boldsymbol{\mathcal{L}_{\mathrm{ADs}}}$ |  | $-9.4$ |  | $-9.5$ |  | $-13.5$
$\boldsymbol{\log}$ $\boldsymbol{\mathcal{L}}$ |  | $-228.51$ |  | $-235.3$ |  | $-232.4$

Table 2: Centaurus A and SBG models with experimental systematic uncertainties
included as nuisance parameters

Figure 6: Same as fig. 2 (upper row) but including experimental systematic
uncertainties as nuisance parameters.

Figure 7: Same as fig. 3 but including experimental systematic uncertainties
as nuisance parameters.

####  4.1.4 Evaluation of the model performances

How well each of the tested models describes the measured data can be
quantified for example by comparing the values of the likelihood function for
the best-fit model parameters. An overview of all likelihood values is given
in fig. 8, including the contributions by the different observables. Here, the
general trends that have been described above are visible, for example, that
the models with strong source evolution are disfavored. This is mostly due to
a poor description of the energy spectrum (blue). Also, when comparing the
models with and without systematics, one can see that the likelihood ratio
improves consistently when systematic shifts are allowed, mainly due to a
better description of the shower maximum depth distributions (grey), but also
the energy spectrum shows an improvement for the cases with strong evolution.
The arrival direction likelihood (green) is almost independent of the source
evolution and the systematics, and is always the largest for the SBG catalog.

For a more quantitative comparison, we use the test statistic calculated as 2
times the likelihood ratio between a model and the respective reference model
with the same evolution and (no) systematics:

$\mathrm{TS}_{\mathrm{tot}}=\sum_{\mathrm{obs}=E,X_{\mathrm{max}},\mathrm{ADs}}2(\log\mathcal{L}^{m=x}-\log\mathcal{L}_{\mathrm{ref}}^{m=x})^{\mathrm{obs}}.$
(4.1)

Hence, the test statistic describes the improvement of adding a specific
catalog to a model compared to just homogeneous sources. The values for the
test statistic of each model are given in table 3. Apparently, the arrival
directions observable contributes most to the total test statistic, and is
hence very important for the differentiation of different models. This is as
expected from a simulation study [17]. But, from this analysis, it also
becomes clear that the exact values of TS should be treated with caution as
they can vary considerably and depend on e.g. the distribution of the arrival
directions in an energy bin.

Figure 8: Overview of the likelihood values for all tested models: for better
comparability, for all models 2 times the likelihood ratio with the reference
model with $m=0.0$ is given. The white bar with black edges gives the total
likelihood, and the four smaller bars inside give the contributions of the
observables and the systematics to the total likelihood.

|  Cen A, $m=0.0$ |  Cen A, $m=3.4$ |  SBG, $m=3.4$ | $\gamma$AGN, $m=5.0$ | $\gamma$AGN+EGMF**, $m=5.0$
---|---|---|---|---|---
|  | \+ syst |  | \+ syst |  | \+ syst |  | \+ syst |  | \+ syst
TStot | 22.8 | 17.3 | 22.2 | 19.1 | 27.6 | 25.6 | 23.9* | 9.8* | 34.3* | 33.2*
TSE | $-0.1$ | $-1.4$ | $-0.4$ | $-1.1$ | $-5.2$ | $-4.5$ | 26.8 | 3.9 | 18.2 | 8.4
TS${}_{X_{\mathrm{max}}}$ | 1.9 | 0.2 | 1.8 | $1.0$ | $6.2$ | $2.0$ | $-0.8$ | 6.4 | 4.4 | 14.7
TSADs | 20.9 | 18.7 | 20.8 | 19.0 | $26.6$ | 27.1 | $-2.1$ | $-3.0$ | 11.7 | 8.6

Table 3: Overview of test statistic values for the tested models. In the case
with systematics, the contribution from eq. 3.6 is taken into account
additionally in TStot. * Note that the test statistic of each model is always
calculated with respect to the reference model with the same source evolution
and (no) systematics. This implies that e.g. the test statistics of the
$\gamma$-AGN model is large only due to the poor fit of the reference model
with $m=5.0$ (compare to fig. 8). ** Note also that the model with EGMF has an
additional fit parameter $\beta_{e}$ compared to the other models (see section
4.2).

Additionally, it is important to note that the test statistic only quantifies
the best-fit result. As indicated by the broad posterior distributions (see
fig. 1), alternative combinations of model parameters can lead to similar
values of the total likelihood function. This means that the SBG model can
also describe the data almost equally well with a smaller signal fraction and
harder spectral index, leading to a better description of the energy spectrum
but simultaneously a worse description of the arrival directions. So, the
values given for the test statistic of individual observables in table 3
depend largely on the best-fit parameter combination.

Due to this, it is often useful to not only compare the test statistic of the
best-fit solution, but also the Bayesian evidence as an additional test. As
introduced in section 3.1.1, it can be used to compare the fit quality of a
whole model, not just the best-fit parameter combination. The values for the
Bayesian evidence $\log b$ for each model are stated in tables 1 and 2. When
comparing the values of the evidence, it is clear that they follow the same
trend as the maximum likelihood value which means that the maximum-likelihood
value / total test statistic is representative of the fit quality of the whole
model. This allows us to focus on the maximum-likelihood values in the
following.

The arrival directions test statistic contributed by each separate energy bin
can be seen in fig. 9. Here, a general trend can be observed, independently of
the systematic effects and of the source evolution. The arrival directions are
well described in the energy bin
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.3$ and the bins around
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.7$ for both the SBG and the
Centaurus A models. Peaks of the arrival direction test statistic at threshold
energies of 40 EeV$\simeq 10^{19.6}~{}\mathrm{eV}$ and 60 EeV$\simeq
10^{19.8}~{}\mathrm{eV}$ have also been observed in Refs. [2, 6]. But, only
with the present analysis which models the energy dependency of each mass
group’s contribution, these can be compared to the modeled spectra of the
strongest sources in the catalog (fig. 2 and fig. 6). Here, one can see that
at $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.3$ the helium contribution has
its peak, while at $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.6$ nitrogen is
predominant, and at $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.8$ silicon
begins to emerge. At these energies, the rigidities correspond to
$R_{\mathrm{He}}=\frac{10^{19.3}}{2}\mathrm{V}\simeq 10~{}\mathrm{EV}$,
$R_{\mathrm{N}}=\frac{10^{19.6}}{7}\mathrm{V}\simeq 6~{}\mathrm{EV}$ and
$R_{\mathrm{Si}}=\frac{10^{19.8}}{14}\mathrm{V}\simeq 5~{}\mathrm{EV}$,
respectively. The similar rigidities of around 6 eV for the nitrogen and
silicon contributions lead to a relatively constant blurring over the whole
energy range which could explain why the test statistic of [6, 2] is so large
(TS${}_{\text{ADs-only}}^{\mathrm{SBG}}\simeq 25.0$) even though the
correlation analysis is performed together for all CRs above an energy
threshold with one fixed blurring. It is important, however, to note that the
energy dependency of the test statistic is subject to large statistical
fluctuations [17], so especially the less pronounced peak at
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.3$ could just be a statistical
fluctuation.

Figure 9: Arrival directions likelihood ratio / test statistic as a function
of the energy bin (marker sits on the bin center) for the tested source
classes and evolutions. The test statistic is always calculated compared to
the reference model with the same evolution and without (left) / with (right)
systematics. The sum over all energy bins corresponds to the values for TSADs
given in table 3.

The energy dependency of the arrival direction test statistic reveals another
interesting finding, which is that the SBG model consistently describes the
measured arrival directions better than the Centaurus A model in the energy
range $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.4-19.6$. For these energies,
the modeled arrival direction pdfs (fig. 4) are very similar for both models
in the region around Centaurus A / NGC 4945, but differ in the southern region
around NGC 253. From table 3 it can be concluded that $\sim 20$ of the total
arrival directions test statistic for the SBG model originates from the
Centaurus A / NGC 4945 region. By removing NGC 253 from the catalog we
explicitly checked that that source contributes around $\sim 4$ to 5 to the
arrival directions test statistic.

###  4.2 The $\boldsymbol{\gamma}$-AGN model

When the $\gamma$-AGN model with source evolution $m=5.0$ is applied to the
data, as expected from the reference model (appendix A) a very hard spectral
index of $\gamma\approx-3.5$ is found in combination with a small maximum
rigidity of $\log_{10}(R_{\mathrm{cut}}/\mathrm{V})\approx 18.1$. The best-fit
signal fraction is around $f_{0}\approx 15\%$. The $\gamma$-AGN catalog is
dominated by a single source between 10 EeV and 100 EeV, the faraway blazar
Markarian 421, which leads to a constant catalog contribution in that energy
range of around 15%.

When comparing the likelihood of the $\gamma$-AGN model with the respective
reference model with $m=5.0$ (fig. 8 and table 3), it is clear that the
inclusion of the catalog sources mainly leads to an improvement of the energy
spectrum description. This is due to the reduced overshooting of the spectrum
at low energies by secondary products from faraway background sources, which
is reduced by the contribution of the closer catalog sources. Overall, the
likelihood of the $\gamma$-AGN model is small, and it does not reach values
larger than the reference models with less evolution. For that reason, the
best-fit parameters will not be discussed in more depth.

Also, it is visible that the arrival directions test statistic is negative.
The best-fit arrival directions are depicted in fig. 10 (upper row). Here, the
dominance of Markarian 421 is clearly visible. Only at the highest energies
around 100 EeV where the statistics of the data set is very small, Centaurus A
starts contributing. The best-fit blurring reaches the maximum value allowed
in the sampling space (see section 3.1.3). This means, if even larger values
of the blurring were accepted, the arrival directions would be modeled as
completely isotropic, leading to
TS${}_{\mathrm{ADs}}^{\gamma-\mathrm{AGN}}=0$. When looking at the energy
dependency of the arrival directions test statistic displayed in fig. 9, one
can see that the negative contribution comes from the smallest energy bins
where only Markarian 421 contributes. Note that the poor description of the
arrival directions with the $\gamma$-AGN model does not depend on the source
evolution of the background sources.

Figure 10: Modeled arrival directions for the $\gamma$-AGN model, left for the
model without and right with distance-dependent blurring. Again, the energy
bins $\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.3$ (upper row),
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})$ =19.6 (middle row) and
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})=19.9$ (lower row) are shown as
examples. Note the different color bar normalizations of the two models.

##### Influence of a distance-dependent blurring:

We tested if a distance-dependent blurring as expected from an extragalactic
magnetic field could improve the agreement between the $\gamma$-AGN model and
the measured data777For the SBG model, we also tested the influence of a
distance-dependent blurring. However, as almost all contributing SBGs are
located at similar distances of $\sim 4$ Mpc, the distance-dependent and
distance-independent blurrings are fully correlated. The likelihood does not
increase by more than $\sim 0.5$ which is compensated by the additional degree
of freedom, so this model is not discussed further here.. For that, we include
an additional fit parameter $\beta_{e}$ in the arrival directions modeling,
and replace eq. 2.14 by the following

$\delta(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}},C^{s})=\frac{\sqrt{\delta_{0}^{2}+\beta_{e}^{2}\
(d^{s}/\mathrm{Mpc})}}{R_{\mathrm{det}}(E^{e}_{\mathrm{det}},A^{k}_{\mathrm{det}})/10\
\mathrm{EV}}.$ (4.2)

The prior for the distance-dependent blurring is flat
$\beta_{e}\in[0^{\circ},11.5^{\circ}]$. This means, a 10 eV particle can at
most obtain a distance-dependent blurring of $\sim 20^{\circ}$ for a source
distance of 3 Mpc. This maximum value was chosen so that the magnetic field
parameters do not become much larger than theoretical limits [42].

The found model parameters for the $\gamma$-AGN model with distance-dependent
blurring do not differ much from the ones without it. The posterior
distribution of $\beta_{e}$ is broad between $\beta_{e}=5^{\circ}$ and the
maximum allowed $\beta_{e}=11.5^{\circ}$, while the distance-dependent
blurring has decreased to $\delta_{0}\simeq 20^{\circ}$. This leads to the
arrival directions displayed in fig. 10. In contrast with the case without
distance-dependent blurring, now the contribution by faraway Markarian 421 is
completely blurred out, and only Centaurus A is visible above
$\log_{10}(E_{\mathrm{det}}/\mathrm{eV})\simeq 19.5$. These arrival directions
now describe the data better than the isotropic reference model, as visible in
fig. 8. Nevertheless, with
TS${}_{\mathrm{ADs}}^{\gamma\mathrm{-AGN+EGMF}}\simeq 10$, it still does not
describe the measured arrival directions as well as the models with only
Centaurus A as a single source do. We tested explicitly if this model can
reproduce the results from Refs. [6, 2] and found that it cannot reach large
enough arrival direction test statistics.

For that reason, we can conclude that the $\gamma$-AGN model where the flux
contribution scales linearly with the $\gamma$-ray intensity, and which thus
favors blazars like Markarian 421, is not supported by the data of the Pierre
Auger Observatory. The measured arrival directions cannot be described
sufficiently well even when including an extragalactic magnetic field with a
strength in the range of theoretical limits. Additionally, as discussed above,
the strong source evolution associated with high-luminosity AGNs leads to an
overproduction of secondaries during propagation, which causes models with
strong evolution to poorly describe the measured energy spectrum.

##  5 Discussion and conclusion

The arrival directions of UHECRs exhibit correlations with catalogs of
powerful extragalactic objects [2, 6]. In this work, these correlations with a
catalog of starburst galaxies and one of jetted active galactic nuclei have
been investigated further using a combined fit of energy spectrum, shower
maximum depth distributions, and arrival directions. The astrophysical model
that is used in the fit includes propagation effects and an energy-dependent
modeling of the contribution of each source depending on its distance and
injection spectrum. Also, a rigidity-dependent magnetic field blurring is
included, which is calculated individually for each source depending on the
element contributions of the respective source on Earth. At the current stage,
no coherent magnetic field deflections are employed in the model, so all
obtained results are under the assumption of mostly turbulent deflections by
the Galactic and extragalactic magnetic fields. In the future, coherent
deflections could also be included in the astrophysical model, hoping also
that a better convergence between different Galactic magnetic field models
will be reached.

The model has several free fit parameters, in particular those describing the
injected spectrum and composition at the sources. Compared to the reference
models with only homogeneous background sources (as used in Refs. [8, 10]),
two additional fit parameters are introduced, the signal fraction of the
nearby sources at fixed energy and the size of the magnetic field blurring
that scales with rigidity. By modeling the arrival directions as a sum of
weighted von Mises–Fisher distributions for each source and each mass number
group in each energy bin, the model naturally produces overdensities which
decrease in size for increasing rigidity. Also, the flux modulations naturally
become stronger with the energy when the contribution from the nearby catalog
sources rises.

The best agreement with the data for models with a background source evolution
following the star-formation rate is found for the SBG model with a hard
spectral index, a signal fraction of $\sim 20\%$ at 40 EeV and a
blurring888Note that in the non-resonant scattering regime, the expected root-
mean-square blurring from turbulent magnetic fields is
$\delta_{\mathrm{rms}}\simeq
20^{\circ}\frac{B}{2.5\,\mathrm{nG}}\frac{10\,\mathrm{EV}}{E/Z}\sqrt{\frac{d_{\mathrm{source}}}{4\,\mathrm{Mpc}}}\sqrt{\frac{l_{\mathrm{c}}}{\mathrm{Mpc}}}$
[76], allowing for conclusions on the corresponding coherence length
$l_{\mathrm{c}}$ and magnetic field strength $B$ for the observed blurring of
$20^{\circ}$ at $R=10$ EV. of $\sim 20^{\circ}$ for a particle with a rigidity
of 10 eV. For this model, a test statistic of
TS${}_{\mathrm{tot}}^{\mathrm{SBG}}=27.6$ is reached. Compared to the analysis
using only the (energy-independent) arrival directions [2, 6] with
TS${}_{\text{ADs-only}}^{\mathrm{SBG}}=25.0$, this corresponds to an increase.
Even when considering the most important systematic effects on the energy and
$X_{\mathrm{max}}$ scales, the total test statistic of the combined fit is
still TS${}_{\mathrm{tot}}^{\mathrm{SBG+syst}}=25.6$. Additionally, different
to Refs. [2, 6] no scan of the energy threshold is necessary as the whole
energy dependency is included in the astrophysical model for the combined fit.
In Ref. [16] the conversion of the total test statistic to a significance was
investigated. We found that the significance can approximately be calculated
using a $\chi^{2}$ distribution with two degrees of freedom for the two
additional fit parameters ($f_{0}$, $\delta_{0}$) compared to the respective
reference model. With this, the 1-sided significance of the SBG model amounts
to $\sim 4.5\sigma$ (or $\sim 4.7\sigma$ when not including the effects of
experimental systematic uncertainties). The total test statistic is dominated
by a region in which two interesting source candidates reside, NGC 4945 and
Centaurus A. In a separate investigation of a model where we picked just
Centaurus A on top of a homogeneous background, a test statistic of around
TS${}_{\mathrm{tot}}^{\mathrm{CenA}}\simeq 20$ is found, which is almost
independent of the background evolution and systematic uncertainties. The
contribution of Centaurus A to the total flux is estimated to be around $3\%$
at 40 EeV (although this fraction could be as large as $25\%$ for the case of
flat evolution), with a dominant nitrogen component reaching Earth above this
energy. The signal fraction and blurring for the model with Centaurus A as a
single source are in agreement with the ones obtained for NGC 4945 in the case
of the SBG model.

According to the results of this work, the $\gamma$-AGN model which exhibits a
test statistic of TS${}_{\text{ADs-only}}^{\gamma\mathrm{-AGN}}=17.9$ in the
arrival direction correlation analysis [2] does not agree with the energy-
dependent arrival directions. This is mostly due to the strong contribution of
the faraway blazar Markarian 421 which is not found in the data. Nevertheless,
as described above, the nearby AGN Centaurus A, which is also part of the
jetted AGN catalog, does agree well with the data. The $\gamma$-AGN model also
does reproduce the energy spectrum adequately, mostly due to the
overproduction of secondaries at energies of 5 EeV to 10 EeV owing to the
strong source evolution considered, something which should also affect other
potential scenarios associated with high-luminosity AGNs. Note also that the
$\gamma$-ray emission, which tends to favor the narrowly beamed jetted blazars
pointing towards Earth, may not be an optimal tracer of UHECR emission for
which the beaming is expected to be much broader due to deflections by cosmic
magnetic fields [77].

Future directions for the present analysis could be the inclusion of
deflections by the regular part of the magnetic field of the Galaxy as
described above. Also, other source catalogs with different flux weighting, as
well as an extension of the analysis to lower energies, where other
populations of sources may also contribute, could be investigated.
Additionally, the analysis will profit from the on-going upgrade of the Pierre
Auger Observatory, AugerPrime [78], which will enhance the sensitivity of the
surface detector to the masses of individual cosmic rays, allowing for better
predictions of their rigidities and possible magnetic field deflections.

##  Appendix A Results for the reference case of only homogeneously
distributed sources

In table 4, the fit parameters for the reference models of only homogeneous
background sources are given. These serve as a basis and comparison for the
results of the models including catalog sources discussed in section 4. The
fit parameters are in agreement with the ones found in other recent works of
the Pierre Auger Collaboration [10]. Note however, that Ref. [10] uses SimProp
instead of CRPropa3 and includes a local overdensity as well as a contribution
by an additional lower-energy component, so no exact agreement of the fit
parameters is expected.

|  Reference, $\boldsymbol{m=0}$ (flat) |  Reference, $\boldsymbol{m=3.4}$ (SFR) |  Reference, $\boldsymbol{m=5.0}$ (AGN-like)
---|---|---|---
| posterior | MLE | posterior | MLE | posterior | MLE
$\gamma$ | $-2.18_{-0.33}^{+0.31}$ | $-2.25$ | $-3.23_{-0.21}^{+0.24}$ | $-3.09$ | $-3.92_{-0.08}^{+0.01}$ | $-4.0$
$\log_{10}(R_{\mathrm{cut}}$/V) |  ++$18.18_{-0.03}^{+0.03}$ |  ++18.17 |  ++$18.09_{-0.02}^{+0.02}$ |  ++18.10 |  ++$18.04_{-0.01}^{+0.01}$ |  ++18.04
$I_{\mathrm{H}}$ | $6.5_{-1.6}^{+2.0}\times 10^{-2}$ | $7.0\times 10^{-2}$ | $9.0_{-9.0}^{+2.2}\times 10^{-3}$ | $<10^{-10}$ | $1.8_{-1.8}^{+0.3}\times 10^{-3}$ | $<10^{-10}$
$I_{\mathrm{He}}$ | $2.0_{-0.3}^{+0.3}\times 10^{-1}$ | $2.0\times 10^{-1}$ | $1.3_{-0.2}^{+0.2}\times 10^{-1}$ | $1.4\times 10^{-1}$ | $4.5_{-1.8}^{+1.8}\times 10^{-2}$ | $<10^{-10}$
$I_{\mathrm{N}}$ | $6.0_{-0.3}^{+0.4}\times 10^{-1}$ | $6.0\times 10^{-1}$ | $7.2_{-0.3}^{+0.3}\times 10^{-1}$ | $7.1\times 10^{-1}$ | $8.1_{-0.3}^{+0.3}\times 10^{-1}$ | $8.8\times 10^{-1}$
$I_{\mathrm{Si}}$ | $1.2_{-0.3}^{+0.2}\times 10^{-1}$ | $1.2\times 10^{-1}$ | $1.2_{-0.2}^{+0.3}\times 10^{-1}$ | $1.3\times 10^{-1}$ | $1.0_{-0.3}^{+0.3}\times 10^{-1}$ | $8.2\times 10^{-2}$
$I_{\mathrm{Fe}}$ | $2.1_{-0.9}^{+0.6}\times 10^{-2}$ | $1.9\times 10^{-2}$ | $2.6_{-0.9}^{+0.7}\times 10^{-2}$ | $2.3\times 10^{-2}$ | $3.8_{-1.1}^{+1.1}\times 10^{-2}$ | $3.6\times 10^{-2}$
$\boldsymbol{\log b}$ | $-270.6\pm 0.2$ |  | $-277.7\pm 0.3$ |  | $-297.0\pm 1.3$ |
$\boldsymbol{D_{E}}\ (N_{J}=14)$ |  | 22.2 |  | 28.1 |  | 63.2
$\boldsymbol{D_{X_{\mathrm{max}}}}\ (N_{{X}_{\mathrm{max}}}=74)$ |  | 126.8 |  | 132.4 |  | 134.1
$\boldsymbol{D}$ |  | 149.0 |  | 160.5 |  | 197.3
log $\boldsymbol{\mathcal{L}_{\mathrm{ADs}}}$ |  | $0$ |  | $0$ |  | $0$
log $\boldsymbol{\mathcal{L}}$ |  | $-250.5$ |  | $-256.2$ |  | $-274.6$

Table 4: Results for the reference models of only homogeneously distributed
sources.

|  Reference, $\boldsymbol{m=0}$ (flat) |  Reference, $\boldsymbol{m=3.4}$ (SFR) |  Reference, $\boldsymbol{m=5.0}$ (AGN-like)
---|---|---|---
| posterior | MLE | posterior | MLE | posterior | MLE
$\gamma$ | $-1.01_{-0.33}^{+0.41}$ | $-1.13$ | $-1.34_{-0.50}^{+0.63}$ | $-1.39$ | $-1.19_{-0.34}^{+0.41}$ | $-1.40$
$\log_{10}(R_{\mathrm{cut}}$/V) | $18.19_{-0.05}^{+0.03}$ | 18.19 | $18.20_{-0.06}^{+0.05}$ | 18.19 | $18.25_{-0.05}^{+0.04}$ | 18.25
$I_{\mathrm{H}}$ | $4.8_{-3.2}^{+2.3}\times 10^{-2}$ | $5.1\times 10^{-2}$ | $1.1_{-1.1}^{+0.1}\times 10^{-2}$ | $<10^{-10}$ | $1.1_{-1.1}^{+0.2}\times 10^{-2}$ | $3.5\times 10^{-8}$
$I_{\mathrm{He}}$ | $2.9_{-0.4}^{+0.4}\times 10^{-1}$ | $2.9\times 10^{-1}$ | $1.1_{-0.3}^{+0.4}\times 10^{-1}$ | $1.3\times 10^{-1}$ | $3.5_{-3.5}^{+0.8}\times 10^{-2}$ | $5.7\times 10^{-4}$
$I_{\mathrm{N}}$ | $5.3_{-0.4}^{+0.4}\times 10^{-1}$ | $5.2\times 10^{-1}$ | $6.4_{-0.7}^{+0.7}\times 10^{-1}$ | $6.4\times 10^{-1}$ | $5.9_{-0.5}^{+0.7}\times 10^{-1}$ | $6.8\times 10^{-1}$
$I_{\mathrm{Si}}$ | $1.2_{-0.4}^{+0.4}\times 10^{-1}$ | $1.2\times 10^{-1}$ | $1.8_{-0.5}^{+0.6}\times 10^{-1}$ | $1.8\times 10^{-1}$ | $2.8_{-0.5}^{+0.5}\times 10^{-1}$ | $2.4\times 10^{-1}$
$I_{\mathrm{Fe}}$ | $2.2_{-1.1}^{+0.9}\times 10^{-2}$ | $1.7\times 10^{-2}$ | $5.7_{-2.2}^{+1.4}\times 10^{-2}$ | $4.8\times 10^{-2}$ | $9.3_{-2.2}^{+1.9}\times 10^{-2}$ | $8.3\times 10^{-2}$
$\nu_{E}/\sigma$ | $-1.43_{-0.67}^{+0.57}$ | $-1.37$ | $0.19_{-0.60}^{+0.39}$ | $0.18$ | $\phantom{+}1.21_{-0.30}^{+0.29}$ | $1.17$
$\nu_{\mathrm{Xmax}}/\sigma$ | $-0.88_{-0.28}^{+0.33}$ | $-0.78$ | $-1.62_{-0.36}^{+0.34}$ | $-1.53$ | $-2.23_{-0.24}^{+0.27}$ | $-2.00$
$\boldsymbol{\log b}$ | $-257.9\pm 0.1$ |  | $-267.0\pm 0.1$ |  | $-280.3\pm 1.6$ |
$\boldsymbol{D_{\mathrm{syst}}}$ |  | 2.5 |  | 2.4 |  | 5.4
$\boldsymbol{D_{E}}\ (N_{J}=14)$ |  | 12.2 |  | 20.8 |  | 21.8
$\boldsymbol{D_{X_{\mathrm{max}}}}\ (N_{{X}_{\mathrm{max}}}=74)$ |  | 107.6 |  | 114.6 |  | 123.4
$\boldsymbol{D}$ |  | 122.3 |  | 137.8 |  | 150.6
log $\boldsymbol{\mathcal{L}_{\mathrm{ADs}}}$ |  | $0$ |  | $0$ |  | $0$
log $\boldsymbol{\mathcal{L}}$ |  | $-237.2$ |  | $-244.9$ |  | $-251.3$

Table 5: Results for the reference models of only homogeneously distributed
sources including experimental systematic uncertainties as nuisance
parameters.

## 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; 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 Kultur und
Wissenschaft 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, grant no. PN19150201/16N/2019 and
PN1906010 within the National Nucleus Program, and projects number TE128, PN-
III-P1-1.1-TE-2021-0924/TE57/2022 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.

## References

  * [1] The Pierre Auger Collaboration, A. Aab et al., _The Pierre Auger Cosmic Ray Observatory_ , _Nucl. Instrum. Methods Phys. Res. A_ 798 (2015) 172.
  * [2] The Pierre Auger Collaboration, P. Abreu et al., _Arrival Directions of Cosmic Rays above 32 EeV from Phase One of the Pierre Auger Observatory_ , _ApJ_ 935 (2022) 170.
  * [3] The Pierre Auger Collaboration, A. Aab et al., _Observation of a Large-scale Anisotropy in the Arrival Directions of Cosmic Rays above ${}^{\textrm{18}}$ eV_, _Science_ 357 (2017) 1266.
  * [4] The Pierre Auger Collaboration, A. Aab et al., _Large-scale cosmic-ray anisotropies above 4 EeV measured by the Pierre Auger Observatory_ , _ApJ_ 868 (2018) 4.
  * [5] The Pierre Auger Collaboration, P. Abreu et al., _Update on the correlation of the highest energy cosmic rays with nearby extragalactic matter_ , _ApJ_ 34 (2010) 314.
  * [6] The Pierre Auger Collaboration, A. Aab et al., _An Indication of anisotropy in arrival directions of ultra-high-energy cosmic rays through comparison to the flux pattern of extragalactic gamma-ray sources_ , _ApJ_ 853 (2018) L29.
  * [7] A. di Matteo et al., _UHECR arrival directions in the latest data from the original Auger and TA surface detectors and nearby galaxies_ , in _Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021)_ , vol. 395, p. 308, July, 2021, DOI.
  * [8] The Pierre Auger Collaboration, A. Aab et al., _Combined fit of spectrum and composition data as measured by the Pierre Auger Observatory_ , _JCAP_ 2017 (2017) 038.
  * [9] E. Guido for the Pierre Auger Collaboration, _Combined fit of the energy spectrum and mass composition across the ankle with the data measured at the Pierre Auger Observatory_ , in _Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021)_ , vol. 395, p. 311, 2021, DOI.
  * [10] The Pierre Auger Collaboration, A. Abdul Halim et al., _Constraining the sources of ultra-high-energy cosmic rays across and above the ankle with the spectrum and composition data measured at the Pierre Auger Observatory_ , _JCAP_ 2023 (2023) 024.
  * [11] B. Eichmann, J.P. Rachen, L. Merten, A. van Vliet and J. Becker Tjus, _Ultra-high-energy cosmic rays from radio galaxies_ , _JCAP_ 2018 (2018) 036.
  * [12] C. Ding, N. Globus and G.R. Farrar, _The Imprint of Large-scale Structure on the Ultrahigh-energy Cosmic-Ray Sky_ , _ApJL_ 913 (2021) L13.
  * [13] A. van Vliet, A. Palladino, A. Taylor and W. Winter, _Extragalactic magnetic field constraints from ultrahigh-energy cosmic rays from local galaxies_ , _MNRAS_ 510 (2022) 1289.
  * [14] D. Allard, J. Aublin, B. Baret and E. Parizot, _What can be learnt from UHECR anisotropies observations - I. Large-scale anisotropies and composition features_ , _A &A_ 664 (2022) A120.
  * [15] B. Eichmann, M. Kachelrieß and F. Oikonomou, _Explaining the UHECR spectrum, composition and large-scale anisotropies with radio galaxies_ , _JCAP_ 07 (2022) 006.
  * [16] T. Bister for the Pierre Auger Collaboration, _A combined fit of energy spectrum, shower depth distribution and arrival directions to constrain astrophysical models of UHECR sources_ , in _Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021)_ , vol. 395, p. 368, 2021, DOI.
  * [17] T. Bister for the Pierre Auger Collaboration, _Sensitivity of the combined fit of energy spectrum, shower depth distributions, and arrival directions at the Pierre Auger Observatory_ , _EPJ Web of Conferences_ 283 (2023) 03008.
  * [18] J. Biteau for the Pierre Auger Collaboration, _The ultra-high-energy cosmic-ray sky above 32 EeV viewed from the Pierre Auger Observatory_ , in _Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021)_ , vol. 395, p. 307, 2021, DOI.
  * [19] M. Ajello et al (The Fermi-LAT Collaboration), _3FHL: The Third Catalog of Hard Fermi-LAT Sources_ , _ApJS_ 232 (2017) 18.
  * [20] M. Ackermann et al. (The Fermi LAT Collaboration), _GeV Observations of Star-forming Galaxies with the Fermi Large Area Telescope_ , _ApJ_ 755 (2012) 164.
  * [21] C. Lunardini, G.S. Vance, K.L. Emig and R.A. Windhorst, _Are starburst galaxies a common source of high energy neutrinos and cosmic rays?_ , _JCAP_ 2019 (2019) 073.
  * [22] G.E. Romero, J.A. Combi, S.E. Perez Bergliaffa and L.A. Anchordoqui, _Centaurus A as a source of extragalactic cosmic rays with arrival energies well beyond the GZK cutoff_ , _Astropart. Phys._ 5 (1996) 279.
  * [23] J.H. Matthews, A.R. Bell, K.M. Blundell and A.T. Araudo, _Fornax A, Centaurus A, and other radio galaxies as sources of ultrahigh energy cosmic rays_ , _MNRAS_ 479 (2018) L76.
  * [24] The Pierre Auger Collaboration, J. Abraham et al., _Correlation of the Highest-Energy Cosmic Rays with Nearby Extragalactic Objects_ , _Science_ 318 (2007) 938.
  * [25] The Pierre Auger Collaboration, A. Aab et al., _Searches for Anisotropies in the Arrival Directions of the Highest Energy Cosmic Rays Detected by the Pierre Auger Observatory_ , _ApJ_ 804 (2015) 15.
  * [26] A.M. Hopkins and J.F. Beacom, _On the normalisation of the cosmic star formation history_ , _ApJ_ 651 (2006) 142.
  * [27] H. Yuksel and M.D. Kistler, _Enhanced Cosmological GRB Rates and Implications for Cosmogenic Neutrinos_ , _PRD_ 75 (2007) 083004.
  * [28] R. Alves Batista, R.M. de Almeida, B. Lago and K. Kotera, _Cosmogenic photon and neutrino fluxes in the Auger era_ , _JCAP_ 2019 (2019) 002.
  * [29] J. Heinze, D. Boncioli, M. Bustamante and W. Winter, _Cosmogenic Neutrinos Challenge the Cosmic Ray Proton Dip Model_ , _ApJ_ 825 (2016) 122.
  * [30] G.B. Gelmini, O. Kalashev and D.V. Semikoz, _Gamma-Ray Constraints on Maximum Cosmogenic Neutrino Fluxes and UHECR Source Evolution Models_ , _JCAP_ 2012 (2012) 044.
  * [31] M. Ahlers, L.A. Anchordoqui and S. Sarkar, _Neutrino diagnostics of ultra-high energy cosmic ray protons_ , _PRD_ 79 (2009) 083009.
  * [32] B. Peters, _Primary cosmic radiation and extensive air showers_ , _Il Nuovo Cimento (1955-1965)_ 22 (1961) 800.
  * [33] D. Allard and R.J. Protheroe, _Interactions of UHE cosmic ray nuclei with radiation during acceleration: consequences for the spectrum and composition_ , _A &A_ 502 (2009) 803.
  * [34] N. Globus, D. Allard, R. Mochkovitch and E. Parizot, _UHECR acceleration at GRB internal shocks_ , _MNRAS_ 451 (2015) 751.
  * [35] M. Unger, G.R. Farrar and L.A. Anchordoqui, _Origin of the ankle in the ultrahigh energy cosmic ray spectrum, and of the extragalactic protons below it_ , _PRD_ 92 (2015) 123001.
  * [36] K. Fang and K. Murase, _Linking high-energy cosmic particles by black-hole jets embedded in large-scale structures_ , _Nature Physics_ 14 (2018) 396.
  * [37] X. Rodrigues, J. Heinze, A. Palladino, A. van Vliet and W. Winter, _Active Galactic Nuclei Jets as the Origin of Ultrahigh-Energy Cosmic Rays and Perspectives for the Detection of Astrophysical Source Neutrinos at EeV Energies_ , _PRL_ 126 (2021) 191101.
  * [38] A. Condorelli, D. Boncioli, E. Peretti and S. Petrera, _Testing hadronic and photohadronic interactions as responsible for ultrahigh energy cosmic rays and neutrino fluxes from starburst galaxies_ , _PRD_ 107 (2023) 083009.
  * [39] M.S. Muzio, G.R. Farrar and M. Unger, _Probing the environments surrounding ultrahigh energy cosmic ray accelerators and their implications for astrophysical neutrinos_ , _PRD_ 105 (2022) 023002.
  * [40] R. Alves Batista et al., _CRPropa 3—a public astrophysical simulation framework for propagating extraterrestrial ultra-high energy particles_ , _JCAP_ 2016 (2016) 038.
  * [41] R. Beck, _Galactic and Extragalactic Magnetic Fields_ , _AIP Conference Proceedings_ 1085 (2008) 83.
  * [42] R. Durrer and A. Neronov, _Cosmological Magnetic Fields: Their Generation, Evolution and Observation_ , _Astron. Astrophys. Rev._ 21 (2013) 62.
  * [43] R. Alves Batista and G. Sigl, _Diffusion of cosmic rays at EeV energies in inhomogeneous extragalactic magnetic fields_ , _JCAP_ 2014 (2014) 031.
  * [44] R. Aloisio and V. Berezinsky, _Diffusive Propagation of Ultra-High-Energy Cosmic Rays and the Propagation Theorem_ , _ApJ_ 612 (2004) 900.
  * [45] S. Mollerach and E. Roulet, _Magnetic diffusion effects on the ultra-high energy cosmic ray spectrum and composition_ , _JCAP_ 2013 (2013) 013.
  * [46] S. Mollerach and E. Roulet, _Ultrahigh energy cosmic rays from a nearby extragalactic source in the diffusive regime_ , _PRD_ 99 (2019) 103010.
  * [47] M. Lemoine, _Extra-galactic magnetic fields and the second knee in the cosmic-ray spectrum_ , _PRD_ 71 (2005) 083007.
  * [48] D. Wittkowski for the Pierre Auger Collaboration, _Reconstructed properties of the sources of UHECR and their dependence on the extragalactic magnetic field_ , in _Proceedings of 35th International Cosmic Ray Conference — PoS(ICRC2017)_ , vol. 301, p. 563, Aug., 2018, DOI.
  * [49] A.J. Koning, S. Hilaire and M.C. Duijvestijn, _TALYS: Comprehensive Nuclear Reaction Modeling_ , _AIP Conference Proceedings_ 769 (2005) 1154.
  * [50] A.J. Koning and D. Rochman, _Modern Nuclear Data Evaluation with the TALYS Code System_ , _Nuclear Data Sheets_ 113 (2012) 2841.
  * [51] R.C. Gilmore, R.S. Somerville, J.R. Primack and A. Domínguez, _Semi-analytic modeling of the EBL and consequences for extragalactic gamma-ray spectra_ , _MNRAS_ 422 (2012) 3189.
  * [52] The Pierre Auger Collaboration, A. Aab et al., _Measurement of the cosmic-ray energy spectrum above 2.5 x 10 ${}^{\textrm{18}}$ eV using the Pierre Auger Observatory_, _PRD_ 102 (2020) 062005.
  * [53] M. de Domenico, M. Settimo, S. Riggi and E. Bertin, _Reinterpreting the development of extensive air showers initiated by nuclei and photons_ , _JCAP_ 07 (2013) .
  * [54] T. Pierog, I. Karpenko, J.M. Katzy, E. Yatsenko and K. Werner, _EPOS LHC: test of collective hadronization with LHC data_ , _PRD_ 92 (2013) 034906.
  * [55] The Pierre Auger Collaboration, A. Aab et al., _Depth of maximum of air-shower profiles at the Pierre Auger Observatory. II. Composition implications_ , _PRD_ 90 (2014) 122006.
  * [56] A. Zonca et al., _healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in Python_ , _J. Open Source Softw._ 4 (2019) 1298.
  * [57] K.M. Górski et al., _HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere_ , _ApJ_ 622 (2005) 759.
  * [58] S. Hackstein, F. Vazza, M. Brüggen, J.G. Sorce and S. Gottlöber, _Simulations of ultra-high energy cosmic rays in the local Universe and the origin of cosmic magnetic fields_ , _MNRAS_ 475 (2018) 2519.
  * [59] R. Jansson and G.R. Farrar, _The Galactic Magnetic Field_ , _ApJ_ 761 (2012) L11.
  * [60] M.S. Pshirkov, P.G. Tinyakov, P.P. Kronberg and K.J. Newton-McGee, _Deriving global structure of the Galactic Magnetic Field from Faraday Rotation Measures of extragalactic sources_ , _ApJ_ 738 (2011) 192.
  * [61] X.H. Sun, W. Reich, A. Waelkens and T. Enßlin, _Radio observational constraints on Galactic 3D-emission models_ , _A &A_ 477 (2008) 573.
  * [62] T.R. Jaffe et al., _Comparing polarized synchrotron and thermal dust emission in the Galactic plane_ , _MNRAS_ 431 (2013) 683.
  * [63] M. Erdmann, G. Müller, M. Urban and M. Wirtz, _The nuclear window to the extragalactic universe_ , _ApJ_ 85 (2016) 54.
  * [64] M. Unger and G.R. Farrar, _Uncertainties in the Magnetic Field of the Milky Way_ , in _Proceedings of 35th International Cosmic Ray Conference — PoS(ICRC2017)_ , vol. 301, p. 558, Aug., 2018, DOI.
  * [65] R. Fisher, _Dispersion on a Sphere_ , _Proc. R. Soc. Lond. A._ 217 (1953) 295.
  * [66] P. Virtanen et al., _SciPy 1.0: fundamental algorithms for scientific computing in Python_ , _Nature Methods_ 17 (2020) 261.
  * [67] N. Pullen and R.J. Morris, _Bayesian Model Comparison and Parameter Inference in Systems Biology Using Nested Sampling_ , _PLoS ONE_ 9 (2014) e88419.
  * [68] R.E. Kass and A.E. Raftery, _Bayes Factors_ , _JASA_ 90 (1995) 773.
  * [69] J. Salvatier, T.V. Wiecki and C. Fonnesbeck, _Probabilistic programming in Python using PyMC3_ , _PeerJ Comput. Sci._ 2 (2016) e55.
  * [70] D. Bergman for the Telescope Array Collaboration, _Telescope Array Combined Fit to Cosmic Ray Spectrum and Composition_ , in _Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021)_ , vol. 395, p. 338, July, 2021, DOI.
  * [71] R. de Almeida for the Pierre Auger Collaboration, _Large-scale and multipolar anisotropies of cosmic rays detected at the Pierre Auger Observatory with energies above 4 EeV_ , in _Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021)_ , vol. 395, p. 335, 2021, DOI.
  * [72] J.K. Lindsey, _Parametric Statistical Inference_ , Oxford University Press, Oxford, New York (May, 1996).
  * [73] S. Onn and I. Weissman, _Generating uniform random vectors over a simplex with implications to the volume of a certain polytope and to multivariate extremes_ , _Ann. Oper. Res._ 189 (2011) 331.
  * [74] A. Yushkov for the Pierre Auger Collaboration, _Mass Composition of Cosmic Rays with Energies above 10 ${}^{\textrm{17.2}}$ eV from the Hybrid Data of the Pierre Auger Observatory_, in _Proceedings of 36th International Cosmic Ray Conference — PoS(ICRC2019)_ , vol. 358, p. 482, 2019, DOI.
  * [75] J. Vícha for the Pierre Auger Collaboration, _Testing Model Predictions of Depth of Air-Shower Maximum and Signals in Surface Detectors using Hybrid Data of the Pierre Auger Observatory_ , _EPJ Web of Conferences_ 283 (2023) 02012.
  * [76] S. Mollerach and E. Roulet, _Progress in high-energy cosmic ray physics_ , _Progress in Particle and Nuclear Physics_ 98 (2018) 85.
  * [77] K. Murase, C.D. Dermer, H. Takami and G. Migliori, _Blazars as ultra-high-energy cosmic-ray sources: implications for TeV gamma-ray observations_ , _ApJ_ 749 (2012) 63.
  * [78] The Pierre Auger Collaboration, A. Aab et al., _The Pierre Auger Observatory Upgrade - Preliminary Design Report_ , 2016, DOI.