Probing for anisotropies at different angular scales in the arrival directions of the cosmic rays detected at the Pierre Auger Observatory

L. Deval,⁴⁹ A. di Matteo,⁴⁹ C. Dobrigkeit,²¹ J.C. D’Olivo,⁶⁴ L.M. Domingues Mendes,^16,67 T. Dominguez,¹ Y. Dominguez Ballesteros,²⁷ Q. Dorosti,⁴¹ R.C. dos Anjos,²⁴ J. Ebr,²⁹ F. Ellwanger,³⁷ R. Engel,^35,37 M. Erdmann,³⁸ A. Etchegoyen,^7,12 C. Evoli,^42,43 H. Falcke,^73,75,74 G. Farrar,⁸⁰ A.C. Fauth,²¹ T. Fehler,⁴¹ F. Feldbusch,³⁶ A. Fernandes,⁶⁷ M. Fernández Alonso,¹⁴ B. Fick,⁷⁹ J.M. Figueira,⁷ P. Filip,^35,7 A. Filipčič,^71,70 B. Flaggs,⁸² A. Franco,⁴⁵ M. Freitas,⁶⁷ T. Fujii,^81,j A. Fuster,^7,12 C. Galea,⁷³ C. Galelli,⁴⁶ B. García,⁶ C. Gaudu,³⁴ P.L. Ghia,³¹ U. Giaccari,⁴⁵ M. Giammarco,^53,43 C. Glaser,³⁹ F. Gobbi,¹⁰ F. Gollan,⁷² G. Golup,¹ P.F. Gómez Vitale,¹¹ J.P. Gongora,¹¹ N. González,⁷ D. Góra,⁶⁵ A. Gorgi,^50,49 M. Gottowik,³⁷ F. Guarino,^56,47 G.P. Guedes,²² Y.C. Guerra,¹⁰ S. Hahn,³⁵ P. Hamal,²⁹ M.R. Hampel,⁷ P. Hansen,³ V.M. Harvey,¹³ A. Haungs,³⁷ M. Havelka,²⁹ T. Hebbeker,³⁸ C. Hojvat,^d J.R. Hörandel,^73,74 P. Horvath,³⁰ M. Hrabovský,³⁰ T. Huege,^37,15 A. Insolia,^54,44 P.G. Isar,⁶⁹ M. Ismaiel,^73,74 P. Janecek,²⁹ V. Jilek,²⁹ K.-H. Kampert,³⁴ L. Karczubová,²⁹ B. Keilhauer,³⁷ V.V. Kizakke Covilakam,⁷ H.O. Klages,³⁷ M. Kleifges,³⁶ A. Klingel,²⁹ J. Köhler,³⁷ F. Krieger,³⁸ M. Kubatova,²⁹ B.L. Lago,¹⁷ N. Leal,⁷ M.A. Leigui de Oliveira,²³ Y. Lema-Capeans,⁷² A. Letessier-Selvon,³² I. Lhenry-Yvon,³¹ L. Lopes,⁶⁷ P. Magnani,^42,43 M. Mallamaci,^57,44 S. Mancuso,^50,49 D. Mandat,²⁹ P. Mantsch,^d A.G. Mariazzi,³ C. Marinelli,^42,43 I.C. Mariş,¹⁴ G. Marsella,^57,44 D. Martello,^52,45 S. Martinelli,³⁷ O. Martínez Bravo,⁶⁰ A. Martínez-Mendez,²⁷ M.A. Martins,²⁹ H.-J. Mathes,³⁷ J. Matthews,^g G. Matthiae,^58,48 E. Mayotte,⁷⁷ S. Mayotte,⁷⁷ P.O. Mazur,^d G. Medina-Tanco,⁶⁴ D. Melo,⁷ A. Menshikov,³⁶ C. Merx,³⁷ S. Michal,²⁹ M.I. Micheletti,⁵ L. Miramonti,^55,46 M. Mogarkar,⁶⁵ S. Mollerach,¹ F. Montanet,^h L. Morejon,³⁴ K. Mulrey,^73,74 R. Mussa,⁴⁹ W.M. Namasaka,³⁴ S. Negi,²⁹ L. Nellen,⁶⁴ G.A. Neves,⁶⁷ K. Nguyen,⁷⁹ G. Nicora,⁹ M. Niechciol,⁴¹ D. Nitz,⁷⁹ D. Nosek,²⁸ A. Novikov,⁸² V. Novotny,²⁸ L. Nožka,³⁰ A. Nucita,^52,45 L.A. Núñez,²⁷ S.E. Nuza,⁴ J. Ochoa,^7,37 M. Olegario,¹⁹ C. Oliveira,²⁰ L. Östman,²⁹ M. Palatka,²⁹ J. Pallotta,⁹ G. Parente,⁷² T. Paulsen,³⁴ M. Pech,²⁹ J. Pękala,⁶⁵ R. Pelayo,⁶¹ D. Perez,¹ C. Pérez Bertolli,⁷² L. Perrone,^52,45 S. Petrera,^42,43 T. Pierog,³⁷ M. Pimenta,⁶⁷ M. Platino,⁷ P. Privitera,⁸¹ C. Priyadarshi,⁶⁵ ⁷ Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), Buenos Aires, Argentina

⁸ International Center of Advanced Studies and Instituto de Ciencias Físicas, ECyT-UNSAM and CONICET, Campus Miguelete – San Martín, Buenos Aires, Argentina

⁹ Laboratorio Atmósfera – Departamento de Investigaciones en Láseres y sus Aplicaciones – UNIDEF (CITEDEF-CONICET), Argentina

¹⁰ Observatorio Pierre Auger, Malargüe, Argentina

¹¹ Observatorio Pierre Auger and Comisión Nacional de Energía Atómica, Malargüe, Argentina

¹² Universidad Tecnológica Nacional – Facultad Regional Buenos Aires, Buenos Aires, Argentina

¹³ Adelaide University, Adelaide, S.A., Australia

¹⁴ Université Libre de Bruxelles (ULB), Brussels, Belgium

¹⁵ Vrije Universiteit Brussels, Brussels, Belgium

¹⁶ Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, RJ, Brazil

¹⁷ Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Petropolis, Brazil

¹⁸ Universidade de São Paulo, Escola de Engenharia de Lorena, Lorena, SP, Brazil

¹⁹ Universidade de São Paulo, Instituto de Física de São Carlos, São Carlos, SP, Brazil

²⁰ Universidade de São Paulo, Instituto de Física, São Paulo, SP, Brazil

²¹ Universidade Estadual de Campinas (UNICAMP), IFGW, Campinas, SP, Brazil

²² Universidade Estadual de Feira de Santana, Feira de Santana, Brazil

²³ Universidade Federal do ABC, Santo André, SP, Brazil

²⁴ Universidade Federal do Paraná, Setor Palotina, Palotina, Brazil

²⁵ Universidade Federal do Rio de Janeiro, Instituto de Física, Rio de Janeiro, RJ, Brazil

²⁶ Universidad de Medellín, Medellín, Colombia

²⁷ Universidad Industrial de Santander, Bucaramanga, Colombia

²⁸ Charles University, Faculty of Mathematics and Physics, Institute of Particle and Nuclear Physics, Prague, Czech Republic

²⁹ Institute of Physics of the Czech Academy of Sciences, Prague, Czech Republic

³⁰ Palacky University, Olomouc, Czech Republic

³¹ CNRS/IN2P3, IJCLab, Université Paris-Saclay, Orsay, France

³² Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Sorbonne Université, Université de Paris, CNRS-IN2P3, Paris, France

³³ Université Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France

³⁴ Bergische Universität Wuppertal, Department of Physics, Wuppertal, Germany

³⁵ Karlsruhe Institute of Technology (KIT), Institute for Experimental Particle Physics, Karlsruhe, Germany

³⁶ Karlsruhe Institute of Technology (KIT), Institut für Prozessdatenverarbeitung und Elektronik, Karlsruhe, Germany

³⁷ Karlsruhe Institute of Technology (KIT), Institute for Astroparticle Physics, Karlsruhe, Germany

³⁸ RWTH Aachen University, III. Physikalisches Institut A, Aachen, Germany

³⁹ TU Dortmund University, Department of Physics, Dortmund, Germany

⁴⁰ Universität Hamburg, II. Institut für Theoretische Physik, Hamburg, Germany

⁴¹ Universität Siegen, Department Physik – Experimentelle Teilchenphysik, Siegen, Germany ⁴²Gran Sasso Science Institute, L’Aquila, Italy

⁴³INFN Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy

⁴⁴INFN, Sezione di Catania, Catania, Italy

⁴⁵INFN, Sezione di Lecce, Lecce, Italy

⁴⁶INFN, Sezione di Milano, Milano, Italy

⁴⁷INFN, Sezione di Napoli, Napoli, Italy

⁴⁸INFN, Sezione di Roma “Tor Vergata”, Roma, Italy

⁴⁹INFN, Sezione di Torino, Torino, Italy

⁵⁰Osservatorio Astrofisico di Torino (INAF), Torino, Italy

⁵¹Politecnico di Milano, Dipartimento di Scienze e Tecnologie Aerospaziali, Milano, Italy

⁵²Università del Salento, Dipartimento di Matematica e Fisica “E. De Giorgi”, Lecce, Italy

⁵³Università dell’Aquila, Dipartimento di Scienze Fisiche e Chimiche, L’Aquila, Italy

⁵⁴Università di Catania, Dipartimento di Fisica e Astronomia “Ettore Majorana“, Catania, Italy

⁵⁵Università di Milano, Dipartimento di Fisica, Milano, Italy

⁵⁶Università di Napoli “Federico II”, Dipartimento di Fisica “Ettore Pancini”, Napoli, Italy

⁵⁷Università di Palermo, Dipartimento di Fisica e Chimica ”E. Segrè”, Palermo, Italy

⁵⁸Università di Roma “Tor Vergata”, Dipartimento di Fisica, Roma, Italy

⁵⁹Università Torino, Dipartimento di Fisica, Torino, Italy

⁶⁰Benemérita Universidad Autónoma de Puebla, Puebla, México

⁶¹Unidad Profesional Interdisciplinaria en Ingeniería y Tecnologías Avanzadas del Instituto Politécnico Nacional (UPIITA-IPN), México, D.F., México

⁶²Universidad Autónoma de Chiapas, Tuxtla Gutiérrez, Chiapas, México

⁶³Universidad Michoacana de San Nicolás de Hidalgo, Morelia, Michoacán, México

⁶⁴Universidad Nacional Autónoma de México, México, D.F., México

⁶⁵Institute of Nuclear Physics PAN, Krakow, Poland

⁶⁶University of Łódź, Faculty of High-Energy Astrophysics,Łódź, Poland

⁶⁷Laboratório de Instrumentação e Física Experimental de Partículas – LIP and Instituto Superior Técnico – IST, Universidade de Lisboa – UL, Lisboa, Portugal

⁶⁸“Horia Hulubei” National Institute for Physics and Nuclear Engineering, Bucharest-Magurele, Romania

⁶⁹Institute of Space Science, Bucharest-Magurele, Romania

⁷⁰Center for Astrophysics and Cosmology (CAC), University of Nova Gorica, Nova Gorica, Slovenia

⁷¹Experimental Particle Physics Department, J. Stefan Institute, Ljubljana, Slovenia

⁷²Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de Santiago de Compostela, Santiago de Compostela, Spain

⁷³IMAPP, Radboud University Nijmegen, Nijmegen, The Netherlands

⁷⁴Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF), Science Park, Amsterdam, The Netherlands

⁷⁵Stichting Astronomisch Onderzoek in Nederland (ASTRON), Dwingeloo, The Netherlands

⁷⁶Case Western Reserve University, Cleveland, OH, USA

⁷⁷Colorado School of Mines, Golden, CO, USA

Introduction

The study of the arrival-direction distribution of ultra-high-energy cosmic rays (UHECRs) is central to understanding their origin. These analyses are classified into two regimes based on the angular scale: large-scale anisotropies and small- to intermediate-scale anisotropies.

At large angular scales, a significant dipolar modulation in right ascension (R.A.) was observed above 8 EeV¹ by the Pierre Auger Collaboration. The modulation was first reported in [1], and has now reached $6.8\sigma$ significance [2]. Its direction lies around 115° away from the Galactic Center, suggesting that UHECRs at these energies are mostly of extragalactic origin. Moreover, the dipole amplitude has an approximately linear growth with energy [2].

At small to intermediate angular scales, no significant anisotropy has been observed in the data collected with the two largest detector arrays for UHECRs: the Pierre Auger Observatory [3] and the Telescope Array [4]. Nonetheless, both collaborations reported the presence of several hints of anisotropies above a few tens of EeV, see for example [5, 6, 7, 8]. The most significant one is an excess observed above 38 EeV within a top-hat circle of 27° centered on the radiogalaxy Centaurus A with a post-trial significance of $4.0\sigma$ [5, 6]. Interestingly, the Pierre Auger Collaboration also observed overdensities in the Centaurus region above 20 EeV, by scanning circular regions crossing the Supergalactic plane with a fixed angular radius of 27° [9]. The scan was repeated for different energy thresholds, and notably, an excess present in the Centaurus region is observed above each energy threshold. Moreover, the central position of the excess region is stable within few degrees as a function of the energy threshold.

The above results were obtained by treating the two angular-scale regimes separately. In this work, we introduce a novel approach that enables a search for overdensities across all angular scales. This analysis is complementary to the power-spectrum analysis which is sensitive to the combination of excesses at a given angular scale all over the sky [2, 10]. With the approach presented here, we look instead for the most significant excess of events and provide a more comprehensive characterization thereof. As a result, we get new insight into the transition from large-scale anisotropies at lower energies to intermediate-scale excesses at higher energies.

We structure this paper as follows. In § 2, we describe the dataset used for the analysis. In § 3, we give a description of the likelihood and the test statistics used in the analysis. In § 4, we report the results obtained with the analysis focusing first, in § 4.1, on a set of energy bins already introduced in previous analyses of the Pierre Auger Collaboration, then, in § 4.2, we perform the analysis introducing additional energy bins. Finally, a discussion of the results and the conclusion are presented in § 5.

The dataset

The Pierre Auger Observatory is the largest facility in the world to detect UHECRs. It is located in Argentina, near the city of Malargüe ($35.2^\circ$ S, $69.5^\circ$ W) at an average altitude of 1,400 m a.s.l. and it covers an area of 3,000 km$^2$. The Observatory is a multi-hybrid system, composed of a surface detector (SD) and a fluorescence detector (FD). The SD is an array of water-Cherenkov detectors located in a triangular grid. The baseline array for the study of the highest-energy cosmic rays has 1600 detectors with a spacing of 1,500 m (SD-1500). The SD has a duty cycle of $\sim100\%$. The FD consists of 24 fluorescence telescopes distributed between 4 sites along the perimeter of the SD-1500 and overlooking the atmosphere above the array. The FD operates only in clear moonless nights, with a duty cycle of $\sim15\%$.

In this work, we adopt the same data selection as it was defined in [2]. We consider only data registered with the SD-1500 between 1st January 2004 and 31st December 2022. For the years 2021 and 2022, when the installation of the AugerPrime upgrade [11] was underway, we consider only data registered with stations with non-upgraded electronics.

The SD-1500 is fully efficient above 4 EeV for events entering the atmosphere with zenith angle $\theta< 80^\circ$ (achieving 85% of sky coverage). According to $\theta$, physical events are reconstructed using different algorithms. The vertical reconstruction is applied to events entering the atmosphere with $\theta< 60^\circ$, while the inclined reconstruction is applied to events with $60^\circ\le\theta< 80^\circ$. Vertical events are included in the analysis if the SD station with the largest signal is surrounded by at least five functioning stations. Inclined events are included if the SD station closest to the reconstructed core is surrounded by at least five functioning stations. Additionally, we require that the reconstructed core falls in an isosceles triangle of nearby functioning stations. Events detected during periods of unreliable data acquisition are removed from the analysis. The total exposure over the whole period is 123,000 km$^2$ sr yr.

The arrival directions of vertical events are derived by fitting a spherical model to the arrival times of the particles contained in the front of the showers at the SD detectors. Given that the shower front of the inclined events is mostly formed by muons, the arrival direction is obtained by fitting the arrival times with a model considering the muon propagation from the production site. A detailed description of the reconstruction procedure for the vertical events can be found in [12], while a description for the inclined events is present in [13]. In both cases, the angular resolution is better than $1^\circ$.

The energy estimation is performed differently for vertical and inclined events. For vertical showers, the proxy for the primary energy is the signal of water-Cherenkov detectors estimated at a reference distance of 1000 m, $S(1000)$. For inclined showers, the proxy is the so-called $N_{19}$: the scaling factor by which the 2D muon density at ground simulated for proton-induced showers at $E = 10^{19}$ eV needs to be multiplied to fit the signal of the water-Cherenkov detectors. Since the particle density at ground depends on the zenith angle due to the attenuation of shower particles in the atmosphere, both $S(1000)$ and $N_{19}$ need to be converted into zenith-independent energy estimators using the constant intensity cut method [14], by which $S(1000)$ ($N_{19}$) is converted to the value expected for a shower arriving with a zenith angle of $38^\circ$ ($68^\circ$). The energy estimator $S_{38}$ ($N_{68}$) is calibrated using a subset of high-quality hybrid events (i.e. events observed simultaneously by SD and FD). Using the hybrid data, we perform a quasi-calorimetric measurement of the shower energy. Therefore, the calibration procedure enables a good energy estimation for all the events detected with the SD. The statistical uncertainty on the energy calibration is $\sim15\%$ for events with energy above 4 EeV and reduces to $\sim7\%$ at the highest energies, above 20 EeV. The systematic uncertainty is $\sim14\%$ [13, 15].

Atmospheric [16] and geomagnetic corrections [17] are applied to the energies of vertical events, to avoid spurious effects in the distribution in R.A. and azimuth angle. Inclined events do not require atmospheric corrections, since atmospheric effects are negligible for muons, while geomagnetic effects are included in the reconstruction procedure.

Description of the analysis

In this work, we characterize the anisotropies in the arrival-direction distribution of the events using a likelihood-ratio test. To construct it, we start by considering a localized source signal distribution superposed over an isotropic background, described by

The first term, $(1-\alpha)$, represents the isotropic fraction, where $\alpha$ is the signal fraction. The second term follows a von Mises–Fisher distribution, which is the generalization of a Gaussian distribution on a sphere. The angular distance measured from the central position of the excess is denoted by $\Theta$. The concentration parameter $\kappa$ controls the steepness of the profile around the reference direction.

The distribution in equation (1) exhibits two distinct asymptotical regimes depending on the value of $\kappa$. In the limit $\kappa\ll1$, the distribution is nearly isotropic and can be expanded in dipole and quadrupole moments

with dipole and quadrupole amplitudes $d=\alpha\kappa$ and $q=\alpha\kappa^2/2$, respectively. We note that $d$ and $q$ represent the dipole and quadrupole amplitudes associated to a modulation centered in the excess position. Since the development is dominated by the dipolar component, $d$ can be compared with previous findings of the Pierre Auger Collaboration [1, 2]. However, the quadrupole term $q$ describes the axisymmetric component ($a_{20}$) with axis in the excess direction, and this is not directly comparable with the amplitudes derived in previous analyses that tested the presence of a quadrupolar modulation [2, 18].

As $\kappa$ increases, the steepness of the von Mises–Fisher term of the distribution becomes more pronounced. For $\kappa\gg1$ it can be approximated, in a small angular region around the center of the excess, by a Gaussian profile with width $\sigma=1/\sqrt{\kappa}$.

In this regime, the distribution can be used to test whether a fraction $\alpha$ of events is clustered within an angular scale $1/\sqrt{\kappa}$ around the reference direction.

Equation (1) provides an accurate description of the expected arrival-direction distribution of the events emitted by a single source located in the reference direction, and propagating through a turbulent magnetic field [19]. Large values of $\kappa$ correspond to quasi–ballistic propagation, whereas small values describe magnetic turbulent diffusion. With this model, we can simultaneously probe dipolar patterns and the transition between large-angular scale anisotropies and small- to intermediate angular scale anisotropies.

For a dataset of $N$ events detected with the Pierre Auger Observatory, the likelihood of the events being distributed according to equation (1) is given by

where $\mathbf{u}_0$ denotes the reference direction, $\mathbf{u}_j$ is the arrival direction of the $j$-th event, and $\omega(\mathbf{u}_j)$ accounts for the directional exposure of the array. The geometrical exposures for vertical events and inclined events are evaluated separately. To keep the analysis as data-driven as possible, $\omega(\mathbf{u}_j)$ is evaluated combining the vertical and inclined geometrical exposures using the ratio of inclined to vertical events observed in each energy bins as the expected exposure ratio for that bin. $\mathcal{N}(\alpha,\kappa,\mathbf{u}_0)$ is a normalization factor depending on the signal fraction, the concentration parameter of the Von Mises–Fisher term and the reference direction

Under the isotropic hypothesis ($\alpha=0$), the likelihood reduces to

The likelihood-ratio between $\mathcal{L}(\alpha,\kappa,\mathbf{u}_0)$ and its value obtained for the isotropic hypothesis $\mathcal{L}(\alpha=0)$ defines the test statistics ($TS$) that we employ in the analysis

As outlined above, the concentration parameter $\kappa$ controls the angular width of the tested distribution. In this work, we compute a profile likelihood as a function of $(\kappa)$, i.e., for each value of $\kappa=10^{-1.0},10^{-0.9},10^{-0.8},\ldots,10^{+3.5}$, we find the $\alpha$ and $\mathbf{u}_0$ that minimize $-2\ln(\mathcal{L}(\alpha,\mathbf{u}_0;\kappa))$. Since the equivalent dipole amplitude is $d=\alpha\kappa$ and the quadrupole amplitude $q=\alpha\kappa^2/2$, with the selected lower limit $\kappa=10^{-1}$, we are able to probe a purely dipolar flux (with a quadrupolar term a factor 20 times smaller than the dipolar one). On the other hand, the upper limit of the $\kappa$ range is selected to match the angular resolution of the Observatory, since $\kappa=10^{3.5}$ corresponds to the equivalent Gaussian width $\sigma=1.0^\circ$.

The minimization is performed using the differential evolution algorithm present in the python library scipy [20]. In the minimization, the declination of $\mathbf{u}_0$ is constrained to be below $+44.8^\circ$, requiring the center of the von Mises–Fisher distribution to lie within the field of view of the Observatory. Although a von Mises–Fisher distribution centered outside the field of view can still produce a non-negligible contribution to the observed event distribution, restricting the center to the field of view avoids border effects, whereby a small number of clustered events close to the edge of the field of view could result in an artificially large signal. Since we profile over $\kappa$, the likelihood described in equation (3) has three degrees of freedom, given by the signal fraction and the celestial coordinates of the reference position ($\mathcal{L}(\alpha,u_0(\mathrm{R.A.},\mathrm{dec});\kappa)$). At any given $\kappa$, for all parameters, the uncertainties are derived from the confidence intervals obtained from the likelihood fit, in which all parameters are varied simultaneously during the minimization.

For every step in the $\kappa$-parameter space, we evaluate the isotropic chance probability $p_{\mathrm{iso}}$ of obtaining the maximum TS. $p_{\mathrm{iso}}$ is evaluated performing 10,000 isotropic Monte-Carlo simulations, and it is defined as the fraction of isotropic simulations yielding a maximum TS at a fixed $\kappa$ equal or larger than the TS obtained on data for the same $\kappa$. Whenever we do not obtain a maximum TS in simulations larger than the TS obtained on data, we evaluate $p_{\mathrm{iso}}$ by doing an extrapolation using a $\chi^2$ distribution with three degrees of freedom. This procedure allows us to probe the most significant excess in the observed arrival-direction distribution across different angular scales, providing a more detailed characterization of how these excesses evolve.

For each energy interval that we consider, we report either $d$ (if $\kappa< 1$) or $\sigma$ (if $\kappa\ge1$) of the excess with the lowest $p_{\mathrm{iso}}$. The values of $d$ are calculated from using the combination of $\alpha$ and $\kappa$ that maximizes the TS, with the center of the excess fixed to the position corresponding to the minimum $p_{\mathrm{iso}}$. This procedure allows us to probe the dipolar modulation centered in the excess position that best describes the data. The uncertainties on $d$ are derived from the profile likelihood obtained by profiling the TS over $\alpha$

where the minimization is performed within the allowed parameter ranges $\alpha\in[0,1]$ and $\kappa\in[0.1,1]$. The confidence interval on $d$ is then determined from the values for which $\Delta\mathrm{TS}_{\mathrm{prof}}(d)=1$.

In this work, we do not account for the small time-dependent non-uniformities in the array exposure, arising from the gradual deployment of the detector and the occasional downtime of individual stations, nor for the effect of the tilt of the array, both of which were included in some previous analyses of large-scale anisotropies by the Pierre Auger Collaboration (see e.g. [2]). Both effects have a minimal impact on the analysis performed here ($\le0.1\%$ in the evaluation of the dipole amplitude), therefore, including them does not change the results.

Results

Results for standard energy bins

We present the results of the analysis described in § 3. In this section, we show the results obtained by dividing the dataset into four energy intervals: (4–8, 8–16, 16–32, $>32$) EeV. These intervals follow the binning adopted in previous studies performed by the Pierre Auger Collaboration for large angular scales [1, 2]. The results corresponding to the lowest $p_{\mathrm{iso}}$, obtained by profiling over $\kappa$ in each energy bin, are listed in Table 1. The values of $p_{\mathrm{iso}}$ marked with an asterisk (*) are derived with the $\chi^2$ extrapolation.

$E/\mathrm{EeV}$	$\kappa$	$\alpha$ (%)	R.A./$^\circ$	dec/$^\circ$	TS	$d$ (%)	$\sigma/^\circ$	$p_{\mathrm{iso}}$
4–8	2.5	$0.7^{+0.2}_{-0.2}$	$112^{+26}_{-33}$	$-65^{+16}_{-13}$	10.1	—	36.2	0.023
8–16	0.10	$67^{+10}_{-10}$	$94^{+9}_{-9}$	$-28^{+10}_{-10}$	49.2	$6.7^{+1.0}_{-1.0}$	—	$1.1 \times10^{-10}$*
16–32	0.10	$100_{-17}$	$95^{+17}_{-16}$	$-53^{+11}_{-10}$	25.2	$10^{+2}_{-2}$	—	$1.3 \times10^{-5}$*
32–$\infty$	0.16	$90^{+10}_{-25}$	$141^{+19}_{-19}$	$-44^{+18}_{-15}$	13.8	$14^{+4}_{-4}$	—	$3.0 \times10^{-3}$

In Table 1, we quote either the equivalent dipole amplitude $d=\alpha\kappa$ (if $\kappa< 1$) or the equivalent Gaussian width $\sigma=1/\sqrt{\kappa}$ (if $\kappa\ge1$). Dipole amplitudes are reported for the (8–16, 16–32, $>32$) EeV intervals. The resulting dipole amplitudes are compatible with those obtained from the 3D dipole reconstruction reported in [2], and the lowest $p_{\mathrm{iso}}$ is found in the (8–16) EeV interval.

In Figure 1, the full results obtained by profiling over $\kappa$ are shown. In the plots, the maximum TS derived for each scanned $\kappa$ is shown together with the 68%, 95%, and 99.7% confidence level (C.L.) bands derived from 10,000 isotropic simulations. We note that the C.L. bands get broader with increasing $\kappa$. Indeed, at larger $\kappa$ (smaller angular widths of the distribution), higher TS can be obtained from isotropic simulations than at smaller $\kappa$, given that more independent regions are probed as a function of increasing $\kappa$, making it more likely to find larger values of TS by chance. Moreover, few clustered events can lead to large TS for small angular scales. With the introduction of $p_{\mathrm{iso}}$, we are able to compute the significance of the excesses observed at different values of $\kappa$. In the (4–8) EeV interval, no TS exceeding the 99.7% C.L. band is observed. In contrast, for (8–16, 16–32, $>32$) EeV, values of TS that exceed the 99.7% C.L. interval are observed. These values correspond to the results reported in Table 1. A secondary excess of the TS outside the 95% C.L. band is observed in ($>32$) EeV at $\kappa= 10^{1.1} \simeq13$. The equivalent Gaussian $\sigma$ is $16^\circ$. The reference direction associated to this excess is (R.A., dec) = ($209^\circ$, $-50^\circ$). The equivalent Gaussian width derived with the minimization can be compared to the top-hat radius $\Psi$ of the Centaurus A targeted analysis through the relation $\Psi= 1.59 \times\sigma$. The Gaussian width obtained is compatible with a top-hat excess with radius $\Psi= 25^\circ$. Therefore, the excess is compatible with the previous finding of the Pierre Auger Collaboration of an overdensity of intermediate angular scale within the $27^\circ$ circular region around Centaurus A, (R.A., dec) = ($201.4^\circ$, $-43.0^\circ$) [5, 6]. No excess with a TS above 95% C.L. is found at concentrations larger than 100 in any of the four energy intervals considered.

In the left column of Figure 2, we show maps in Galactic coordinates of the lowest $p_{\mathrm{iso}}$ obtained by maximizing equation (6) while profiling over $\kappa$ and the reference direction $\mathbf{u}_0$. In these maps, we compute the profile likelihood as a function of $\kappa$ and $\mathbf{u}_0$. For each reference direction, using the center positions of a HEALPix grid [21] defined with the parameter nSide = 64 and for each $\kappa$, we find the $\alpha$ that minimizes $-2\ln(\mathcal{L}(\alpha;\mathbf{u}_0,\kappa))$. We then evaluate $p_{\mathrm{iso}}$ associated to each scanned value of $\kappa$. In each pixel of the maps, we show the lowest $p_{\mathrm{iso}}$ derived.

In the right column of the same figure, the corresponding maps with the best-fit $\kappa$ are depicted. In this second set of maps, the opacity of each point is scaled according to the ratio between the local $\log_{10} p_{\mathrm{iso}}$ value and the minimum $\log_{10} p_{\mathrm{iso}}$ in the map. This procedure suppresses large-$p_{\mathrm{iso}}$ areas, with white areas corresponding to values of $p_{\mathrm{iso}} \simeq1$. Dashed curves represent the boundaries of the field of view of the Pierre Auger Observatory. For the (8–16) EeV and (16–32) EeV intervals, the minimum $p_{\mathrm{iso}}$ is associated with low values of $\kappa$. For ($>$32) EeV, two distinct regions appear: one at low $\kappa$, around the dipole direction, and a second at intermediate $\kappa$, located within the $27^\circ$ circular region around Centaurus A (we note that the position of Centaurus A in Galactic coordinates is $(l,b)=(309.5^\circ, 19.4^\circ)$).

Results for additional energy bins

In this section, we analyze the data using a sliding window spanning a factor of 2 in energy. The $i$-th energy interval ($E_{\min,i}$–$E_{\max,i}$) EeV is defined as

both rounded to the nearest EeV; and $E_{\max,10} = +\infty$. A total of 11 energy bins are selected with the sliding window procedure, the last one being a cumulative bin. The full set of intervals is (4–8, 5–10, 6–13, 8–16, 10–20, 13–25, 16–32, 20–40, 25–51, 32–64, $>40$) EeV. We underline that this sliding window in energy includes the three non-open energy intervals adopted in previous anisotropy studies on large angular scales, i.e. (4–8, 8–16, 16–32) EeV.

The results derived by performing the scan in $\kappa$ yielding the lowest $p_{\mathrm{iso}}$ for each of the energy bins defined with the sliding window procedure are listed in Table 2. As in Table 1, we report the value of $\kappa$, the best-fit estimators, and the corresponding TS for the outcome yielding the lowest $p_{\mathrm{iso}}$ in each energy interval. The equivalent dipole amplitudes or Gaussian widths are quoted depending on whether $\kappa< 1$ or $\kappa\ge1$.

$E$ / EeV	$\kappa$	$\alpha$ (%)	R.A. /$^\circ$	dec /$^\circ$	TS	$d$ (%)	$\sigma$ /$^\circ$	$p_{\mathrm{iso}}$
4–8	2.5	$0.7^{+0.2}_{-0.2}$	$112^{+26}_{-33}$	$-65^{+16}_{-13}$	10.1	—	36.2	0.023
5–10	2.0	$1.4^{+0.3}_{-0.3}$	$129^{+12}_{-13}$	$-40^{+13}_{-12}$	18.5	—	40.6	$4.6 \times10^{-4}$
6–13	0.10	$39^{+8}_{-8}$	$98^{+12}_{-12}$	$-28^{+14}_{-13}$	27.4	$3.8^{+0.8}_{-0.8}$	—	$4.9 \times10^{-6*}$
8–16	0.10	$67^{+10}_{-10}$	$94^{+9}_{-9}$	$-28^{+10}_{-10}$	49.2	$6.7^{+1.0}_{-1.0}$	—	$1.1 \times10^{-10*}$
10–20	0.10	$69^{+12}_{-12}$	$95^{+11}_{-11}$	$-30^{+12}_{-11}$	34.1	$6.9^{+1.2}_{-1.2}$	—	$1.9 \times10^{-7*}$
13–25	0.10	$87^{+13}_{-18}$	$87^{+15}_{-15}$	$-50^{+12}_{-10}$	26.9	$8.7^{+1.8}_{-1.7}$	—	$5.9 \times10^{-6*}$
16–32	0.10	$100_{-17}$	$95^{+17}_{-16}$	$-53^{+11}_{-10}$	25.2	$10^{+2}_{-2}$	—	$1.3 \times10^{-5*}$
20–40	0.13	$95^{+5}_{-22}$	$100^{+18}_{-18}$	$-54^{+12}_{-12}$	21.2	$12^{+3}_{-3}$	—	$8.9 \times10^{-5*}$
25–51	0.13	$100_{-19}$	$108^{+17}_{-17}$	$-45^{+14}_{-13}$	19.4	$13^{+4}_{-3}$	—	$2.3 \times10^{-4}$
32–64	0.16	$95^{+5}_{-26}$	$143^{+18}_{-18}$	$-39^{+18}_{-15}$	14.4	$15^{+5}_{-4}$	—	$2.3 \times10^{-3}$
40–$\infty$	16	$2.7^{+0.7}_{-0.6}$	$210^{+7}_{-7}$	$-48^{+4}_{-4}$	22.9	—	14.4	$2.4 \times10^{-4}$

Dipole amplitudes are derived for the energy intervals (6–13, 8–16, 10–20, 13–25, 16–32, 20–40, 25–51, 32–64) EeV. We note that the dipole amplitude increases with energy, in agreement with the trend reported in [2]. The corresponding $p_{\mathrm{iso}}$ decreases rapidly with energy, reaching the lowest value for (8–16) EeV, and then it gradually increases at higher energies.

For ($>40$) EeV, the value of $\kappa\simeq16$ corresponding at the lowest $p_{\mathrm{iso}}$ describes an excess of intermediate angular scale ($\sigma= 14.4^\circ$). The central position of the excess is located within the $27^\circ$ region around Centaurus A.

The full sets of values of $p_{\mathrm{iso}}$ obtained after profiling over $\kappa$ are shown in Figure 3. The energy intervals are grouped into four plots. We note that the profiling over $\kappa$ in the (5–

10. EeV interval yields to results with $p_{\mathrm{iso}}$ outside the 99.7% C.L. starting at $\kappa= 10^{-1}$, and the lowest $p_{\mathrm{iso}}$ is observed for $\kappa= 2$. With the exception of the (4–8), (>32), and (>40) EeV interval, all energy bins exhibit excesses outside of the 99.7% C.L. band at $\kappa< 1$.

In the interval (32–64) EeV, the most significant excess observed is at $\kappa\simeq0.16$, outside the 99.7% C.L. band. Moreover, a second excess outside the 95% C.L. band is present at $\kappa\simeq20$. This excess at $\kappa\simeq20$ describes an intermediate-angular-scale overdensity located in the Centaurus region.

The combination reverses in the energy interval (>40) EeV. In this bin, the most significant configuration corresponds to an overdensity at intermediate angular scales with $\kappa\simeq16$. However, an excess outside the 99.7% C.L. band is also present at $\kappa< 1$. The additional plots, analogous to the ones presented in Figure 2 and Figure 4, for the energy intervals (25–50), (32–64), and (>40) EeV are reported in appendix § 6.

In Figure 3, an interesting feature appears in the highest-energy intervals, with the presence of two spatially separated excesses, one on large angular scales and another one on intermediate angular scales, outside the 95% C.L. band. We tested whether one of the excesses arises as a statistical fluctuation of the other. To study this hypothesis, we take into consideration the energy interval ($>32$) EeV, where the most significant excess is found on large angular scales, with $\kappa= 0.16$ and $p_{\mathrm{iso}} = 3.2 \times10^{-3}$ (corresponding to $d = 0.14^{+0.04}_{-0.04}$) with center position (R.A., dec) = $(141^\circ, -44^\circ)$, and an excess on intermediate angular scales is also present with $\kappa= 13$ and $p_{\mathrm{iso}} = 6.1 \times10^{-3}$ (corresponding to $\sigma= 16.1^\circ$) and center position (R.A., dec) = $(209^\circ, -50^\circ)$. The distance between the centers of the two excesses is $\Theta= (45^{+14}_{-14})^\circ$. We perform a series of simulations of events distributed according to the likelihood present in equation (3). The parameters $(\kappa,\alpha,u_0)$ of the likelihood are selected according to the most significant excess found in the interval ($>32$) EeV, and the number of events present in the simulations follows the number of events present in ($>32$) EeV. We then maximize the TS while profiling over $\log_{10} k$, from $\log_{10} k = 0.8$ to $\log_{10} k = 2$ in steps of 0.1. We evaluate the frequency of obtaining an intermediate angular scale excess that is as significant as the intermediate-scale excess found in the interval ($>32$) EeV and as spatially distant from the large-scale excess as in the interval ($>32$) EeV. In other words, we evaluate the frequency of having $p_{\mathrm{iso}} < 6.1 \times10^{-3}$ and $\Theta/\sigma_{\Theta} > 3.3$, where $\sigma_{\Theta} = 14^\circ$ is the uncertainty on the angular separation. From simulations, the frequency of finding an excess on intermediate angular scales sufficiently significant and distant enough from the large-scale excess position observed on data is $2.8 \times10^{-3}$. Completing the analogous procedure on isotropic simulations, we find that the frequency of obtaining an excess on intermediate angular scales with enough separation is $2.6 \times10^{-3}$.

We complete a similar procedure to study whether the large-scale excess is a statistical fluctuation of the intermediate-scale one. We perform simulations considering the presence of an intermediate angular scale excess comparable to the ones observed on data, distributed according to the likelihood present in equation (3). The parameters $(\kappa,\alpha,u_0)$ of the likelihood are selected according to the most significant excess of intermediate angular scale found in the interval ($>32$) EeV. We then evaluate the frequency of obtaining an excess on large angular scales as significant as the one observed on data and enough distant from the intermediate-scale excess. The presence of a large-scale excess is derived profiling over $\log_{10} \kappa$, from $\log_{10} \kappa= -1$ to $\log_{10} \kappa= 0$ in steps of 0.1. In this case, we evaluate the frequency of having $p_{\mathrm{iso}} < 3.2 \times10^{-3}$ and $\Theta/\sigma_{\Theta} > 3.3$. We repeat the same procedure performing isotropic simulations. We find that the frequency is 0.014, while the frequency from isotropic simulations is 0.011, therefore the probability of finding an excess at large angular scale is similar whether an intermediate-scale excess is present or not. Based on these simulations, there is no evidence to suggest that the large-scale (intermediate-scale) excess above 32 EeV arises from a statistical fluctuation of the intermediate-scale (large-scale) excess.

In Figure 4, maps of the flux $\Phi$ normalized to its mean value $\langle\Phi\rangle$ are depicted for the energy intervals (25–51) EeV, (32–64) EeV, and ($>40$) EeV. To better compare the excesses observed at large and intermediate angular scales, in the left column, the flux maps are smoothed by a von Mises–Fisher distribution with concentration parameter $\kappa$ correspondent to $\sigma= 45^\circ$ ($\kappa\simeq1.5$), while in the right column we use $\sigma= 15^\circ$ ($\kappa\simeq15$). The two selections reproduce the smearing compatible with large-scale excesses $\sigma= 45^\circ$ and intermediate-scale excesses $\sigma= 15^\circ$. In each map, the black star indicates the position of the excess direction associated to the lowest $p_{\mathrm{iso}}$ obtained maximizing TS while profiling over $\kappa$. The dashed circles indicate the regions centered on the reference positions of the excesses for which the integral of the von Mises–Fisher distribution, using the parameters corresponding to the lowest value of $p_{\mathrm{iso}}$, contains 39.3% of the flux from the von Mises–Fisher distribution. This choice allows a direct comparison with the region enclosed by the $1\sigma$ contour of a two-dimensional Gaussian distribution. Areas with large values of $\Phi/\langle\Phi\rangle$ are present at the edge of the field of view in all the maps with $\sigma= 15^\circ$. These are border effects given by the low statistics of events present in those regions. In the maps with smoothing $\sigma= 45^\circ$, $\Phi/\langle\Phi\rangle$ increases with energy. This is in agreement with the increase with energy of the dipole amplitude observed on the data. In the maps with $\sigma= 15^\circ$, we observe the presence of several areas with $\Phi/\langle\Phi\rangle$ above 1. In particular, an overdensity is present in the Centaurus region for all energy intervals considered. The flux ratio of this region increases with energy, and has its maximum in the ($>40$) EeV interval.

Conclusion and outlook

With this analysis, we performed a uniform search for overdensities over different angular scales in the arrival-direction distribution of the events detected during 19 years of data acquisition of the Pierre Auger Observatory. In particular, the scope was to thoroughly characterize the most significant overdensities present in the dataset using a von Mises–Fisher distribution.

Using the standard energy intervals (4–8, 8–16, 16–32, >32) EeV defined in [1, 2], we find that the most significant overdensities correspond to large-scale excesses in all energy bins above 8 EeV. Moreover, the dipolar amplitudes of these excesses in the intervals (8–16, 16–32, >32) EeV, as well as their increase with energy, are compatible with the results of previous analyses by the Pierre Auger Collaboration.

Using a sliding window in energy, we further find that the most significant excesses are on the large angular scale, for all energy bins starting with (6–13) EeV and up to (>40) EeV, where the most significant finding is an overdensity of intermediate angular scale ($\kappa\simeq16$, $\sigma= 14.4^\circ$) located in the Centaurus region. We note that overdensities located in the Centaurus region that exceed the 95% C.L. band on intermediate angular scales are also present for the interval (32–64) and (>32) EeV.

Based on the simulations, there is no evidence that the coexistence of the two excesses at different angular scales, with characteristics similar to those observed in the data above 32 EeV, can be explained by one excess arising as a statistical fluctuation of the other. Different physical scenarios can give rise to a signal like this, including the superposition of the contributions from a local dominant source and from a population of farther away sources; as well as a mixed mass composition of cosmic rays at given energy, with particles of different charge dominantly contributing to the different excesses. Further studies, possibly including the mass composition information to be gathered from AugerPrime [11], will be very useful to shed light on these possibilities. Also the data that will be collected with the TA$\times$4 [22] upgraded detectors or with future UHECR observatories, such as GRAND [23], POEMMA [24], and GCOS [25] will be helpful for these investigations.

Maps for additional energy bins

Maps in equatorial coordinates

Acknowledgments

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

Argentina – Comisión Nacional de Energía Atómica; Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT); Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET); Gobierno de la Provincia de Mendoza; Municipalidad de Malargüe; NDM Holdings and Valle Las Leñas; in gratitude for their continuing cooperation over land access; Australia – the Australian Research Council; Belgium – Fonds de la Recherche Scientifique (FNRS); Research Foundation Flanders (FWO), Marie Curie Action of the European Union Grant No. 101107047; Brazil – Ministério da Ciência, Tecnologia e Inovação (MCTI); Czech Republic – GACR 24-13049S, CAS LQ100102401, MEYS LM2023032, 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 and CZ.02.01.01/00/22_008/0004632; France – Centre de Calcul IN2P3/CNRS; Centre National de la Recherche Scientifique (CNRS); Institut National de Physique Nucléaire et de Physique des Particules (IN2P3/CNRS); Germany – Bundesministerium für Forschung, Technologie und Raumfahrt (BMFTR); Deutsche Forschungsgemeinschaft (DFG); Ministerium für Finanzen Baden-Württemberg; Helmholtz Alliance for Astroparticle Physics (HAP); Hermann von Helmholtz-Gemeinschaft Deutscher Forschungszentren e.V.; Ministerium für Kultur und Wissenschaft des Landes Nordrhein-Westfalen; Ministerium für Wissenschaft, Forschung und Kunst des Landes Baden-Württemberg; Italy – Istituto Nazionale di Fisica Nucleare (INFN); Istituto Nazionale di Astrofisica (INAF); Ministero dell’Università e della Ricerca (MUR); CETEMPS Center of Excellence; Ministero degli Affari Esteri (MAE), ICSC Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union NextGenerationEU, reference code CN_00000013; México – Consejo Nacional de Ciencia y Tecnología (CONACYT-SECHTI) No. CB-A1-S-46703, Universidad Nacional Autónoma de México (UNAM) PAPIIT-IN114924; Benemérita Universidad Autónoma de Puebla (BUAP), VIEP and Laboratorio Nacional de Supercómputo del Sureste de México (LNS); and Benemérita Universidad Autónoma de Chiapas (UNACH); The Netherlands – Ministry of Education, Culture and Science; Netherlands Organisation for Scientific Research (NWO); Dutch national e-infrastructure with the support of SURF Cooperative; Poland – Ministry of Science and Higher Education, grant No. 2026/WK/05; National Science Centre, grant No. 2022/45/B/ST9/02163; 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 Education and Research, contract no. 30N/2023 under Romanian National Core Program LAPLAS VII, and grant no. PN 23 21 01 02; Slovenia – Slovenian Research and Innovation Agency, grants P1-0031, I0-0033; Spain – Ministerio de Ciencia, Innovación y Universidades/Agencia Estatal de Investigación MICIU/AEI /10.13039/501100011033 (PID2022-140510NB-I00, PCI2023-145952-2, CNS2024-154676, and María de Maeztu CEX2023-001318-M), Xunta de Galicia (CIGUS Network of Research Centers, Consolidación ED431C-2025/11 and ED431F-2022/15) and European Union ERDF; 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, and NSF-2013199; The Grainger Foundation; Astrophysics Centre for Multi-messenger studies in Europe (ACME) EU Grant No 10131928; and UNESCO.