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

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. The model takes into account a rigidity-dependent magnetic field blurring and an energy-dependent evolution of the catalog contribution shaped by interactions during propagation. 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° 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σ (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 γ-ray emission are, however, disfavored as they cannot adequately describe the measured arrival directions.

at the Pierre Auger Observatory.The model takes into account a rigidity-dependent magnetic field blurring and an energy-dependent evolution of the catalog contribution shaped by interactions during propagation.
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 • 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σ (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 γ-ray emission are, however, disfavored as they cannot adequately describe the measured arrival directions.A Results for the reference case of only homogeneously distributed sources 27

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 km 2 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 EeV with a significance exceeding 5σ.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 ∼ 40 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σ (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 [12,13,15].
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.

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.

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σ (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 γ-ray flux (γ-AGNs) with a significance of 3.2σ.
For the SBGs, the γ-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 γ-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 database1 [18].The sources in the jetted γ-AGN catalog are selected from the 3FHL catalog [19].The selection contains only sources with a γ-ray flux measured with Fermi-LAT > 3.3 × 10 −11 cm −2 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 γ-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 γ-AGN catalog, the dominating source candidate is the powerful blazar Markarian 421 at a large distance of ∼ 134 Mpc.Also, the closest radio galaxy Centaurus A with a distance of ∼ 3.5 Mpc is part of the γ-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 [2,5,24,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 ≈ 2.2 (see above).The evolution of the emissivity of the background sources as a function of redshift z is modeled as (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 ∼ 1 is negligibly small [8,28], 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 ( 1 H), helium ( 4 He), nitrogen ( 14 N), silicon ( 28 Si), and iron ( 56 Fe).The generation rate at the sources Q inj (defined as the number of nuclei ejected per unit of energy, volume and time) as a function of the injected energy E inj and the mass number A inj of each representative element is modeled as a power-law with a broken exponential cutoff function which sets in when the maximum energy Here, γ is the spectral index, R 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 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 as2 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 , γ, R cut , I H , I He , I N , I Si , I Fe ) by a fit to the observed data, these describe the state of the spectrum and composition after leaving the source environment.

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 effect3 .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 [50,51] is used.For describing the EBL, the Gilmore model [52] 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.3km/s/Mpc, Ω m = 0.315, Ω Λ = 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 inj /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 ≈ 2.2).Simulations are done for the five representative elements described in section 2.1, so that in total this amounts to 5 × 118 × 150 = 88,500 simulations of 10 4 injected particles each.Upon detection, the arriving UHECRs are grouped into detected mass number bins A det ∈ (1, 2-4, 5-22, 23-38, >39) and detected energy bins E det ∈ (10 18.0 , 10 18.02 ..., 10 20.98 , 10 21.0 ) eV.In the following, the propagation database will be denoted as , 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 det /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 det ∈ (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.

Modeling the contributions from individual sources
The propagation database can now be used to map the injected flux from eq. ( 2 Calculation of the homogeneous background flux: From this, the background arrival histogram µ back for the homogeneous background sources can be calculated as 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 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 δ(d s − d l )): The histogram µ sig (E e det , A k 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 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 exp (E e det , A k det , C s ) is defined and then multiplied with the arrival histogram µ sig to get the exposure-weighted histogram μsig .The weights w 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: (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 det /eV) = 19.5 − 19.6 (≈ 40 EeV), corresponding to the running index e = 15.Accordingly, the two arrival histograms are weighted in the following way: The catalog contribution as a function of the detected energy bin can be calculated from the two summands of eq.(2.8) as follows: As the catalog sources are on average much closer than the homogeneous background, especially for a strong evolution, f S (f 0 , E e det ) usually rises with energy.

Simulated observables
The histogram µ tot (E e det , A k 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 µ tot (E e det , A k det ) have to be compared to the measured ones on Earth.In our case, these observables are the energy spectrum, shower maxima (X 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.

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 ′ from [53] which transforms the energy spectrum predicted by the model to the reconstructed one including resolution and bias effects.The spectrum J(E e det ) is then calculated as with the total vertical (zenith angle θ ≤ 60 • ) exposure E θ≤60 • , the arrival histogram µ e,k from eq. (2.8), and the width of the energy bins ∆E e 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).

Depth of the shower maximum
The Gumbel distributions g( Xx max |E e det , A k det ) [54] are used to model the expected X 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 [55] hadronic interaction model.The X max bins denoted by Xx max are given as Xx max ∈ (550, 570, ..., 1050) g/cm −2 .The shower maximum depth distributions are influenced by the acceptance A and the resolution 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. [56].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: (2.11)For the comparison of the measured X max with the predicted model values, a histogram is produced: (2.12) Here, c normalizes the right side of the equation over all X x max bins to one.

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 [57] is used which is based on HEALPix4 [58].Considering nside = 64 it divides the sky into n pix = 49,152 pixels of the same angular size.The binning approximately corresponds to the detector resolution of 0.9 • [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 E θ≤80 • up to zenith angle θ = 80 • in each healpy pixel 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: .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 [59] 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.[60][61][62][63].These do, however, differ substantially regarding the predictions of the arrival directions [64,65].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 association 5 .
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 [66] F (x|µ, κ), the equivalent of a 2-dimensional Gaussian on the sphere.The concentration parameter κ describes the width of the distribution.It is linked with the blurring angle δ via κ = 1/δ 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 δ(E e det , A k det ) as with the free model parameter δ 0 corresponding to the spread of protons with energies of 10 EeV (or that of nitrogen with 70 EeV).Here, the detected rigidity 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 det and A det by using the representative atomic numbers Z det = 1, 2, 7, 14 and 26 for the detected mass number bins A det ∈ {1}, {2, 3, 4}, {5, ..., 22}, {23, ..., 38} and {39, ..., 56}, respectively.The von Mises-Fisher distributions are then stored in a histogram F (pix p |E e det , A k det , C s ) using again the same healpy pixels pix p as for the background.The final arrival probability map for the catalog sources is then calculated by matrix multiplication: with c e normalizing p S(E e det , 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 θ ≤ 60 • are used, see section 3.2): Finally, the probability densities for background and catalog sources are added using the energy-dependent signal fraction function f S (f 0 , E e det ) from eq. ( 2.9): (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 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, δ 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.

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 ν E/Xmax 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% [53].A systematic shift of the energy scale, parameterized by the nuisance parameter ν 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.[8,10], 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 highacceptance region.The X max scale uncertainty is influenced by effects like the calibration and the atmospheric conditions, leading to an energy-dependent uncertainty between ±6 g/cm 2 and ±10 g/cm 2 [56].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 max scale.The X max scale uncertainty is included in the astrophysical model by shifting the X max values used for the evaluation of the Gumbel distributions in eq. ( 2.11) by ν Xmax , 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 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 max , a negative shift of the X max scale (implying that the true composition is actually heavier than measured) is represented here by shifting the mean modeled X 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.

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.

Inference of model parameters
The astrophysical model contains several free parameters (γ, R cut , I A , f 0 , δ 0 , Q 0 ) ∈ Θ as well as the nuisance parameters ν E and ν Xmax .

Fitting techniques
The inference of the model parameters is performed with two methods.The first one is a gradient-based minimizer based on scipy [67].It determines the best fit according to the likelihood function L(Θ) (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(Θ|d, M), where d denotes the data and M the model.For the estimation of the posterior distributions two ingredients are needed, the likelihood function and the prior distributions p(Θ|M), according to Bayes theorem: Here, p(d|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 [68,69].For the sampling of the posterior distributions, we use the Sequential Monte Carlo sampler form PyMC [70] which can also estimate the Bayesian evidence.Note that the use of Bayesian methods for the identification of UHECR sources was explored before in Refs.[71][72][73] with promising outcomes.

Likelihood function
The total log-likelihood function is given as a sum of the single likelihood functions for each independent observable: where "ADs" is used for the arrival directions.The likelihood for the energy spectrum L E is given as a Poissonian [8] log Here, n e is the measured number of events per energy bin e, and µ e = e ′ R e,e ′ µ e ′ 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 det /eV) = 18.7 in previous works [8,74].In Ref. [10] it was shown that an additional sub-ankle component can still have an influence around log 10 (E det /eV) = 18.7, and that the extragalactic component only completely dominates the flux above ∼ 10 EeV.For that reason, the minimum energy is set to log 10 (E det /eV) = 19.0 in this work.Nevertheless, it is enforced that the flux predicted by the astrophysical model between log 10 (E det /eV) = 18.7 and log 10 (E det /eV) = 19.0does not exceed the measured one by using a one-sided Poissonian likelihood for these energy bins.
The shower maxima X max information is binned into the histogram µ ẽ,x following eq.(2.12).Due to the small statistics of the X max measurements with the FD at the highest energies, the bins above log 10 (E det /eV) = 19.6 are combined into one, consequently denoted by ẽ instead of e.Since the energy spectrum information is already incorporated by the energy likelihood and should not have an influence twice, the X max likelihood is given as a multinomial instead of a Poissonian For the X max distributions, the minimum energy is also set to log 10 (E det /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 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 (pdf e,p ) n e,p . (3.5) For the arrival directions, a minimum energy of log 10 (E det /eV) = 19.2 is used which is higher than that for energy and X 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 eV ≈ 16 EeV, this analysis is in line with the other recent arrival directions analyses we performed [2,75], and can use the same data set for better comparison.
As the experimental systematic uncertainties on the energy and X max scale, ν E and ν Xmax , are expressed in units of standard deviations, the nuisance parameters ν are expected to follow a Gaussian likelihood with mean 0 and standard deviation 1 When the experimental systematic uncertainties are considered in the fit, L syst has to be included in the total likelihood (eq.(3.4)).
For the Poissonian energy and the multinomial X max likelihood, a deviance D can be defined which characterizes the goodness of fit [76].It is given as twice the negative likelihood ratio between the fitted model and the saturated model L 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 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 max , and systematics, and the likelihood value will be given for the arrival directions.

Prior distributions
For the spectral index γ ∈ [−4, 3] and the logarithmic rigidity cutoff log 10 (R cut /V) ∈ [18.0, 20.5], flat priors are used as in previous work [8].The lower border for γ 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 ∈ [0, 50%].Before, we also tested the range f 0 ∈ [0, 100%], but found that models with very large signal fractions cannot recreate the findings of Refs.[2,6], so the maximum allowed signal fraction was accordingly set to 50%.The blurring prior is also flat and can have values δ 0 ∈ [0 • , 86 • ].The upper border value (1.5 rad) corresponds, for reasonable UHECR compositions (e.g.nitrogen at ≈ 40 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 [77].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 ∈ [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σ and 4σ is used to keep the nuisance parameters in a reasonable range.

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 • from Ref. [53] are used.The exposure of the data set is 60,426 km 2 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 [78].Above 10 19.0 eV, the data set contains 1022 events which are binned into a 2-dimensional histogram in energy and X 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 • ) and inclined events (zenith angle larger than 60 • and up to 80 • ).The exposure amounts to 95,700 km 2 sr yr for the vertical subset and 26,300 km 2 sr yr for the inclined subset.

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

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 interval6 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 γ and the maximum rigidity log 10 (R cut /V) is visible as in Ref. [8], as well as a correlation between the signal fraction f 0 and the blurring δ 0 as expected from Ref. [6].The integral composition fractions are typically almost uncorrelated with the other parameters.
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 ≳ 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 flat evolution is smaller than that for SFR by around ∆D ≈ 10.The maximum rigidity is relatively small with log 10 (R cut /V) ≈ 18.2 V.   4) has been set to 0 for better readability of the values of log b, log L, and log L ADs .
The spectral index is very hard already for flat evolution, but can become as hard as γ ≈ −3 for the SFR case, as also found in [10].The hard injection spectrum allows for a good description of the pronounced features of the energy spectrum, as well as the small X max variance.As described above, γ is not the spectral index of the spectrum at acceleration, but rather after leaving the source environment.Thus, in-source interactions or magnetic field confinement in the source could lead to hard best-fit values of γ even for shock acceleration scenarios for which a softer injection ∝ E −2 is expected [35,38,79].Other explanations could be a low-energy flux suppression by the EGMF [48,49], or the effect of systematic uncertainties (see below).

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 ≈ 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 ∼ 3% at 40 EeV to ∼ 10% at 100 EeV.The uncertainty on the signal fraction is large as described above, so contributions up to ∼ 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 ≃ 10 eV at 70 EeV.For this, the model predicts a blurring of around δ 0 = 20 • .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 ∼ 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 max scale uncertainty as a nuisance parameter in the model, as will be shown in sec.4.1.3.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.
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 Cen A, m = 0.0 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 ≈ 20%.In fig.2, it is again visible how the catalog contribution becomes relevant at higher energies, where it reaches values of up to ≈ 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 • (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 • and 30 • , the reduced χ 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.
Cen A, m = 0.0 Auger, 20 Cen A, m = 3.4 Auger, 20 Auger, 20  Observatory inside the 20 • 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 • 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 EeV ≃ (f source /0.05) × 10 39 erg/s, with f 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 × ω NGC4945 flux 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.

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 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σ 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 max scale.Here, a negative shift is preferred, which is around −1σ for m = 0 and −1.5σ for m = 3.4, in accordance with results for the reference models (see appendix A).This shifts the X max scale up by around 10 − 15 g/cm 2 (see section 2.5), which is in line with the findings from Ref. [80].According to these results, the true X 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 γ ≈ −1 for both evolutions.21.9The 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 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 max scale.

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: 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.As is apparent from the table, the arrival directions observable provides the largest contribution to the total test statistic.This is understandable, as the reference model already provides a proper fit of the energy spectrum and X max data [10], so the subdominant contribution by the nearby source candidates only has a minor impact on these observables.For the arrival directions, however, the improvement from fully isotropic arrival directions in the reference model, to the anisotropic ones provided by the model including source candidates (Fig. 4), is substantial.Hence, the energy spectrum and X max distributions are necessary for constraining the source emission, while the arrival directions are most important for the differentiation of different source candidates.This is as expected from a simulation study [17].But, from that 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.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 det /eV) = 19.3 and the bins around log 10 (E det /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≃ 10 19.6 eV and 60 EeV≃ 10 19.8 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 det /eV) = 19.3 the helium contribution has its peak, while at log 10 (E det /eV) = 19.6 nitrogen is predominant, and at log 10 (E det /eV) = 19.8silicon begins to emerge.At these energies, the rigidities correspond to R He = 10 19.3  2 V ≃ 10 EV, R N = 10 19.6 7 V ≃ 6 EV and R Si = 10 19.8 14 V ≃ 5 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 [2,6] is so large (TS SBG ADs-only ≃ 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 det /eV) = 19.3 could just be a statistical fluctuation.The sum over all energy bins corresponds to the values for TS ADs 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 det /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 ∼ 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 ∼ 4 to 5 to the arrival directions test statistic.

The γ-AGN model
When the γ-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 γ ≈ −3.5 is found in combination with a small maximum rigidity of log 10 (R cut /V) ≈ 18.1.The best-fit signal fraction is around f 0 ≈ 15%.The γ-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 γ-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 γ-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 γ−AGN ADs = 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 γ-AGN model does not depend on the source evolution of the background sources.
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 γ-AGN model and the measured data7 .For that, we include an additional fit parameter β e in the arrival directions modeling, and replace eq.(2.14) by the following

.2)
The prior for the distance-dependent blurring is flat β e ∈ [0 • , 11.5 • ].This means, a 10 EV particle can at most obtain a distance-dependent blurring of ∼ 20 • 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 γ-AGN model with distance-dependent blurring do not differ much from the ones without it.The posterior distribution of β e is broad between β e = 5 • and the maximum allowed β e = 11.5 • , while the distance-dependent blurring has decreased to δ 0 ≃ 20 • .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 det /eV) ≃ 19.5.These arrival directions now describe the data better than the isotropic reference model, as visible in fig.8. Nevertheless, with TS γ−AGN+EGMF ADs ≃ 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.[2,6] and found that it cannot reach large enough arrival direction test statistics.
For that reason, we can conclude that the γ-AGN model where the flux contribution scales linearly with the γ-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.

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, for the first time, 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 energydependent 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, at least in the directions of the strongest sources.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 become stronger as the contribution from the nearby catalog sources rises with the energy.This enables the modeling of the energy evolution of the arrival directions, which can then be compared to the data in energy bins, and not just above an energy threshold as in previous analyses [8,10].
The largest test statistic of TS SBG 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 CenA tot ≃ 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 γ-AGN model, which exhibits a test statistic of TS γ−AGN ADs-only = 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 due to its large UHECR flux weight, calculated from the γ-ray flux in the model.From this, it can be concluded that the γ-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 [82].The γ-AGN model also does not describe 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.Nevertheless, as described above, a model based solely on the nearby AGN Centaurus A, which is also part of the jetted AGN catalog, does reproduce the data well.
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 [83], 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.

A Results for the reference case of only homogeneously distributed sources
In table 4 and table 5, 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.Table 5. Results for the reference models of only homogeneously distributed sources including experimental systematic uncertainties as nuisance parameters.

Figure 1 .
Figure 1.Posterior distributions for the SBG model.The pink cross marks the MLE, the black vertical line the posterior mean.

Figure 2 .Figure 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 4 .
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 det /eV) =19.3 (upper row ), log 10 (E det /eV) =19.6 (middle row ) and log 10 (E det /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.

Figure 5 .
Figure 5. Modeled spectra on Earth inside a circular region of angular size 20 • (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• 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 • region have not been fitted explicitly (see section 3.1.2),so this serves as a cross check of the model.

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

Figure 8 .
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.

Figure 9 .
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 TS ADs given in table3.

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

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 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 L ref ADs for fully isotropic arrival directions (as in the reference models, see table

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

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 TS tot .*Notethat 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 γ-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 β e compared to the other models (see section 4.2).