Arrival Directions of Cosmic Rays above 32 EeV from Phase One of the Pierre Auger Observatory

Cosmic rays are observed up to the astounding energies of more than 10²⁰ eV, making them the most energetic particles known in the universe. However, the origin of these particles remains elusive. The search for the sources of ultra-high-energy cosmic rays (UHECRs), at energies above a few EeV (1 EeV ≡ 10¹⁸ eV), is challenging since they are almost all charged particles and thus deflected by the magnetic fields permeating the interstellar, intra-halo, and intergalactic media (see, e.g., Alves Batista et al. 2019, for an overview). These magnetic fields are difficult to study and their modeling is far from being complete. However, above a few tens of EeV, the deflections could be small enough for cosmic rays to retain some directional information on the position of their sources, at least for nuclei with a sufficiently small charge (e.g., Erdmann et al. 2016; Farrar & Sutherland 2019).

The cosmological volume within which UHECR sources should be sought is fortunately limited. Cosmic rays at EeV energies can interact with the photon backgrounds populating intergalactic space, through the so-called Greisen–Zatsepin–Kuzmin (GZK) effect (Greisen 1966; Zatsepin & Kuz'min 1966). In particular, protons are expected to undergo photopion production and nuclei photodissociation interactions. The mean free path for energy losses depends on the cosmic-ray mass and energy. At 100 EeV, the loss length is of the order of 200–300 Mpc for protons and iron and 3–6 Mpc for intermediate nuclei such as helium and nitrogen (Allard 2012; see also Figure 6 from Addazi et al. 2022 for a recent overview). Such short distances mean that the sources of the highest-energy cosmic rays must be in the local universe.

The recent detection by the Pierre Auger Collaboration of a dipolar anisotropy in the arrival directions of UHECRs with energies above 8 EeV is evidence that the majority of UHECR sources are not in the Milky Way (Pierre Auger Collaboration2017a). The direction of the dipole points ∼120° away from the Galactic center and is instead consistent at the 2σ level with the local distribution of stellar mass (2MASS Redshift Survey; Huchra et al. 2012), after accounting for the deflections expected in the Galactic magnetic field (Jansson & Farrar 2012). Even without relying on magnetic deflections, the case for a density of UHECR sources following local extragalactic structures is further strengthened by the consistency at the 1σ confidence level (C.L.) between the directions of the UHECR anti-dipole and of the Local Void at equatorial coordinates (α, δ) = (294°, 15°) or Galactic coordinates (l, b) = (51°, −3°) (Biteau 2021). Combined with the growth of the dipole amplitude with energy expected from the shrinking horizon out to which extragalactic sources remain visible (Pierre Auger Collaboration 2018a), the properties of the large-scale anisotropy discovered by the Pierre Auger Collaboration provide a growing body of evidence against a Galactic origin of these cosmic rays. Which (classes of) extragalactic sources host UHECR accelerators nonetheless remains an open question.

In this article, we update previous searches for anisotropies at the highest energies (Pierre Auger Collaboration 2015a, 2018b) with an unprecedentedly large data set. In particular, we exploit the entire Phase 1 of the Pierre Auger Observatory, i.e., the phase preceding the AugerPrime upgrade (Pierre Auger Collaboration2016a). Important progress has been made on estimating the mass distribution of UHECRs using only the surface detector of the Observatory with its full duty cycle (see, e.g., Pierre Auger Collaboration 2016b, 2017b; Ave et al. 2017; Pierre Auger Collaboration 2021a, 2021b). However, the proposed methods are still not ready to be employed in arrival-direction studies, e.g., by selecting only candidate light nuclei which would be less deflected by magnetic fields, should such a subsample exist in the data set. In the following, we then consider, as in previous works, only the energy and arrival direction of each event recorded with the Pierre Auger Observatory over 17 yr of operation.

The data set includes more than 2600 events with energies E ≥ 32 EeV and zenith angles up to 80°, as described in Section 2. The release of this data set complements the publication of the arrival directions of events at energies between 4 and 8 EeV and above 8 EeV made available by the Pierre Auger Collaboration (2017a). ¹⁰⁰ The choice of an energy threshold at 32 EeV for the present release anticipates upcoming publications focused on lower energy bins, namely 8–16 EeV and 16–32 EeV, as investigated by, e.g., the Pierre Auger Collaboration (2018a) and Pierre Auger Collaboration (2020a), where ∼1500 and ∼2000 events were studied above 32 EeV, respectively. In Section 3, we describe a first set of analyses that are not based on specific source models, i.e., a blind search for excesses in the sky, an autocorrelation study, and the search for correlations with the Galactic and supergalactic planes as well as the Galactic center. Section 4 is devoted to the comparison of UHECR arrival directions with the expected flux pattern from specific classes of galaxies traced by their electromagnetic emission, from radio wavelengths to gamma-rays. Finally, Section 5 is devoted to a more in-depth study of the Centaurus region, which has intrigued the UHECR community since the early days of the Pierre Auger Observatory (Pierre Auger Collaboration 2007).

To encourage further studies of the Phase 1 high-energy data set, this article is accompanied by supplementary materials. These include the data set itself in Appendix A and the dedicated analysis software in Appendix B. Appendix C describes the catalogs of galaxies used here.

The Pierre Auger Observatory (Pierre Auger Collaboration2015b) is located in Argentina near the town of Malargüe. Stable data acquisition began on 2004 January 1. The Observatory is composed of a surface detector (SD) made of 1660 water-Cherenkov stations distributed on a triangular grid overlooked by a fluorescence detector (FD). The FD consists of 27 telescopes at four locations on the perimeter of the SD array.

Here, we analyse the events with reconstructed energies larger than 32 EeV recorded with the SD array from 2004 January 1 to 2020 December 31. The SD is used to sample secondary particles in air showers and has full efficiency above 4 EeV with ∼100% duty cycle.

Events recorded with the SD are reconstructed differently based on their arrival direction in local coordinates: events with zenith angles, θ, less than 60° are called vertical events, while events arriving with zenith angles from 60° to 80° are called inclined events. Vertical events are included when the SD station with the largest signal is surrounded by at least four active stations. This a priori condition is complemented by the a posteriori requirement that the reconstructed core of the shower falls within an elementary isosceles triangle of active stations. These requirements ensure that the footprint of the shower is well-contained within the array, with ample data for an accurate reconstruction (Pierre Auger Collaboration 2010a). Inclined events, on the other hand, are selected if the station closest to the reconstructed core position is surrounded by at least five active stations. Note that other analyses performed by the Pierre Auger Collaboration at lower energies may use a tighter selection. For example, the UHECR spectrum of Pierre Auger Collaboration (2020b) is measured by requiring that all six active stations around the one with the highest signal are active. We are able to use a relaxed selection as the high-energy events included here all have large footprints, with an average of 17.7 triggered stations. We inspected each event and verified that the reconstruction was robust even with inactive stations in the core region. With respect to previous analyses, the identification of active stations that were not triggered has been improved to ensure a better selection. This was done through an a posteriori check of the consistency of the signal distribution at the ground: if a station is not triggered in a region of the array where the signal is more than twice that of the full trigger efficiency, which occurs for 11 events in the data set, the station is classified as non-active at the moment of the event (Pierre Auger Collaboration 2010a).

The selection results in 2040 events with θ < 60° and 595 with θ ≥ 60° above 32 EeV. ¹⁰¹ The exposure can be computed in a geometrical way since we are operating above the energy threshold for full efficiency for both data samples (3 EeV for vertical and 4 EeV for inclined). The geometrical exposure for the selection and time span considered is 95,700 km² sr yr for the vertical sample and 26,300 km² sr yr for the inclined data set.

The reconstruction procedure for vertical events is described in detail in Pierre Auger Collaboration (2020c). The arrival direction is determined by fitting a spherical model to the arrival times of particles comprising the shower front. For inclined events, the reconstruction procedure is described in Pierre Auger Collaboration (2014). The arrival direction is, in this case, obtained by fitting the arrival times with a front which takes into account muon propagation from its production point. For both data sets, the angular resolution, defined as the 68% containment radius, is better than 1° at all energies considered here.

The energy estimate is based on different observables for the two samples. The signal at a reference distance of 1000 m from the shower core, S(1000), is used for the vertical sample. The inclined reconstruction uses as an estimator N₁₉, which represents the muon content of the shower with respect to a reference simulated proton shower with energy E = 10¹⁹ eV. For both samples, a correction is applied to take into account the absorption that showers undergo at different zenith angles. This correction is performed through a data-driven procedure called constant intensity cut, which is described in Pierre Auger Collaboration (2020b). The constant intensity cut method is used to convert S(1000) and N₁₉ for each shower to the value they would have if the same shower had arrived from a reference zenith angle of 38° and 68° for vertical and inclined events, respectively. The corrected energy estimators, S₃₈ and N₆₈, are then calibrated using hybrid events, i.e., events observed with both the FD and the SD. Since the FD analysis enables a quasi-calorimetric measurement of the shower energy, the calibration procedure results in a reliable energy estimation for the whole SD data set without using air-shower simulations. The systematic uncertainty in the energy calibration is ∼14% while the energy resolution for the SD at the energies considered here is ∼7% (Pierre Auger Collaboration 2014, 2020d).

We checked the consistency between the vertical and inclined data sets by comparing the ratio of number of events in the two samples, N_incl/N_vert = 0.292 ± 0.014, and the value expected from the ratio of geometrical exposures, accounting for the finite energy resolution of each data stream, $\tfrac{{\omega }_{\mathrm{incl}}/{c}_{\mathrm{incl}}(\geqslant 32\,\mathrm{EeV})}{{\omega }_{\mathrm{vert}}/{c}_{\mathrm{vert}}(\geqslant 32\,\mathrm{EeV})}=0.278$. In the latter ratio, ω is the geometrical exposure for each data set, which does not depend on energy, and c(≥32 EeV) accounts for the net spillover of events from low to higher energies (see the unfolding procedure described in Pierre Auger Collaboration 2020b). The ratios are in agreement at the 1σ C.L., showing that the vertical and inclined samples can be used together. To keep the analysis as data-driven as possible, we use the ratio of events observed above 32 EeV as the expected exposure ratio when constructing simulated data sets above any energy threshold. It should be noted that at the highest energies probed here, E ≥ 80 EeV, a deficit of inclined events is observed at a significance level of 2.5σ. A further discussion of this deficit, which does not affect the results presented below, is provided in Appendix A together with the information on how to access the data.

An earlier wide-ranging search with the Observatory for small- and intermediate-scale anisotropy was reported by Pierre Auger Collaboration (2015a). Searches for localized excesses in top-hat windows of angular radius Ψ across the entire field of view of the Observatory, or around the Galactic center, Centaurus A, and candidate host galaxies identified in multiwavelength surveys, were performed by comparing the expected and observed numbers of events within the window. Similar analyses were performed along the Galactic and supergalactic planes, by counting the number of events within an angle Ψ from these structures, and an autocorrelation study exploited the number of pairs of events separated by less than Ψ. The analyses were repeated above energy thresholds ranging from 40 to 80 EeV. An additional scan on the maximum distance of the sources was performed for analyses against catalogs of candidate host galaxies. Both scans in energy threshold and maximum distance were motivated by the limited horizon from which UHECRs can reach Earth, although the determination of its observational value remains hindered by uncertainties on UHECR composition.

In this Section, we update the results presented by Pierre Auger Collaboration (2015a), with the exception of the search for correlations with catalogs, which is performed in Section 4.

3.1. Search for Localized Excesses

The first analysis is a blind search for excesses over the fraction of the sky covered with the Observatory. The number of UHECRs detected in circular windows on the sky (N_obs) is compared to that expected, in the same window, from an isotropic distribution of events (${N}_{\exp }$). This search is performed over the entire field of view, which covers about 85% of the sky. The search windows are centered on a HEALPix grid (HEALPix v3.70; Górski et al. 2005), defined by the parameter nSide = 64, which sets the size of the pixels to be of the order of the angular resolution of the Observatory. Events are counted within search windows of radius Ψ, ranging from 1° to 30° in 1° steps. Similarly, the search is performed by selecting events above energy thresholds, E_th, ranging from 32 to 80 EeV in 1 EeV steps. For each window and energy threshold, we estimate the binomial probability of obtaining by chance N_obs or more events from an isotropic distribution of data. The computation of ${N}_{\exp }$ is performed by simulating events with coordinates distributed according to the sum of the vertical and inclined exposures, weighted in proportion to the observed number of events at energies above 32 EeV (see Section 2). For each realization of the simulated data set, the number of events is of the same size as observed across the field of view. Simulated events follow the same energy distribution as the observed events. Performing the analysis on simulated isotropic data sets allows us to take into account the trial factors for having tested different directions, radii and energy thresholds. We consider the post-trial probability as the fraction of these simulations with an equal or lower local p-value than the best one obtained with the observed data set.

We also compute the local Li–Ma significance (equation (17) in Li & Ma 1983) for each point in the sky, where the ON-region is centered on each point of the HEALPix grid and the OFF-region is defined as the remainder of the field of view. The local significance map is displayed in Galactic coordinates in Figure 1. The most significant excess, with 5.4σ local significance, is found above an energy threshold of 41 EeV within a top-hat window of 24° radius centered on equatorial coordinates (α, δ) = (1963, − 466), which corresponds to Galactic coordinates (l, b) = (3054, 162). At this position of the parameter space, 153 events are observed when 97.7 are expected from isotropy. The local p-value in this position is 3.7 × 10⁻⁸, resulting in a post-trial p-value of 3%.

Zoom In Zoom Out Reset image size

3.2. Autocorrelation

Another model-independent approach to assess the clustering of events is the search for autocorrelation, i.e., counting pairs of events separated by a given angular distance. This approach is particularly effective if the events form multiple clusters on similar angular scales in different directions in the sky.

Following Pierre Auger Collaboration (2015a), we count the number of event pairs, N_obs, above energy thresholds ranging from 32 to 80 EeV, that are separated by less than an angle Ψ ranging from 1° to 30° in steps of 025 up to 5° and of 1° above. We compute the expected number of pairs, ${N}_{\exp }$, by analysing simulated isotropic event sets of the same size as the observed data set. For each Ψ and E_th, we consider the local p-value as the fraction of simulated data sets, f(E_th, Ψ), for which ${N}_{\exp }\geqslant {N}_{\mathrm{obs}}$. The values of f are shown in Figure 2(a) and the best results are shown in Table 1.

Zoom In Zoom Out Reset image size

Table 1. Results of the Search for Autocorrelation and Correlation with Astrophysical Structures

Search	Eth [EeV]	Angle, Ψ [deg]	Nobs	${N}_{\exp }$	Local p-value, ${f}_{\min }$	Post-trial p-value
Autocorrelation	62	3.75	93	66.4	2.5 × 10−3	0.24
Supergalactic plane	44	20	394	349.1	1.8 × 10−3	0.13
Galactic plane	58	20	151	129.8	1.4 × 10−2	0.44
Galactic center	63	18	17	10.1	2.6 × 10−2	0.57

Note. The energy threshold, E_th, and the search angle, Ψ, which minimize the local p-value, are based on the number of observed and expected events/pairs. The post-trial p-value accounts for the scan in energy threshold and search angle, Ψ.

Download table as:  ASCIITypeset image

3.3. Correlations with Structures

In Pierre Auger Collaboration (2015a), we presented the results of cross-correlation studies with three flux-limited catalogs: the 2MASS Redshift Survey of near-infrared galaxies (Huchra et al. 2012), the Swift-BAT 70-month catalog of active galactic nuclei (AGNs) observed in hard X-rays (Baumgartner et al. 2013), and the catalog of radio-emitting galaxies from van Velzen et al. (2012). Such cross-correlation analyses inherently assume all galaxies under investigation to have an equal weight (standard-candle approach) and do not easily account for the inverse-square law of the UHECR flux, nor for its attenuation resulting from energy losses induced by propagation. These limitations were addressed by Pierre Auger Collaboration (2018b) through a likelihood-ratio test that expanded upon the maximum-likelihood test presented by Pierre Auger Collaboration (2010b). We also tested two additional catalogs based on gamma-ray observations from Fermi-LAT. The full-sky gamma-ray survey of Fermi-LAT has shown star-forming galaxies and jetted AGNs to be the main contributors to the extragalactic gamma-ray background at GeV energies, although their relative contributions remains uncertain (see, e.g., Ajello et al. 2015; Roth et al. 2021).

4.1. From Catalogs to UHECR Sky Models

4.2. Likelihood-ratio Analysis

As in Pierre Auger Collaboration (2018b), the correlation of UHECR arrival directions with the flux pattern expected from the catalogs is evaluated against isotropy using a likelihood-ratio analysis. The model as a function of direction u is computed in equal-area bins on the sphere using HEALPix v3.70 with the parameter nSide=64, as in Section 3.1.

The null hypothesis under investigation, H₀, is that of an isotropic flux distribution. Accounting for the directional exposure of the array, ω( u ), the isotropic model for the UHECR count density reads

$$\begin{eqnarray}&&{n}^{{H}_{0}}({\boldsymbol{u}})=\displaystyle \frac{\omega ({\boldsymbol{u}})}{{\sum }_{i}\omega ({{\boldsymbol{u}}}_{i})},\end{eqnarray} \tag{ 1 }$$

which is normalized so that the sum over the HEALPix pixels indexed over i and of direction u _i is equal to one.

The alternative hypothesis, H₁, in which H₀ is nested, is considered as the sum of the isotropic component and a component derived from the tested catalog. The amplitude of the latter component is the variable signal fraction, α. The isotropic remainder accounts for faint or distant galaxies not included in the catalogs or for a heavy nuclear component deflected away on large angular scales. The model for the UHECR count density under H₁ reads

$$\begin{eqnarray}&&{n}^{{H}_{1}}({\boldsymbol{u}})=(1-\alpha )\times {n}^{{H}_{0}}({\boldsymbol{u}})+\alpha \times \displaystyle \frac{{\sum }_{j}{s}_{j}({\boldsymbol{u}};{\rm{\Theta }})}{{\sum }_{i}{\sum }_{j}{s}_{j}({{\boldsymbol{u}}}_{{\boldsymbol{i}}};{\rm{\Theta }})},\end{eqnarray} \tag{ 2 }$$

where the index j runs over the galaxies in the catalog. The contribution to the UHECR flux from each galaxy, s_j ( u ; Θ), is modeled as a von Mises–Fisher distribution centered on the direction of the galaxy with smearing angle Θ. The amplitude of its contribution is proportional to the electromagnetic flux of the galaxy, ϕ_j , accounting for attenuation as a function of luminosity distance, a(d_j ), so that

$$\begin{eqnarray}&&{s}_{j}({\boldsymbol{u}};{\rm{\Theta }})=\omega ({\boldsymbol{u}})\times {\phi }_{j}a({d}_{j})\times \exp \left(\displaystyle \frac{{\boldsymbol{u}}\cdot {{\boldsymbol{u}}}_{j}}{2(1-\cos {\rm{\Theta }})}\right).\end{eqnarray} \tag{ 3 }$$

The von Mises–Fisher distribution is maximum in the direction of the galaxy of interest, u _j , effectively leaving aside coherent deflections which remain under-constrained by current models of the Galactic magnetic field (Erdmann et al. 2016). The smearing angle Θ, equivalent to a 2D Gaussian extent in the small-angle limit, is assumed to be the same for all galaxies in a given catalog. This parameter accounts for the average angular dispersion in intervening magnetic fields. As a note, normalization of the von Mises–Fisher distribution in Equation (3) is omitted, as it is the same for every galaxy and because the overall anisotropic component is normalized in the same way on the celestial sphere (see Equation (2)).

The likelihood-ratio test between the nested models H₀ and H₁ defines the test statistic (TS) as $\mathrm{TS}=2\mathrm{ln}({{ \mathcal L }}_{1}/{{ \mathcal L }}_{0})$, where the likelihood scores of the null and alternative hypothesis, ${{ \mathcal L }}_{0}$ and ${{ \mathcal L }}_{1}$, are obtained as the product over the events of the models ${n}^{{H}_{0}}$ and ${n}^{{H}_{1}}$, respectively. The evaluation of the TS is performed by grouping events by HEALPix bins. With an observed event count k_i in the direction u _i , the TS is evaluated as

$$\begin{eqnarray}&&\mathrm{TS}=2\sum _{i}{k}_{i}\times \mathrm{ln}\displaystyle \frac{{n}^{{H}_{1}}({{\boldsymbol{u}}}_{i})}{{n}^{{H}_{0}}({{\boldsymbol{u}}}_{i})}.\end{eqnarray} \tag{ 4 }$$

TS is maximized as a function of the two free parameters in the analysis (the search radius, Θ, and the signal fraction, α) above successive energy thresholds. Maximization can be achieved by scanning the 2D parameter space by steps of 0.2% in signal fraction and 0.2° in search radius. This approach provides an accurate estimate that is independent from any specific maximization algorithm. Alternatively, maximization with the Minuit package provides a fast estimate for simulated data sets, with an accuracy of the TS better than 0.1 units for event counts larger than 100. Above a fixed energy threshold, the TS is observed through Monte-Carlo simulations to follow a χ² distribution with two degrees of freedom under the null hypothesis (Wilks 1938). The 1 and 2σ C.L. on the best-fit parameters are set by isoTS contours differing from the maximum TS value by 2.3 and 6.2 units, respectively.

The scan in energy threshold is accounted for, as in Section 3, by estimating the post-trial p-value through isotropic Monte-Carlo simulations. The post-trial p-value, which accounts for the energy scan, differs from the local p-value expected from Wilks' theorem by a penalty factor that is well-approximated by a linear function of TS: pen = 1 + (0.30 ± 0.01) × TS. This empirical penalty factor is estimated from simulated isotropic data sets analyzed against each catalog and the uncertainty on the linear coefficient is estimated from the variance across the four tested catalogs. The penalty factor reaches a value of ∼ 10 for TS = 30.

4.3. Results

The radii and signal fractions maximizing the TS above fixed energy thresholds ranging in 32–80 EeV are displayed in Figure 3 for the four catalogs. The TS follows a double hump structure as a function of energy, with a first peak at energies above ∼40 EeV and a second peak at energies above ∼60 EeV. The latter peak corresponds to the maximum signal fraction for all catalogs, ranging in 11%–19%. Lower signal fractions ranging in 6%–16% are inferred from the global TS maximum, at energies above ∼40 EeV. As shown in the upper axis in Figure 3, the four times larger number of events in the first peak (1,387 above 40 EeV versus 331 above 60 EeV) yields a more significant deviation from isotropy above 40 EeV.

Zoom In Zoom Out Reset image size

The amplitude of variations of the best-fit parameters as a function of energy threshold can be evaluated against the statistical uncertainties on these parameters, as shown in Figure 4. As the search is performed above successive energy thresholds in steps of 1 EeV, successive energy bins have a non-negligible overlap. For reference, we estimate that there is a total of five to six independent energy bins, by identifying the successive reference energy thresholds above which the number of events is less that half that above a previous reference energy. Such a procedure suggests reference energy thresholds at E ≳ 32, 40, 50, 60, 70, and 80 EeV, with boundaries distant by more than ${\rm{\Delta }}{\mathrm{log}}_{10}E\,=\,0.06$, which corresponds to an energy resolution of ±7% relevant in the range covered here (Pierre Auger Collaboration 2020b). As illustrated by the set of figures above energy thresholds ranging in 32–80 EeV (see online material attached to Figure 4), the reconstructed parameters do not show significant variations with energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (49 panels)

Tarball of figure set images

For the sake of completeness, we provide the best-fit parameters and maximum TS obtained above energy thresholds corresponding to the global maximum at E ≳ 40 EeV, in the upper part of Table 2, as well as those obtained above the secondary maximum identified at E ≳ 60 EeV, in the lower part of the same table. The most significant departure from isotropy is identified for all four catalogs at energy thresholds in the range 38–40 EeV, with post-trial p-values of 8.3 × 10⁻⁴, 7.9 × 10⁻⁴, 4.2 × 10⁻⁴, and 3.2 × 10⁻⁵ for jetted AGNs traced by their γ-ray emission, galaxies traced by their near-infrared emission, all AGNs traced by their X-ray emission, and starburst galaxies traced by their radio emission, respectively. As in Pierre Auger Collaboration (2018b), we do not penalize for the test of the four catalogs, which all provide similar UHECR flux patterns. As a note, the infrared sample of galaxies contains a large fraction (more than 75%) of each of the three other catalogs and only jetted AGN and starburst catalogs can be considered as strictly distinct galaxy samples.

Table 2. Best-fit Results Obtained with the Four Catalogs at the Global (Upper) and Secondary (Lower) Maximum

Catalog	Eth [EeV]	Fisher Search Radius, Θ [deg]	Signal Fraction, α [%]	TSmax	Post-trial p-value
All galaxies (IR)	40	${16}_{-6}^{+11}$	${16}_{-7}^{+10}$	18.0	7.9 × 10−4
Starbursts (radio)	38	${15}_{-4}^{+8}$	${9}_{-4}^{+6}$	25.0	3.2 × 10−5
All AGNs (X-rays)	39	${16}_{-5}^{+8}$	${7}_{-3}^{+5}$	19.4	4.2 × 10−4
Jetted AGNs (γ-rays)	39	${14}_{-4}^{+6}$	${6}_{-3}^{+4}$	17.9	8.3 × 10−4
All galaxies (IR)	58	${14}_{-5}^{+9}$	${18}_{-10}^{+13}$	9.8	2.9 × 10−2
Starbursts (radio)	58	${18}_{-6}^{+11}$	${19}_{-9}^{+20}$	17.7	9.0 × 10−4
All AGNs (X-rays)	58	${16}_{-6}^{+8}$	${11}_{-6}^{+7}$	14.9	3.2 × 10−3
Jetted AGNs (γ-rays)	58	${17}_{-5}^{+8}$	${12}_{-6}^{+8}$	17.4	1.0 × 10−3

Note. The energy threshold, E_th, Fisher search radius, Θ, and signal fraction, α, which maximize the TS, TS${}_{\max }$, for each of the catalogs. The post-trial p-value accounts for the energy scan and search over α and Θ.

Download table as:  ASCIITypeset image

As discussed in Section 4.1, all four sky models tested here are based on improved versions of the catalogs used by Pierre Auger Collaboration (2018b), although with a mild impact on the significance of the results and no noticeable change in the best-fit parameters. The maximum TS is obtained at the same point of the parameter space using the catalogs of infrared galaxies, starburst galaxies, and X-ray AGNs from Pierre Auger Collaboration (2018b), with TS values of 16.0, 23.1, and 18.0, respectively, differing by less than 2 units from the results in Table 2. The most important change is observed for the gamma-ray catalog of jetted AGNs: the maximum TS (13.5) is obtained above ∼60 EeV with the earlier catalog version based on the 2FHL catalog (E_γ > 50 GeV), while it is obtained above ∼40 EeV with the current version based on the 3FHL catalog (E_γ > 10 GeV). The change can be understood from the lower energy threshold of the 3FHL catalog, which reduces the relative flux of blazars beyond 100 Mpc (Mkn 421 and Mkn 501) with respect to the flux of local radio galaxies (Cen A, NGC 1275, and M 87).

A visual inspection of the sky models displayed in Appendix C highlights the main similarity between the four catalogs, namely a hotspot expected in the Auger field of view in the direction of the group of galaxies composed of the radio galaxy Centaurus A, the Seyfert galaxy NGC 4945, and the starburst galaxy M 83. These three galaxies, at distances of about 4 Mpc, constitute one of the pillars of the so-called Council of Giants (McCall 2014) surrounding the Milky Way and Andromeda galaxy. Inspection of the two AGN models, tracing accretion through X-ray emission and jet activity through γ-ray emission, does not suggest bright secondary hotspots in other sky regions at the highest energies (E ≳ 60 EeV), as attenuation of the UHECR flux dramatically reduces the contribution from more distant galaxies. On the other hand, both the infrared model of stellar mass and the radio model of enhanced star-forming activity suggest hotspots in the directions of other members of the Council of Giants: the starburst galaxies NGC 253 and M 82, which are the only two starburst galaxies currently detected at TeV energies. ¹⁰⁶ While M 82 lies in the blind region of the Pierre Auger Observatory, which can only be observed with the Telescope Array (Telescope Array Collaboration 2018), the contribution from NGC 253 is responsible for the larger departure from isotropy obtained with the starburst model with respect, e.g., to the X-ray AGN model (see Appendix C). The infrared model instead yields a smaller TS than both the X-ray AGN and starburst models. Within the infrared model, the region of the Virgo cluster (at d ∼ 20 Mpc) would be brighter than the Centaurus region, which is in tension with the UHECR observations. Following the same procedure as Pierre Auger Collaboration (2018b), we performed a quantitative comparison between the four models to determine whether one of them is favored by the data against the others. The infrared, X-ray, and γ-ray models fit the data at E ≥ 38–40 EeV poorer than the starburst model with C.L. ≲ 3σ. No firm evidence for a catalog preference is identified.

The deviation from isotropy suggested with all four galaxy catalogs is driven by a hotspot in the Centaurus region. This region shows an enhanced flux in all four sky models, arising mainly from Centaurus A for the two AGN models, NGC 4945 for the starburst model, and from both galaxies in the infrared model. The peak direction of the UHECR hotspot, as identified through the blind search described in Section 3.1, points 29 away from the main contributor to the starburst model, NGC 4945, and 51 away from the main contributor to the AGN models, Centaurus A.

Centaurus A, being the closest radio galaxy at 3.68 ± 0.05 Mpc, has been the target of searches for UHECR excess by the Pierre Auger Collaboration for more than a decade (Pierre Auger Collaboration 2007). We update such searches by performing the same analysis described in Section 3.3 using as target the position of Centaurus A, (α, δ) = (2014, − 430). The map of the local p-values as a function of energy threshold and top-hat search angle is shown in Figure 5. The most significant excess is found at E_th = 38 EeV in a circle of top-hat radius of Ψ = 27°, where the number of observed events is N_obs = 215 while ${N}_{\exp }=152.0$ events would be expected from isotropy. The minimum local p-value, which is estimated as in Section 3 from the binomial probability to observe N_obs or more events from an isotropic distribution, is 2.1 × 10⁻⁷. After penalization for the scan in energy and search angle, the post-trial p-value is 4.5 × 10⁻⁵, similar to that obtained with the likelihood-ratio test for starburst galaxies against isotropy.

Zoom In Zoom Out Reset image size

The best-fit parameters of the search in the direction of Centaurus A are unsurprisingly similar to those of the blind search. The lower post-trial p-value with respect to the blind search results from the direction being fixed a priori, as suggested by the early-day searches from the Pierre Auger Collaboration (Pierre Auger Collaboration 2007, 2010b). The top-hat angular scale inferred from the blind search and from the search at the position of Centaurus A, Ψ = 24°–27°, can be compared to the Fisher search radius inferred from the catalog-based searches through the relation Ψ = 1.59 × Θ. ¹⁰⁷ The catalog-based searches yield Θ = 14° − 16° that corresponds to Ψ = 22° − 25°, i.e., a range of values that is consistent with those inferred from the other searches.

Both the catalog-based searches and search in the Centaurus region point to a most significant signal at an energy threshold close to 40 EeV. This energy range encompasses the flux suppression of the energy spectrum above the toe, at E₃₄ = 46 ± 3 ± 6 EeV (Pierre Auger Collaboration 2020b). The evolution of the signal with energy displayed in Figure 3 appears to be mainly driven by the event distribution in the Centaurus region, as illustrated in Figure 6. The pre-trial p-value in the Centaurus region is obtained by profiling the local p-value against the search radius and penalizing for this free parameter. The profile as a function of energy threshold is compared to the TS of the starburst catalog. The latter is chosen as example, noting that the results obtained with other catalogs show a similar dependence on energy threshold (see Figure 3).

Zoom In Zoom Out Reset image size

Constraints from maximum shower-depths up to a few tens of EeV and from the broadband spectrum above the ankle energy suggest that UHECRs are accelerated in proportion to their charge, following so-called Peters' cycles (Pierre Auger Collaboration 2017c, 2020d). The cosmic-ray composition above the toe in the energy spectrum is then expected to be dominated by UHECRs near a maximum magnetic rigidity, R_cut. Accounting for both systematic uncertainties on the energy and maximum shower-depth scales, we inferred in Pierre Auger Collaboration (2017c) a maximum rigidity of ${\mathrm{log}}_{10}({R}_{\mathrm{cut}}/{\rm{V}})={18.72}_{-0.03}^{+0.04}$ with our reference model. Adopting this value as the typical rigidity of UHECRs above the toe, a lower bound on the charge of the bulk of UHECRs above a given energy threshold can be estimated as ${Z}_{\min }={E}_{\mathrm{th}}/{R}_{\mathrm{cut}}$, as displayed along the top axis of Figure 6. The uncertainties on the points illustrate those on the maximum rigidity in the reference scenario. It should be noted that the UHECR composition at the highest energies remains poorly constrained with Phase 1 data and can only be conjectured from a model-dependent approach at this stage.

At rigidities close to R_cut = 5 EV, i.e., ${\mathrm{log}}_{10}({R}_{\mathrm{cut}}/{\rm{V}})\approx 18.7$, UHECR propagation in the magnetic field of the Milky Way enters into a semi-ballistic regime (Erdmann et al. 2016). Excesses identified in the UHECR sky could thus be used both to track back putative sources and possibly to constrain the configuration and strength of the Galactic magnetic field (see Boulanger et al. 2018, and references therein). The angular scale inferred from the catalog-based search, as well as that from the blind search and search in the Centaurus region, are consistent with the average angular dispersion expected in the Milky Way of the Auger mix of nuclear species (Pierre Auger Collaboration 2018b). Nonetheless, the lack of a significant preference for a specific class of galaxies and the strength of the anisotropy signal, reaching at best post-trial p-values of (3 − 5) × 10⁻⁵, still limit the identification of the host galaxies of UHECR accelerators and UHECR constraints on the Galactic magnetic field.

Although only pieces of evidence for anisotropy on intermediate angular scales can be claimed with the Phase 1 high-energy data set, the continued operation of the array may enable the reach of the 5σ discovery threshold. The latter corresponds to a post-trial p-value of 2.9 × 10⁻⁷ or 5.7 × 10⁻⁷ considering a search for both excesses and deficits (2-sided test) or just for excesses (1-sided test). The growth of the signal in the Centaurus region, quantified by the excess of events with respect to the isotropic expectation, and the growth of the TS of the starburst model are displayed as a function of accumulated exposure in Figure 7. These analyses yield post-trial significances of 3.9–4.2σ for a 1- or 2-sided test applied to the Phase 1 high-energy data set. Both the TS and the excess of events are expected to grow linearly with exposure and the fluctuations observed around such a linear behavior are consistent with those expected from simulations. The model-independent search in the Centaurus region shows the smallest fluctuations and may be the most robust approach to forecasting the evolution of the signal. Assuming a fixed top-hat angular scale of Ψ = 27° and a continued growth of the excess at a rate of 5.2 ± 1.2 events per 10,000 km² yr sr, the 5σ (1-sided) discovery threshold would be expected for a total accumulated exposure of 165,000 ± 15,000 km² yr sr (68% C.L.), which would be within reach by the end of 2025 (±2 calendar years) adopting an approach similar to that developed in the present study.

Zoom In Zoom Out Reset image size

We have presented the measurement and analysis of arrival directions of the highest-energy events detected at the Pierre Auger Observatory during its first phase of operation. With a total of 2635 UHECR events above 32 EeV and an accumulated exposure of 122,000 km² sr yr, no indication for anisotropies on angular scales ranging from one to thirty degrees emerges from autocorrelation studies or from blind searches over the entire sky. This lack of significant deviation from isotropy can be attributed a posteriori to the small amplitude of the anisotropic signal evidenced here, to the vastness of the parameter space that has been probed, in addition to the limited number of events at the highest energies. More focused searches along the Galactic center and Galactic plane do not reveal any excesses. The fluxes along these structures and the associated statistical uncertainty are Φ_GC(≥40 EeV, Ψ = 25°) = (10.9 ± 1.1) × 10⁻³ km⁻² yr⁻¹ sr⁻¹ and Φ_GP(≥40 EeV, Ψ = 25°) = (9.8 ± 0.7) × 10⁻³ km⁻² yr⁻¹ sr⁻¹, respectively. These values can be compared to the average flux over the field of view of the Observatory of Φ_ISO(≥40 EeV) = (11.3 ± 0.4) × 10⁻³ km⁻² yr⁻¹ sr⁻¹. A study along the supergalactic plane, not distinguishing among the various galaxies forming this structure, similarly yields Φ_SGP(≥ 40 EeV, Ψ = 25°) = (9.8 ± 0.6) × 10⁻³ km⁻² yr⁻¹ sr⁻¹.

Accounting for the attenuation of the UHECR mix inferred from lower energy observations, the sky viewed from the Pierre Auger Observatory is better modeled with a ∼10% flux excess in the directions of nearby galaxies observed in the radio, near-infrared, X-ray, and gamma-ray bands. A 1-sided test for an excess disfavors isotropy at the 3.3–4.2σ level, depending on the catalog. A model-independent analysis of the Centaurus region, which contains the most prominent active and star-forming galaxies expected to contribute at these energies, reveals an excess that is significant at the 4.1σ C.L.

The average flux above 40 EeV in a 25° top-hat region centered on Centaurus A can be estimated as Φ_Cen(≥40 EeV, Ψ = 25°) = (15.9 ± 1.3) × 10⁻³ km⁻² yr⁻¹ sr⁻¹. In comparison, regions centered on the Virgo cluster and on the starburst galaxy NGC 253 show fluxes of Φ_Virgo(≥40 EeV, Ψ = 25°) = (12.2 ± 1.8) × 10⁻³ km⁻² yr⁻¹ sr⁻¹ and Φ_{NGC 253}(≥40 EeV, Ψ = 25°) = (12.8 ± 1.2) × 10⁻³ km⁻² yr⁻¹ sr⁻¹. As illustrated by the model sky maps in Appendix C, the regions of NGC 253 and of the Virgo cluster could be expected to be as bright as and brighter than the Centaurus region if the UHECR emission rate was simply traced by the star formation rate and stellar mass, respectively. At the present stage, although the starburst catalog enables the identification of the most significant deviation from isotropy (4.2σ) and the jetted AGN catalog the least significant deviation (3.3σ), no firm preference for correlation with a specific class of galaxies can be stated. It should further be noted that such a preferred correlation would not necessarily suggest causation in the form of the identification of the origin of UHECRs, as regular and turbulent magnetic fields traversed by these charged particles could alter the anisotropic pattern observed on Earth (e.g., Kotera & Lemoine 2008; Erdmann et al. 2016; Farrar & Sutherland 2019; Bell & Matthews 2022).

Though the most significant deviation from isotropy is found at energies around ∼40 EeV for almost all the analyses, the excess is also hinted at for all catalogs and the Centaurus region at energies around ∼60 EeV, as shown in Figure 8 (see online material). Indeed, it was in this higher energy range that the first indication of anisotropy was found in early Auger data (Pierre Auger Collaboration 2007). An interpretation of the energy evolution of the signal on intermediate angular scales could be drawn in terms of the maximum energy achieved for higher-charge nuclei. In a Peters' cycle scenario such as discussed in Section 5, the evidence for anisotropy above ∼40 EeV would be interpreted as stemming from CNO nuclei, which would suggest Z ≈ 10–12 nuclei to be responsible for the departure from isotropy above ∼60 EeV. The estimate of the maximum rigidity used here is based on the combined fit of spectra and maximum depth of shower performed by Pierre Auger Collaboration (2017c). The direct inclusion in such analyses of arrival-direction information will enable us to test more directly this scenario. If this scenario of local extragalactic sources is extrapolated to lower energies, one could expect a contribution from He nuclei (see, e.g., Lemoine & Waxman 2009) in the energy range where a significant dipole, but no significant quadrupole, has been reported using data from the Observatory. The strength of such an anisotropic contribution could nonetheless be further diluted in the contribution from more distant sources. We foresee that an in-depth comparison could be drawn by studying the evolution of the large-scale dipolar and quadrupolar components as a function of energy. ¹⁰⁸ Alternatively, a more model-dependent but also more-constrained approach could exploit full-sky, flux-limited catalogs encompassing galaxies out to the cosmic-ray horizon at the ankle energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (49 panels)

Tarball of figure set images

At this stage, it is not possible to make claims on which are the sources of the highest-energy particles known in the universe. This is in part due to the deflection they suffer in magnetic fields. Identifying the sources of UHECRs indeed runs parallel to deducing the properties of Galactic and extragalactic magnetic fields, and constraints on one of these will enhance our understanding of the other. An important step will be taken through the inclusion of composition-sensitive observables in arrival-direction studies. This will be done either through searches for anisotropy in the moments of such composition observables or by their use, event by event, to select only candidate light nuclei. Future studies using the Observatory offer the promise to do so by means of the AugerPrime upgrade, currently being completed, which will enhance mass discrimination with the 100% duty cycle of the surface detector.

The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment and effort from the technical and administrative staff in Malargüe. We are very grateful to the following agencies and organizations for financial support:

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

We gratefully acknowledge constructive feedback from the anonymous reviewer as well as exchanges with Alberto Dominguez on the 3FHL catalog and exchanges with Cecilia Lunardini, Kimberly Emig, and Rogier Windhorst on their starburst catalog. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. We acknowledge the usage of the HyperLeda database (http://leda.univ-lyon1.fr).

The data set used here consists of 2635 events above 32 EeV collected at the Pierre Auger Observatory from 2004 January 1 to 2020 December 31. The data set is formatted as shown in Table 3, which lists the twenty highest-energy events. For each event, we report the year in which the event was detected, the Julian day of the year, and the time of detection in UTC seconds. The arrival directions are expressed in local coordinates, (θ, ϕ), the zenith and azimuth angle (measured counterclockwise from the east), respectively, and in equatorial coordinates (J2000), (α, δ), the right ascension (R.A.) and declination (decl.), respectively. Finally, the reconstructed energy, in EeV, and the integrated exposure accumulated up to the time of detection are reported in the last two columns. The full list of 2635 events, with the same information as in Table 3, is available at https://doi.org/10.5281/zenodo.6504276 together with the code, the structure of which is described in Appendix B.

The energies and arrival directions of the events may have changed with respect to those already released in previous works, such as Pierre Auger Collaboration (2015a). These changes are due to the refinements in the reconstruction reported in Section 2 and to updates in the energy scale and calibration which were improved over the years. Similarly, a subsample of the vertical events used here is included in the recent data release from the Collaboration (Pierre Auger Collaboration 2021). The latter were derived with the other reconstruction software used in the Collaboration, which enables independent cross checks and shows good consistency with the reconstruction software used here (Pierre Auger Collaboration 2020c).

As mentioned in Section 2, the ratio of the number of inclined and vertical events is energy dependent. Anisotropy itself could impact the ratio of inclined and vertical events, as the two exposures differ over the sky. This effect is however small: the excess reported in Section 5 would imply an expected ratio of N_incl/N_vert = 0.273 instead of 0.278 for an isotropic distribution. Above 32 EeV, a non-significant excess of inclined events is observed with respect to expectations from the exposure ratio and finite energy resolution (N_incl/N_vert = 0.292 ± 0.014). Above 80 EeV, there are 10 events with θ ≥ 60° and 86 with θ < 60°, which corresponds to a ratio of N_incl/N_vert = 0.116 ± 0.039. The deficit of inclined events is most significant above ∼90 EeV, which results in a post-trial significance (under the assumption of isotropy) at the level of ∼2.5σ, when penalized for a search as a function of energy. Such a discrepancy or a stronger one would have a 1.3% probability of being found as a statistical fluctuation under the hypothesis that the energy calibrations of both data streams are correct. For completeness, we also consider the hypothesis that the deficit of inclined events at the highest energies is at least partly due to a systematic underestimation of inclined energies (or overestimation of vertical ones), as different reconstruction techniques are used for the two sets. We tested for this effect empirically by selecting events with zenith angles between 57° < θ < 63° that are reconstructed by both the vertical and inclined reconstructions and for which six active stations surround the one closest to the core position. There are 161 such events and a power-law relation of the form ${E}_{\mathrm{vert}}=A\cdot {E}_{\mathrm{incl}}^{B}$ was fitted to extract the parameters (A, B) that would convert the energies obtained from the inclined reconstruction to the energies obtained from the vertical reconstruction. The results are such that E_vert = 80 EeV would correspond to E_incl = 76.1 ± 1.6 EeV. We applied the change to the energies of events in the inclined data set and performed, as a cross-check, the likelihood analysis with the starburst catalog (as in Section 4.1) and the Centaurus-region analysis (as in Section 5). In both cases, we found the same results presented with the standard data set. This cross-check demonstrates that the possible systematic uncertainties induced by the difference in energy calibration of the vertical and inclined reconstructions do not affect the results presented in this paper.

The structure of the code used to produce the results of this paper is presented in Figure 9. The main analyses are contained in two folders, called Targeted_Blind for Sections 3 and 5, and Catalog_Based for Section 4. The add-ons and utilities needed to run the analyses are contained in the folders Data, Utilities, and Visuals. A brief description of each folder follows.

• 1.  
  The Data folder contains the file of the data set used in all of the analyses, named AugerApJS2022_Yr_JD_UTC_Th_Ph_RA_Dec_E_Expo.dat. Additionally, the folder contains a C++ script named exposure.cc, which computes the directional exposure of the Observatory for both vertical and inclined events, which is integrated over the duration of the acquisition period; the exposure script produces two files, exposure.root and exposure.fits, which contain the exposure as a function of decl. in TF1-root format and in healpixmap-fits form (RING scheme, Galactic coordinates), respectively. The Time_exposure file provides the evolution of exposure with time, as displayed in the upper axis of Figure 7. The DataPath.h contains the declaration of the data set file to be used by all the analyses for easy user intervention.
• 2.  
  The Utilities folder contains files with auxiliary classes and functions used by all other parts of the code, in particular coordinate-conversion utilities and HEALPix map manipulation.
• 3.  
  The Targeted_Blind folder contains the code for the targeted (Sections 5 and 3.3), blind (Section 3.1), and autocorrelation (Section 3.2) analyses. The first folder level contains:

  • (a)  
    Dedicated utilities (utils.h and utils.cc);
  • (b)  
    Three subfolders: Blind, Targeted, and Autocorrelation containing the respective analyses.

  Each of the three subfolders contains a script performing the computation of the local p-value, 1_pvalue.cc, and a script penalizing for the search over the parameter space, 2_penalization.cc. The results of the Blind code are stored in .txt files for easy readout, while the results of the Autocorrelation code are stored as a .root file. In the Targeted folder, 1_pvalue.cc produces the outputs fig5.root, fig2gc.root, fig2gp.root, and fig2sgp.root, which contain the local p-value in bins of energy threshold and search radius for the Centaurus region, Galactic center, Galactic plane, and supergalactic plane analyses, respectively; the .root files in the Autocorrelation and Targeted folders also contain copies of the support histograms used to calculate the local p-value which are used by the penalization algorithm. The script 2_penalization.cc produces the post-trial p-value in single-value text form, as well as the file fig6tmp.txt, which stores the pre-trial p-values as a function of energy threshold penalized only for the scan in angle (see Figure 6).
• 4.  
  The Catalog_Based folder contains the code for the likelihood-ratio analysis. The first folder level comprises:

  • (a)  
    Dedicated utilities (utils.h and utils.cc);
  • (b)  
    The folder Catalogs, which contains the raw catalogs of galaxies in the subfolder Multiwavelength, as described in Appendix C. The raw catalog files are input, above a fixed energy threshold, to the script propagate.cc, in conjunction with the Auger composition model contained in the file AnaCRP3.root, to produce the attenuated models used in the analysis. These attenuated models can be produced above all energy thresholds by running the script propagate.py, with outputs stored in the subfolder ModelsUHECR. The latter is organized as different folders for each catalog;
  • (c)  
    The analysis routines: 1_pre_trial.cc, which produces the results stored in the files fig3.root, showing the TS, signal fraction, and search radius as a function of the threshold energy, and fig4.root, showing the TS as a function of the signal fraction and search radius with 68% C.L. contours for each catalog; 2_penalization.cc produces the post-trial p-values.
• 5.  
  The Visuals folder contains scripts that produce the figures shown in the paper. The python script skymaps.py produces the sky maps in Hammer–Aitoff view: the Li–Ma significance map, fig1.pdf, the flux maps above successive energy threshold stored in fig8.pdf and the model maps stored in fig10.pdf. The script show_figures.py produces fig2.pdf, fig3.pdf, fig4.pdf, and fig5.pdf from their respective root files. The script evolution.cc produces fig6.pdf, a plot of the pre-trial p-values from the Centaurus-region analysis and likelihood-ratio analysis against starburst galaxies as a function of the threshold energy; it also produces fig7.pdf, a plot of the evolution of TS of the starburst analysis and of the excess in the Centaurus region as a function of the exposure accumulated at the Observatory.

Zoom In Zoom Out Reset image size

The best-fit sky models above 40 EeV obtained with the four catalogs described in Section 4.1 are shown in Figure 10. These sky maps do not include any isotropic component and display only the flux expected from galaxies included in the catalogs, which is smeared on the best-fit Fisher angular scale above 40 EeV obtained with each catalog. A further top-hat smoothing on an angular scale of Ψ = 25° is performed for the sake of comparison with Figure 8.

Zoom In Zoom Out Reset image size

The models shown in Figure 10 are based on the UHECR flux expected from each galaxy in proportion to its electromagnetic flux. The multiwavelength information on the galaxies is made available in the Multiwavelength subfolder of the catalog-based study, as described in Appendix B, and is available online at https://doi.org/10.5281/zenodo.6504276. The Multiwavelength folder contains one file per catalog, with tabulated values detailed in Tables 4, 5, 6, and 7. The first column of each of these tables provides the name of the source as referenced by the authors of the source catalog. The second column provides a counterpart name that is consistent across all four catalogs. The third column provides the type of galaxy, extracted either from the source catalog or from the HyperLEDA database. The fourth and fifth columns provide the equatorial coordinates of the galaxy. The sixth and seventh columns display the distance modulus and associated uncertainty extracted from the modbest entry of the HyperLEDA database. The eighth and ninth columns display the corresponding luminosity distance in Mpc as well as the relative uncertainty on this quantity. The electromagnetic flux of each galaxy is provided in column 10, except in Table 4 where the K-band magnitude is provided. Whenever available, the uncertainty on the quantity provided in column 10 is shown in column 11. Finally, a flag is provided in the last columns of Tables 5, 6, and 7. This flag indicates whether the galaxy was also included in the main samples studied by Pierre Auger Collaboration (2018b) (Y), in one of the cross-check samples (X), or not included in earlier versions of these catalogs (N). The flag column of Table 6 indicates the origin of the redshift estimate, either from HyperLEDA or from NED for the 23 X-ray AGNs that are not listed in HyperLEDA.