Features of the Energy Spectrum of Cosmic Rays above 2.5 × 10 18 eV Using the Pierre Auger Observatory

Features of the Energy Spectrum of Cosmic Rays above 2.5×10 18 eV Using the Pierre Auger Observatory A. Aab,1P. Abreu,2M. Aglietta,3,4 J. M. Albury,5I. Allekotte,6A. Almela,7,8 J. Alvarez Castillo,9J. Alvarez-Muñiz,10 R. Alves Batista,1G. A. Anastasi,11,4 L. Anchordoqui,12 B. Andrada,7S. Andringa,2C. Aramo,13 P. R. Araújo Ferreira,14 H. Asorey,7P. Assis,2G. Avila,15,16 A. M. Badescu,17 A. Bakalova,18 A. Balaceanu,19 F. Barbato,20,13 R. J. Barreira Luz,2 K. H. Becker,21 J. A. Bellido,5C. Berat,22 M. E. Bertaina,11,4 X. Bertou,6P. L. Biermann,23 T. Bister,14 J. Biteau,24 A. Blanco,2J. Blazek,18 C. Bleve,22 M. Boháčová,18 D. Boncioli,25,26 C. Bonifazi,27 L. Bonneau Arbeletche,28 N. Borodai,29 A. M. Botti,7J. Brack,30 T. Bretz,14 F. L. Briechle,14 P. Buchholz,31 A. Bueno,32 S. Buitink,33 M. Buscemi,34,35 K. S. Caballero-Mora,36 L. Caccianiga,37,38 L. Calcagni,39 A. Cancio,8,7 F. Canfora,1,40 I. Caracas,21 J. M. Carceller,32 R. Caruso,34,35 A. Castellina,3,4 F. Catalani,41 G. Cataldi,42 L. Cazon,2M. Cerda,15 J. A. Chinellato,43 K. Choi,10 J. Chudoba,18 L. Chytka,44 R. W. Clay,5A. C. Cobos Cerutti,45 R. Colalillo,20,13 A. Coleman,46 M. R. Coluccia,47,42 R. Conceição,2A. Condorelli,48,26 G. Consolati,38,49 F. Contreras,15,16 F. Convenga,47,42 C. E. Covault,50,51 S. Dasso,52,53 K. Daumiller,54 B. R. Dawson,5J. A. Day,5R. M. de Almeida,55 J. de Jesús,7,54 S. J. de Jong,1,40 G. De Mauro,1,40 J. R. T. de Mello Neto,27,56 I. De Mitri,48,26 J. de Oliveira,55 D. de Oliveira Franco,43 V. de Souza,57 E. De Vito,47,42 J. Debatin,58 M. del Río,16 O. Deligny,24 H. Dembinski,54 N. Dhital,29 C. Di Giulio,59,60 A. Di Matteo,4M. L. Díaz Castro,43 C. Dobrigkeit,43 J. C. D’Olivo,9Q. Dorosti,31 R. C. dos Anjos,61 M. T. Dova,39 J. Ebr,18 R. Engel,58,54 I. Epicoco,47,42 M. Erdmann,14 C. O. Escobar,62 A. Etchegoyen,7,8 H. Falcke,1,63,40 J. Farmer,64 G. Farrar,65 A. C. Fauth,43 N. Fazzini,62 F. Feldbusch,66 F. Fenu,11,4 B. Fick,67 J. M. Figueira,7A. Filipčič,68,69 T. Fodran,1M. M. Freire,70 T. Fujii,64,71 A. Fuster,7,8 C. Galea,1C. Galelli,37,38 B. García,45 A. L. Garcia Vegas,14 H. Gemmeke,66 F. Gesualdi,7,54 A. Gherghel-Lascu,19 P. L. Ghia,24 U. Giaccari,1M. Giammarchi,38 M. Giller,72 J. Glombitza,14 F. Gobbi,15 F. Gollan,7G. Golup,6 M. Gómez Berisso,6P. F. Gómez Vitale,15,16 J. P. Gongora,15 N. González,7I. Goos,6,54 D. Góra,29 A. Gorgi,3,4 M. Gottowik,21 T. D. Grubb,5F. Guarino,20,13 G. P. Guedes,73 E. Guido,4,11 S. Hahn,54,7 R. Halliday,50 M. R. Hampel,7 P. Hansen,39 D. Harari,6V. M. Harvey,5A. Haungs,54 T. Hebbeker,14 D. Heck,54 G. C. Hill,5C. Hojvat,62 J. R. Hörandel,1,40 P. Horvath,44 M. Hrabovský,44 T. Huege,54,33 J. Hulsman,7,54 A. Insolia,34,35 P. G. Isar,74 J. A. Johnsen,75 J. Jurysek,18 A. Kääpä,21 K. H. Kampert,21 B. Keilhauer,54 J. Kemp,14 H. O. Klages,54 M. Kleifges,66 J. Kleinfeller,15 M. Köpke,58 G. Kukec Mezek,69 B. L. Lago,76 D. LaHurd,50 R. G. Lang,57 M. A. Leigui de Oliveira,77 V. Lenok,54 A. Letessier-Selvon,78 I. Lhenry-Yvon,24 D. Lo Presti,34,35 L. Lopes,2R. López,79 R. Lorek,50 Q. Luce,58 A. Lucero,7A. Machado Payeras,43 M. Malacari,64 G. Mancarella,47,42 D. Mandat,18 B. C. Manning,5J. Manshanden,80 P. Mantsch,62 S. Marafico,24 A. G. Mariazzi,39 I. C. Mariş,81 G. Marsella,47,42 D. Martello,47,42 H. Martinez,57 O. Martínez Bravo,79 M. Mastrodicasa,25,26 H. J. Mathes,54 J. Matthews,82 G. Matthiae,59,60 E. Mayotte,21 P. O. Mazur,62 G. Medina-Tanco,9D. Melo,7A. Menshikov,66 K.-D. Merenda,75 S. Michal,44 M. I. Micheletti,70 L. Miramonti,37,38 D. Mockler,81 S. Mollerach,6F. Montanet,22 C. Morello,3,4 M. Mostafá,83 A. L. Müller,7,54 M. A. Muller,43,84,27 K. Mulrey,33 R. Mussa,4M. Muzio,65 W. M. Namasaka,21 L. Nellen,9P. H. Nguyen,5M. Niculescu-Oglinzanu,19 M. Niechciol,31 D. Nitz,67,85 D. Nosek,86 V. Novotny,86 L. Nožka,44 A. Nucita,47,42 L. A. Núñez,87 M. Palatka,18 J. Pallotta,88 M. P. Panetta,47,42 P. Papenbreer,21 G. Parente,10 A. Parra,79 M. Pech,18 F. Pedreira,10 J. P¸ekala,29 R. Pelayo,89 J. Peña-Rodriguez,87 J. Perez Armand,28 M. Perlin,7,54 L. Perrone,47,42 C. Peters,14 S. Petrera,48,26 T. Pierog,54 M. Pimenta,2V. Pirronello,34,35 M. Platino,7B. Pont,1M. Pothast,40,1 P. Privitera,64 M. Prouza,18 A. Puyleart,67 S. Querchfeld,21 J. Rautenberg,21 D. Ravignani,7M. Reininghaus,54,7 J. Ridky,18 F. Riehn,2 M. Risse,31 P. Ristori,88 V. Rizi,25,26 W. Rodrigues de Carvalho,28 G. Rodriguez Fernandez,59,60 J. Rodriguez Rojo,15 M. J. Roncoroni,7M. Roth,54 E. Roulet,6A. C. Rovero,52 P. Ruehl,31 S. J. Saffi,5A. Saftoiu,19 F. Salamida,25,26 H. Salazar,79 G. Salina,60 J. D. Sanabria Gomez,87 F. Sánchez,7E. M. Santos,28 E. Santos,18 F. Sarazin,75 R. Sarmento,2 C. Sarmiento-Cano,7R. Sato,15 P. Savina,47,42,24 C. Schäfer,54 V. Scherini,42 H. Schieler,54 M. Schimassek,58,7 M. Schimp,21 F. Schlüter,54,7 D. Schmidt,58 O. Scholten,90,33 P. Schovánek,18 F. G. Schröder,46,54 S. Schröder,21 A. Schulz,54 S. J. Sciutto,39 M. Scornavacche,7,54 R. C. Shellard,91 G. Sigl,80 G. Silli,7,54 O. Sima,19,51 R. Šmída,64 P. Sommers,83 J. F. Soriano,12 J. Souchard,22 R. Squartini,15 M. Stadelmaier,54,7 D. Stanca,19 S. Stanič,69 J. Stasielak,29 P. Stassi,22 A. Streich,58,7 M. Suárez-Durán,87 T. Sudholz,5T. Suomijärvi,24 A. D. Supanitsky,7J. Šupík,44 Z. Szadkowski,92 A. Taboada,58 A. Tapia,93 C. Timmermans,40,1 O. Tkachenko,54 P. Tobiska,18 C. J. Todero Peixoto,41 B. Tom´e,2G. Torralba Elipe,10 A. Travaini,15 P. Travnicek,18 C. Trimarelli,25,26 M. Trini,69 M. Tueros,39 R. Ulrich,54 M. Unger,54 M. Urban,14 L. Vaclavek,44 M. Vacula,44 J. F. Vald´es Galicia,9I. Valiño,48,26 L. Valore,20,13 A. van Vliet,1E. Varela,79 B. Vargas Cárdenas,9A. Vásquez-Ramírez,87 PHYSICAL REVIEW LETTERS 125, 121106 (2020) Editors' Suggestion Featured in Physics 0031-9007=20=125(12)=121106(10) 121106-1 Published by the American Physical Society D. Veberič,54 C. Ventura,56 I. D. Vergara Quispe,39 V. Verzi,60 J. Vicha,18 L. Villaseñor,79 J. Vink,94 S. Vorobiov,69 H. Wahlberg,39 A. A. Watson,95 M. Weber,66 A. Weindl,54 L. Wiencke,75 H. Wilczyński,29 T. Winchen,33 M. Wirtz,14 D. Wittkowski,21 B. Wundheiler,7A. Yushkov,18 O. Zapparrata,81 E. Zas,10 D. Zavrtanik,69,68 M. Zavrtanik,68,69 L. Zehrer,69 A. Zepeda,96 M. Ziolkowski,31 and F. Zuccarello34,35 (The Pierre Auger Collaboration)* 1IMAPP, Radboud University Nijmegen, Nijmegen, The Netherlands 2Laboratório de Instrumentação e Física Experimental de Partículas—LIP and Instituto Superior T´ecnico—IST, Universidade de Lisboa—UL, Lisboa, Portugal 3Osservatorio Astrofisico di Torino (INAF), Torino, Italy 4INFN, Sezione di Torino, Torino, Italy 5University of Adelaide, Adelaide, S.A., Australia 6Centro Atómico Bariloche and Instituto Balseiro (CNEA-UNCuyo-CONICET), San Carlos de Bariloche, Argentina 7Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), Buenos Aires, Argentina 8Universidad Tecnológica Nacional—Facultad Regional Buenos Aires, Buenos Aires, Argentina 9Universidad Nacional Autónoma de M´exico, M´exico, D.F., M´exico 10Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de Santiago de Compostela, Santiago de Compostela, Spain 11Universit`a Torino, Dipartimento di Fisica, Torino, Italy 12Department of Physics and Astronomy, Lehman College, City University of New York, Bronx, NY, USA 13INFN, Sezione di Napoli, Napoli, Italy 14RWTH Aachen University, III. Physikalisches Institut A, Aachen, Germany 15Observatorio Pierre Auger, Malargüe, Argentina 16Observatorio Pierre Auger and Comisión Nacional de Energía Atómica, Malargüe, Argentina 17University Politehnica of Bucharest, Bucharest, Romania 18Institute of Physics of the Czech Academy of Sciences, Prague, Czech Republic 19“Horia Hulubei”National Institute for Physics and Nuclear Engineering, Bucharest-Magurele, Romania 20Universit`a di Napoli “Federico II”, Dipartimento di Fisica “Ettore Pancini”, Napoli, Italy 21Bergische Universität Wuppertal, Department of Physics, Wuppertal, Germany 22Universit´e Grenoble Alpes, CNRS, Grenoble Institute of Engineering Universit´e Grenoble Alpes, LPSC-IN2P3, 38000 Grenoble, France, France 23Max-Planck-Institut für Radioastronomie, Bonn, Germany 24Universit´e Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France, France 25Universit`a dell’Aquila, Dipartimento di Scienze Fisiche e Chimiche, L’Aquila, Italy 26INFN Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy 27Universidade Federal do Rio de Janeiro, Instituto de Física, Rio de Janeiro, RJ, Brazil 28Universidade de São Paulo, Instituto de Física, São Paulo, SP, Brazil 29Institute of Nuclear Physics PAN, Krakow, Poland 30Colorado State University, Fort Collins, CO, USA 31Universität Siegen, Fachbereich 7 Physik—Experimentelle Teilchenphysik, Siegen, Germany 32Universidad de Granada and C.A.F.P.E., Granada, Spain 33Vrije Universiteit Brussels, Brussels, Belgium 34Universit`a di Catania, Dipartimento di Fisica e Astronomia, Catania, Italy 35INFN, Sezione di Catania, Catania, Italy 36Universidad Autónoma de Chiapas, Tuxtla Guti´errez, Chiapas, M´exico 37Universit`a di Milano, Dipartimento di Fisica, Milano, Italy 38INFN, Sezione di Milano, Milano, Italy 39IFLP, Universidad Nacional de La Plata and CONICET, La Plata, Argentina 40Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF), Science Park, Amsterdam, The Netherlands 41Universidade de São Paulo, Escola de Engenharia de Lorena, Lorena, SP, Brazil 42INFN, Sezione di Lecce, Lecce, Italy 43Universidade Estadual de Campinas, IFGW, Campinas, SP, Brazil 44Palacky University, RCPTM, Olomouc, Czech Republic 45Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), and Universidad Tecnológica Nacional—Facultad Regional Mendoza (CONICET/CNEA), Mendoza, Argentina 46University of Delaware, Department of Physics and Astronomy, Bartol Research Institute, Newark, DE, USA 47Universit`a del Salento, Dipartimento di Matematica e Fisica “E. De Giorgi”, Lecce, Italy 48Gran Sasso Science Institute, L’Aquila, Italy 49Politecnico di Milano, Dipartimento di Scienze e Tecnologie Aerospaziali, Milano, Italy PHYSICAL REVIEW LETTERS 125, 121106 (2020) 121106-2 50Case Western Reserve University, Cleveland, OH, USA 51also at Radboud Universtiy Nijmegen, Nijmegen, The Netherlands 52Instituto de Astronomía y Física del Espacio (IAFE, CONICET-UBA), Buenos Aires, Argentina 53Departamento de Física and Departamento de Ciencias de la Atmósfera y los Oc´eanos, FCEyN, Universidad de Buenos Aires and CONICET, Buenos Aires, Argentina 54Karlsruhe Institute of Technology, Institut für Kernphysik, Karlsruhe, Germany 55Universidade Federal Fluminense, EEIMVR, Volta Redonda, RJ, Brazil 56Universidade Federal do Rio de Janeiro (UFRJ), Observatório do Valongo, Rio de Janeiro, RJ, Brazil 57Universidade de São Paulo, Instituto de Física de São Carlos, São Carlos, SP, Brazil 58Karlsruhe Institute of Technology, Institute for Experimental Particle Physics (ETP), Karlsruhe, Germany 59Universit`a di Roma “Tor Vergata”, Dipartimento di Fisica, Roma, Italy 60INFN, Sezione di Roma “Tor Vergata”, Roma, Italy 61Universidade Federal do Paraná, Setor Palotina, Palotina, Brazil 62Fermi National Accelerator Laboratory, USA 63Stichting Astronomisch Onderzoek in Nederland (ASTRON), Dwingeloo, The Netherlands 64University of Chicago, Enrico Fermi Institute, Chicago, IL, USA 65New York University, New York, NY, USA 66Karlsruhe Institute of Technology, Institut für Prozessdatenverarbeitung und Elektronik, Karlsruhe, Germany 67Michigan Technological University, Houghton, MI, USA 68Experimental Particle Physics Department, J. Stefan Institute, Ljubljana, Slovenia 69Center for Astrophysics and Cosmology (CAC), University of Nova Gorica, Nova Gorica, Slovenia 70Instituto de Física de Rosario (IFIR)—CONICET/U.N.R. and Facultad de Ciencias Bioquímicas y Farmac´euticas U.N.R., Rosario, Argentina 71now at Hakubi Center for Advanced Research and Graduate School of Science, Kyoto University, Kyoto, Japan 72University of Łódź, Faculty of Astrophysics, Łódź, Poland 73Universidade Estadual de Feira de Santana, Feira de Santana, Brazil 74Institute of Space Science, Bucharest-Magurele, Romania 75Colorado School of Mines, Golden, CO, USA 76Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Nova Friburgo, Brazil 77Universidade Federal do ABC, Santo Andr´e, SP, Brazil 78Laboratoire de Physique Nucl´eaire et de Hautes Energies (LPNHE), Universit´es Paris 6 et Paris 7, CNRS-IN2P3, Paris, France 79Benem´erita Universidad Autónoma de Puebla, Puebla, M´exico 80Universität Hamburg, II. Institut für Theoretische Physik, Hamburg, Germany 81Universit´e Libre de Bruxelles (ULB), Brussels, Belgium 82Louisiana State University, Baton Rouge, LA, USA 83Pennsylvania State University, University Park, PA, USA 84Universidade Federal de Alfenas, Poços de Caldas, Brazil 85Karlsruhe Institute of Technology, Karlsruhe, Germany 86Charles University, Faculty of Mathematics and Physics, Institute of Particle and Nuclear Physics, Prague, Czech Republic 87Universidad Industrial de Santander, Bucaramanga, Colombia 88Centro de Investigaciones en Láseres y Aplicaciones, CITEDEF and CONICET, Villa Martelli, Argentina 89Unidad Profesional Interdisciplinaria en Ingeniería y Tecnologías Avanzadas del Instituto Polit´ecnico Nacional (UPIITA-IPN), M´exico, D.F., M´exico 90KVI—Center for Advanced Radiation Technology, University of Groningen, Groningen, The Netherlands 91Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, RJ, Brazil 92University of Łódź, Faculty of High-Energy Astrophysics,Łódź, Poland 93Universidad de Medellín, Medellín, Colombia 94Universiteit van Amsterdam, Faculty of Science, Amsterdam, The Netherlands 95School of Physics and Astronomy, University of Leeds, Leeds, United Kingdom 96Centro de Investigación y de Estudios Avanzados del IPN (CINVESTAV), M´exico, D.F., M´exico (Received 20 May 2020; accepted 31 July 2020; published 16 September 2020) We report a measurement of the energy spectrum of cosmic rays above 2.5×1018 eV based on 215 030 events. New results are presented: at about 1.3×1019 eV, the spectral index changes from Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. PHYSICAL REVIEW LETTERS 125, 121106 (2020) 121106-3 2.51 0.03ðstatÞ0.05ðsystÞto 3.05 0.05ðstatÞ0.10ðsystÞ, evolving to 5.10.3ðstatÞ0.1ðsystÞ beyond 5×1019 eV, while no significant dependence of spectral features on the declination is seen in the accessible range. These features of the spectrum can be reproduced in models with energy-dependent mass composition. The energy density in cosmic rays above 5×1018 eV is ½5.66 0.03ðstatÞ1.40ðsystÞ × 1053 erg Mpc−3. DOI: 10.1103/PhysRevLett.125.121106 Although cosmic rays having energies above 1019 eV were first detected nearly 60 years ago [1,2] and are being investigated by the two largest-ever built detectors, the Pierre Auger Observatory [3] and the Telescope Array [4], the question of their origin remains unanswered. Only recently has the belief that such particles are of extragalactic origin been demonstrated experimentally with the discovery of significant directional anisotropies above 8×1018 eV [5]. These data are well described by a dipole pattern, the amplitude of which increases from ≃6% to ≃10% as the energy rises to ≃4×1019 eV [6]. An important observable for an understanding of ultrahigh energy cosmic rays (UHECRs) is the energy spectrum. We report a measurement above 2.5×1018 eV based on 215 030 events, over 10 times that used in [7]. Over 16 000 events have energies beyond 1019 eV. This spectral determination is unique in making no assumptions about the mass composition or the hadronic physics. Full details are reported in [8]. UHECRs can only be studied through the detection of the showers of particles (extensiveair showers) they createin the atmosphere. A calorimetric estimate of the energy carried by the primary particle is possible using telescopes to collect the fluorescence light emitted by atmospheric nitrogen excited by the shower. The on time of this technique is limited to nights with low-background light while, by contrast, an array of particle detectors deployed on the ground can be operated with a duty cycle close to 100%. The traditional method of assessing the energy of the primary cosmic ray from observations made with the particle detectors requires assumptions about its mass and the hadronic processes that control the cascade development. This is clearly unsatisfactory as the mass is unknown and the center-of-mass energy reached at the LHC corresponds only to that of a proton of ≃1017 eV colliding with a nitrogen nucleus. Also, details of pion interactions, which play a key role in shower development, are lacking. The presence of unknown processes could also lead to hidden systematic uncertainties. To circumvent these limitations, the energies are obtained by making use of a subset of events detected simultaneously by the fluorescence detector (FD) and the particle surface detectors (SD). This “hybrid”approach allows a calorimetric estimate of the energy for events recorded during periods when the FD cannot be operated. A spectrum can thus be derived that is free from assumptions about primary mass or hadronic physics. The Pierre Auger Observatory is such a hybrid system [3]. It is sited near the city of Malargüe, Argentina, at latitude 35.2° S with a mean atmospheric overburden of 875 g cm−2[3]. The SD comprises 1600 water-Cherenkov detectors deployed on a 1500 m triangular grid, covering about 3000 km2. The array is overlooked from four stations, each containing six telescopes used to detect the emitted fluorescence light. Comprehensive atmospheric monitoring, particularly of the aerosol content and the cloud cover, is undertaken [3,9]. The SD samples the shower particles that reach the ground. Signals in the individual detectors are quantified in terms of their response to a muon traveling vertically and centrally through it (a vertical equivalent muon or VEM). The signals are used to determine the impact point of the shower axis, the arrival direction, and the shower size. For the latter, the signal at 1000 m from the shower axis, Sð1000Þ, is used. For the grid spacing of 1500 m, this is the distance that minimizes the uncertainty arising in Sð1000Þ from the imperfect knowledge of the functional form describing the falloff of signal with distance from the shower axis in individual events [10]. Showers detected by the SD arrive from a range of zenith angles, and they are attenuated according to how much atmosphere is traversed. Accordingly, for each event, Sð1000Þis adjusted to a reference value S38, the magnitude that it would have had, had the cosmic ray arrived at the median zenith angle of 38°. The long-established procedure for making this correction, the Constant Intensity Cut [11] method, relies on the quasi-isotropy of cosmic rays in zenith angle given the small anisotropy contrasts in celestial coordinates [8]. The large number of events has made it possible to refine the original approach and quantify the change in shower absorption as a function of energy. Such an evolution is anticipated as, at a given zenith angle, the ratio of the muon to electromagnetic components falls as the energy increases, even for an energy-independent composition. For showers detected by the FD, it is possible to measure the deposition of energy lost to ionization of the atmosphere using a fit to a modified Gaisser-Hillas profile [12]. The integration of the profile yields a calorimetric measure of this loss. The energy of the primary particle EFD is then obtained by the addition of an energy-dependent correction of less than 14%, driven by data [13], to allow for the “invisible energy,”carried into the ground by muons and PHYSICAL REVIEW LETTERS 125, 121106 (2020) 121106-4 neutrinos. The resolution in EFD is well described by σFDðEÞ=E ≃7.4% over the whole energy range [14]. Hybrid events are thus used to develop a calibration curve such that every estimate of S38 can be assigned a valuation of EFD. Here 3338 hybrid events surviving rigorous quality cuts [8] are used to obtain a relationship between S38 and EFD of the form EFD ¼ASB 38, where A¼ð1.860.03Þ×1017 eV and B¼1.0310.004.No zenithal dependence of Aor Bhas been found, further validating the use of the constant intensity method [8]. Such a simple dependence is sufficient to describe the data in full detail. The energies of the hybrid events range from 2.5×1018 eV to 8×1019 eV. The most energetic event, detected at all fluorescence stations, has an energy EFD ¼ð8.50.4Þ×1019 eV, derived from a weighted average of the four independent estimates of the calorimetric energy. For this event S38 ¼354 VEM so that the energy deduced from the calibration curve is ESD ≡ ASB 38 ¼ð7.90.6Þ×1019 eV. The systematic uncertainty in the energy assignment is about 14% over the whole energy range [15]. This benefits from the high-precision AIRFLY Collaboration measurement of the fluorescence yield [16] and from an accurate data-driven estimation of the invisible energy [13]. Other contributions to the uncertainty are related to the estimation of the Aand B parameters, the characterization of the atmosphere, the reconstruction of the longitudinal profile and the FD calibration, which provides the largest contribution. To derive the energy spectrum, we use events recorded by the SD with the largest-signal station not located on the boundary of the array, with zenith angle θ<60° and energy ≥2.5×1018 eV. These selection criteria not only ensure adequate sampling of the shower but also allow the evaluation of the aperture of the SD in a purely geometrical manner in the regime where the array trigger is fully efficient and independent of the mass or energy of the primary particle [17]. The resulting SD dataset consists of 215030 events recorded between 1 January 2004 and 31 August 2018, from an exposure Eof ð60 400 1810Þkm2sr yr. The determination of E, dependent only on the acceptance angle, the surface area and the live time of the array, is discussed in detail in [17]. The procedure for extracting the spectrum from the observations, fully discussed in [8], is summarized here. The energy spectrum, typically a power law (∝E−γ) with spectral index γin a given energy interval, is estimated as Ji¼ciNi=ðEΔEiÞ, with Nithe number of observed events in differential bins of width Δlog10 Ei¼0.1and cithe correction factors required to eliminate the biases caused by the finite energy resolution. The size of the bins is such that it corresponds approximately to the energy resolution in the lowest energy bin, which starts at 2.5×1018 eV. The correction factors are needed because, as the spectrum is steep, the finite resolution causes migration between bins, particularly from lower to higher energies, artificially enhancing the flux. At the lowest energies, the correction depends also on the behavior of the detection efficiency in the energy region where the array is not fully efficient as well as on the bias in the energy due to triggerselection effects. A forward-folding approach is used to determine the correction factors. It consists of finding the model of the energy spectrum folded for detector effects that best describes the data, and then using this model to calculate the values of ci. The SD efficiency can be estimated from the fraction of hybrid events that also satisfy the SD trigger conditions, because above 1018 eV, the hybrid trigger efficiency is 100% independent of primary mass [18]. The energy resolution of ESD, and the bias in its estimate, are found from a study of the distributions of ESD=EFD. The resolution improves from ≈20% at 2×1018 eV to ≈7% at 2×1019 eV and is constant thereafter. The bias is zero above 2.5×1018 eV and increases smoothly going to lower energies and larger zenith angles: at 1018 eV it is ≈10% at 0° and ≈30% at 60°. Thanks to the hybrid measurements, the correction factors are estimated avoiding any reliance on model and primary mass assumptions. The factors are maximal at the lowest energies, ≈8%, and less than 5% at the highest energies available. Further details are given in [8]. The model of the energy spectrum that we used for over a decade is a series of two power laws followed by a slow suppression. With the current exposure, this model turns out to describe the data poorly, as the reduced deviance is found to be 35.6=15 [8]. Consequently, we adopt a more complex function with a sequence of four power laws with smooth transitions [19], JðEÞ¼J0E 1018.5eV−γ1Y 3 i¼11þE Eij1=ωij ðγi−γjÞωij ; with j¼iþ1and ωij ¼0.05. The ωij factors control the widths of the energy intervals over which the slope transitions occur [8]. This model describes the data with a reduced deviance 17.0=12, which allows us to disfavor the previous parametrization with 3.9σconfidence [8]. The resulting differential energy spectrum and the fitted function are shown in Fig. 1. The normalization is J0¼ð1.315 0.004 0.400Þ×10−18 km−2sr−1yr−1eV−1. The ankle is described by a rollover at E12 ¼ ð5.00.10.8Þ×1018 eV, marking a hardening of the spectrum from γ1¼3.29 0.02 0.10 to γ2¼2.51 0.03 0.05.AtE23 ¼ð13 12Þ×1018 eV, the spectrum softens from γ2to γ3¼3.05 0.05 0.10. Finally, the spectrum softens further above a suppression energy of E34 ¼ð46 36Þ×1018 eV with γ4¼5.10.30.1, confirming with higher precision previous reports of the strong attenuation of the flux at the highest energies [7,20,21]. The feature at E23, calling for a two-step PHYSICAL REVIEW LETTERS 125, 121106 (2020) 121106-5 suppression, is a new observation. For all parameters and observables presented in the text, the first error is statistical and the second systematic. From the measured energy spectrum one can infer the differential energy density per dex (dex indicates decade in log10 E, following the convention of [22]), obtained as lnð10Þð4π=cÞE2JðEÞ. It provides a measurement of the energy density of the local Universe attributable to cosmic rays. Above the ankle, a range in which UHECRs are of extragalactic origin [5], the integration over energy results in ð5.66 0.03 1.40Þ×1053 erg Mpc−3. This translates into constraints on the luminosity of the sources, as discussed below. A detailed examination of the systematic uncertainties of the energy spectrum is reported in [8]. The uncertainty in the flux amounts to 30%–40% near 2.5×1018 eV, 25% at 1019 eV, and 60% at the highest energies. The uncertainties include contributions from the absolute energy scale (the largest), the exposure, the unfolding procedure, and the Sð1000Þreconstruction. No indication of further systematic uncertainties has been found from a comparison of the spectra calculated over different time periods, seasons, and ranges of zenith angle. The wide declination range covered, from δ¼−90°to δ¼þ24.8°, allows a search for dependencies of energy spectra on declination. For this, we have divided the sky into three declination bands of equal exposure. In each band, the estimation of the spectrum is made as for the whole field of view, but using unfolding-correction factors relevant to the band in question. We report in Table Ithe parameters characterizing the spectral features for each declination range. They are seen to be in statistical agreement. There is thus no obvious dependence with declination over the energy range covered. A trend for the intensity to be slightly higher in the Southern Hemisphere is observed [8], consistent with the anisotropy observations [6]. We therefore claim a second new result, namely that the energy spectrum does not vary as a function of declination in the range accessible at the Auger Observatory other than in the mild excess from the Southern Hemisphere expected in line with the known energy-dependent anisotropies above 8×1018 eV. A comparison of the spectrum with that of Telescope Array measured in the Northern Hemisphere is discussed in [8] and references therein. Astrophysical implications of the features of the energy spectrum.—We now examine the validity of models proposed to explain features of UHECRs using the new information given here and the data on mass composition and arrival directions recently reported [5,6,23–28].If UHECRs are produced throughout the Universe, to reach Earth they must cross the background photon fields permeating the extragalactic space. In particular, the cosmic microwave background photons induce pion production with protons colliding at around 5×1019 eV and photodisintegration of heavier nuclei at a roughly similar threshold, leading to the expectation of a spectral steepening (the Greisen-Zatsepin-Kuzmin (GZK) effect [29]). Depending on the energy and chemical composition of the UHECRs, higher-energy background photons, such as infrared light, may also be responsible of interactions producing the flux steepening. A popular framework has been that what is observed comes from universal sources, uniformly distributed, that [eV]E 19 10 20 10 ] 1 eV -1 sr -1 yr -2 [km 2 E×)E(J 17 10 18 10 19 10 ] -1 dex -3 Energy density [erg Mpc 52 10 53 10 54 10 83143 47500 28657 17843 12435 8715 6050 4111 2620 1691 991 624 372 156 83 24 9 6 [eV]E 19 10 20 10 ] 2 eV -1 sr -1 yr -2 [km 3 E×)E(J 37 10 sys σ ± stat σ ±spectral index 0.10± 0.02 ± = 3.29 1 γ 0.05± 0.03 ± = 2.51 2 γ 0.10± 0.05 ± = 3.05 3 γ 0.1± 0.3 ± = 5.1 4 γ 1 2 3 4 FIG. 1. Top: energy spectrum scaled by E2with the number of detected events in each energy bin. In this representation the data provide an estimation of the differential energy density per decade. Bottom: energy spectrum scaled by E3fitted with a sequence of four power laws (red line). The numbers (i¼1;…;4) enclosed in the circles identify the energy intervals where the spectrum is described by a power law with spectral index γi. The shaded band indicates the statistical uncertainty of the fit. Upper limits are at the 90% confidence level. TABLE I. Spectral parameters in three different declination ranges. The energies E12,E23, and E34 are given in units of 1018 eV and the normalization parameter J0in units of 1018 km−2sr−1yr−1eV−1. Uncertainties are statistical. ½−90.0°;−42.5°½−42.5°;−17.3°½−17.3°;þ24.8° J01.329 0.007 1.306 0.007 1.312 0.006 γ13.26 0.03 3.31 0.03 3.30 0.03 γ22.53 0.04 2.54 0.04 2.44 0.05 γ33.10.13.00.13.00.1 γ45.20.44.40.35.70.6 E12 5.10.24.90.25.20.2 E23 14 2142121 E34 47 4374514 PHYSICAL REVIEW LETTERS 125, 121106 (2020) 121106-6 accelerate only protons. As a consequence the ankle is then explicable by energy losses of protons through pair production across greater distances [30–32] so that the ankle region would be proton dominated. However, recent results [26] strongly contradict this expectation: in the ankle region, ð3–5Þ×1018 eV, a pure proton composition, or one of only protons and helium, is excluded at the 6.4σlevel. A second consequence [33] concerns the energy E1=2, at which the integral intensity falls by a factor of 2 with respect to a power-law extrapolation from lower energies. The prediction in this framework is that E1=2¼5.3×1019 eV, though this number may change with fluctuations of source luminosities and densities that shape the GZK feature [31,34], and with the maximum achievable energy in the sources. The value found here, ð2.20.10.3Þ×1019 eV, is at variance withthe prediction because of the new feature of the spectrum at ≈1019 eV, which is absent in the popular paradigm that is thus disfavored. Relaxing the universality of the source spectra, the steepening at ≈1019 eV could stem from the distinctive spectrum of a local source that emits protons and contributes significantly to the total intensity. At these energies, diffusive propagation of protons from a nearby source is excluded by limits set on extragalactic magnetic fields from rotation measures [35]. Approximately, protons would thus arrive to the Galaxy as a uniform, parallel beam that may subsequently be focused or defocused while propagating in the Galactic magnetic field. As seen from the Earth, the image of the source is expected to be shifted and broadened, with the effect growing with decreasing energy. Also, multiple broad images may be produced if uncorrelated regions of the magnetic field are experienced by the particles [36–38]. Such a scenario would thus imply the observation of an anisotropy at intermediate angular scales, the size of which depends on the model of turbulence for the magnetic field [39]. Spectral differences would also consequently appear in some parts of the sky. The softening at ≈1019 eV, in particular, would not be expected in every declination range. The absence of such dependence accordingly disfavors the interpretation that the steepening is due to a source in the local Universe emitting protons. Furthermore, the interplay between the luminosity of a given source and its flux attenuation with distance requires fine-tuning to make viable a scenario in which several sources would emit protons with a distinctive spectrum while at the same time no directional effect would be imprinted upon the observed intensity. By contrast, our results fit a scenario in which several nuclear components contribute to the total intensity and in which the electromagnetic fields permeate source environments where nuclei are accelerated to a maximum energy proportional to their charge (Z). This scenario, e.g., [40–43], provides a natural framework to explain the tendency toward heavier masses with increasing energy as inferred from recent works [23–25]. To illustrate the main physics aspects without distraction by the many details a full model scenario would require, we consider here, as in [43], several nuclear components injected at the sources with a power-law spectrum and with the maximal energy of the sources modeled with an exponential cutoff. For simplicity, the sources are assumed to be stationary and uniform in a comoving volume. We show in Fig. 2the best reproduction of the data by simultaneously fitting the energy spectrum above 5×1018 eV and the distribution of the depths of the shower maximum (Xmax), which is mass sensitive (using EPOS LHC [44] as a model of hadronic interactions in their interpretation). The abundance of nuclear elements at the sources is dominated by intermediate-mass nuclei accelerated to ≈5Z×1018 eV and escaping from the source environments with a very hard spectrum. In this scenario, the steepening observed above ≈5×1019 eV results from the combination of the maximum energy of acceleration of the heaviest nuclei at the sources and the GZK effect. The steepening at ≈1019 eV reflects the interplay between the flux contributions of the helium and carbon-nitrogen-oxygen components injected at the source with their distinct cutoff energies, shaped by photodisintegration during the propagation. We note that the ratio E34=E23 is 3.40.3, matching the mass-to-charge ratio of CNO to He, as expected from the benchmark scenario shown in Fig. 2. Some cautionary comments on the illustrative model considered here are in order. The presence of a subdominant light component at the highest energies is not excluded by our data, see, e.g., [45]. Also, viable source scenarios can be found without resorting to a mixed composition with a rigidity-dependent maximum energy if, for instance, predominately heavy (Si to Fe) nuclei are accelerated and photodisintegrate in the source environment [46] or en route to Earth [47,48]. Scenarios with a predominantly 19 10 20 10 [eV]E 51 10 52 10 53 10 54 10 ] -1 dex -3 Energy density [erg Mpc Total = 1A 4≤A≤2 22≤A≤5 38≤A≤23 56≤A≤39 FIG. 2. Energy density obtained with the best fit parameters of the benchmark scenario used for illustration, as described in the text. The dashed curve shows the energy range that is not used in the fit and where an additional component is needed for describing the spectrum. PHYSICAL REVIEW LETTERS 125, 121106 (2020) 121106-7 light composition [32,49] can fit our Xmax data as well as those of Telescope Array [50] in the ankle energy range, but these scenarios are at odds with measurements of the correlation of particle densities at ground and Xmax [26]. At ultrahigh energies, a significant readjustment of current hadronic interaction models would be required [51] to fit our data with a p=He-dominated model while the data of Telescope Array, because of limited statistical power above 1019 eV, cannot yet be used to draw reliable conclusions about composition in this energy range [52]. Interactions of the accelerated nuclei in the environment of the sources may give rise to copious fluxes of nucleons below the ankle energy, produced through photodisintegration. Neutrons escaping from the magnetic confinement regions may then explain the observed flux of protons deduced from Xmax measurements [24,25] in this energy range, due to neutron decay during propagation [46,53–57]. To make up the all-particle spectrum and to fit the composition data below ≈5×1018 eV, an additional component is further required (see, e.g., [58–60] for discussions). This could be the high-energy tail from the sources emitting the bulk of Galactic cosmic rays of lower energies or, as in the “B-component scenario”[61], further explored in [62,63], this additional high-energy component is produced by different sources in the Galaxy. Finally, within this scenario, our data constrain the luminosity density that continuously emitting sources must inject into extragalactic space in UHECRs to supply the observed energy density. This amounts to ≈6×1044 erg Mpc−3yr−1above 5×1018 eV at a redshift of zero, in line with the value of 2×1044 erg Mpc−3yr−1 that can be inferred dividing the measured energy density by the typical cosmic-ray energy loss time Oð1Gpc=cÞ (3.3 Gyr) [64]. Classes of extragalactic sources that match such rates in the gamma-ray band include active galactic nuclei and starburst galaxies [65]. The flux pattern from these objects also provides an indication of anisotropy in UHECR arrival directions [27,28]. 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 yT´ecnicas (CONICET); Gobierno de la Provincia de Mendoza; Municipalidad de Malargüe; NDM Holdings and Valle Las Leñas; in gratitude for their continuing cooperation over land access; Australia—the Australian Research Council; Brazil—Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq); Financiadora de Estudos e Projetos (FINEP); Fundação de Amparo `a 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´erio 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´egional Ile-de-France; D´epartement Physique Nucl´eaire et Corpusculaire (PNC-IN2P3/ CNRS); D´epartement 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´exico—Consejo Nacional de Ciencia y Tecnología (CONACYT) No. 167733; Universidad Nacional Autónoma de M´exico (UNAM); PAPIIT DGAPA-UNAM; Netherlands—Ministry of Education, Culture and Science; Netherlands Organisation for Scientific Research (NWO); Dutch national e-infrastructure with the support of SURF Cooperative; Poland—Ministry of Science and Higher Education, Grant No. DIR/WK/2018/11; National Science Centre, Grants No. 2013/08/M/ST9/00322, No. 2016/23/B/ST9/01635, and No. HARMONIA 5–2013/10/M/ST9/00062, UMO-2016/22/M/ST9/00198; Portugal—Portuguese national funds and FEDER funds within Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia (COMPETE); Romania—Romanian Ministry of Education and Research, the Program Nucleu within MCI (PN19150201/16N/2019 and PN19060102) and project PN-III-P1-1.2-PCCDI-2017-0839/19PCCDI/2018 within PNCDI III; Slovenia—Slovenian Research Agency, Grants No. P1-0031, No. P1-0385, No. I0- 0033, No. N1-0111; Spain—Ministerio de Economía, Industria y Competitividad (FPA2017-85114-P and FPA2017-85197-P), Xunta de Galicia (ED431C 2017/ 07), Junta de Andalucía (SOMM17/6104/UGR), Feder Funds, RENATA Red Nacional Temática PHYSICAL REVIEW LETTERS 125, 121106 (2020) 121106-8 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. *[email protected] http://www.auger.org [1] J. Linsley, L. Scarsi, and B. Rossi, Phys. Rev. Lett. 6, 485 (1961). [2] J. Linsley, Phys. Rev. Lett. 10, 146 (1963). [3] Pierre Auger Collaboration, Nucl. Instrum. Methods Phys. Res., Sect. A 798, 172 (2015). [4] Telescope Array Collaboration, Nucl. Instrum. Methods Phys. Res., Sect. A 689, 87 (2012). [5] Pierre Auger Collaboration, Science 357, 1266 (2017). [6] Pierre Auger Collaboration, Astrophys. J. 868, 4 (2018). [7] Pierre Auger Collaboration, Phys. Rev. Lett. 101, 061101 (2008). [8] Pierre Auger Collaboration, companion paper, Phys. Rev. D 102, 062005 (2020). [9] Pierre Auger Collaboration, Astropart. Phys. 35, 591 (2012). [10] D. W. Newton, J. Knapp, and A. A. Watson, Astropart. Phys. 26, 414 (2007). [11] J. Hersil, I. Escobar, D. Scott, G. Clark, and S. Olbert, Phys. Rev. Lett. 6, 22 (1961). [12] Pierre Auger Collaboration, J. Cosmol. Astropart. Phys. 03 (2019) 018. [13] Pierre Auger Collaboration, Phys. Rev. D 100, 082003 (2019). [14] B. Dawson (Pierre Auger Collaboration), Proc. Sci., ICRC2019 (2019) 231, https://pos.sissa.it/358/231/pdf. [15] V. Verzi (Pierre Auger Collaboration), in Proceedings of the 33rd International Cosmic Ray Conference, Rio de Janeiro, Brazil (Centro Brasileiro de Pesquisas Físicas (CBPF), Rio de Janeiro, Brazil, 2013). [16] AIRFLY Collaboration, Astropart. Phys. 42, 90 (2013);28, 41 (2007);Nucl. Instrum. Methods Phys. Res., Sect. A 597, 50 (2008). [17] Pierre Auger Collaboration, Nucl. Instrum. Methods Phys. Res., Sect. A 613, 29 (2010). [18] Pierre Auger Collaboration, Astropart. Phys. 34, 368 (2011). [19] P. Lipari, Astropart. Phys. 97, 197 (2018) and references therein. [20] R. U. Abbasi et al. Phys. Rev. Lett. 100, 101101 (2008). [21] Telescope Array Collaboration, Astrophys. J. 768,L1 (2013). [22] C. W. Allen, Observatory 71, 157 (1951), http://adsabs .harvard.edu/full/1951Obs....71..157A. [23] Pierre Auger Collaboration, Phys. Rev. D 90, 122005 (2014). [24] Pierre Auger Collaboration, Phys. Rev. D 90, 122006 (2014). [25] A. Yushkov (Pierre Auger Collaboration), Proc. Sci., ICRC2019 (2019) 482, https://pos.sissa.it/358/482/pdf. [26] Pierre Auger Collaboration, Phys. Lett. B 762, 288 (2016). [27] Pierre Auger Collaboration, Astrophys. J. Lett. 853, L29 (2018). [28] L. Caccianiga (Pierre Auger Collaboration), Proc. Sci., ICRC2019 (2019) 206, https://pos.sissa.it/358/206/pdf. [29] K. Greisen, Phys. Rev. Lett. 16, 748 (1966); G. T. Zatsepin and V. A. Kuz’min, Pis’ma Zh. Eksp. Teor. Fiz. 4, 114 (1966) [JETP Lett. 4, 78 (1966)], http://www.jetpletters.ac .ru/ps/1624/article_24846.pdf. [30] A. M. Hillas, Phys. Lett. 24, 677 (1967). [31] V. Berezinsky, A. Gazizov, and S. Grigorieva, Phys. Lett. B 612, 147 (2005). [32] V. Berezinsky, A. Gazizov, and S. Grigorieva, Phys. Rev. D 74, 043005 (2006). [33] V. Berezinsky and S. Grigorieva, Astron. Astrophys. 199,1 (1988), http://articles.adsabs.harvard.edu/pdf/1988A%26A. ..199....1B. [34] M. Blanton, P. Blasi, and A. Olinto, Astropart. Phys. 15, 275 (2001). [35] M. S. Pshirkov, P. G. Tinyakov, and F. R. Urban, Phys. Rev. Lett. 116, 191302 (2016). [36] D. Harari, S. Mollerach, E. Roulet, and F. Sanchez, J. High Energy Phys. 03 (2002) 045. [37] G. R. Farrar, C.R. Phys. 15, 339 (2014). [38] A. Keivani, G. R. Farrar, and M. Sutherland, Astropart. Phys. 61, 47 (2015). [39] G. R. Farrar and M. Sutherland, J. Cosmol. Astropart. Phys. 05 (2019) 004. [40] D. Allard, A. Olinto, and E. Parizot, Astron. Astrophys. 473, 59 (2007). [41] R. Aloisio, V. Berezinsky, and P. Blasi, J. Cosmol. Astropart. Phys. 10 (2014) 020. [42] A. M. Taylor, M. Ahlers, and D. Hooper, Phys. Rev. D 92, 063011 (2015). [43] Pierre Auger Collaboration, J. Cosmol. Astropart. Phys. 04 (2017) 038. [44] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, Phys. Rev. C 92, 034906 (2015). [45] M. S. Muzio, M. Unger, and G. R. Farrar, Phys. Rev. D 100, 103008 (2019). [46] M. Unger, G. R. Farrar, and L. A. Anchordoqui, Phys. Rev. D92, 123001 (2015). [47] D. Hooper and A. M. Taylor, Astropart. Phys. 33, 151 (2010). [48] A. M. Taylor, M. Ahlers, and F. A. Aharonian, Phys. Rev. D 84, 105007 (2011). [49] R. Aloisio and V. Berezinsky, J. Exp. Theor. Phys. 128,52 (2019). [50] W. Hanlon (Telescope Array Collaboration), Proc. Sci., ICRC2019 (2019) 280, https://pos.sissa.it/358/280/pdf. [51] G. R. Farrar, arXiv:1902.11271. [52] W. Hanlon (Telescope Array Collaboration), EPJ Web Conf. 210, 01008 (2019). [53] N. Globus, D. Allard, and E. Parizot, Phys. Rev. D 92, 021302(R) (2015). [54] D. Biehl, D. Boncioli, A. Fedynitch, and W. Winter, Astron. Astrophys. 611, A101 (2018). PHYSICAL REVIEW LETTERS 125, 121106 (2020) 121106-9