Measurement of the Fluctuations in the Number of Muons in Extensive Air Showers with the Pierre Auger Observatory

Measurement of the Fluctuations in the Number of Muons in Extensive Air Showers with the Pierre Auger Observatory A. Aab,1P. Abreu,2M. Aglietta,3,4 J. M. Albury,5I. Allekotte,6A. Almela,7,8 J. Alvarez-Muñiz,9R. Alves Batista,1 G. A. Anastasi,10,4 L. Anchordoqui,11 B. Andrada,7S. Andringa,2C. Aramo,12 P. R. Araújo Ferreira,13 H. Asorey,7P. Assis,2 G. Avila,14 A. M. Badescu,15 A. Bakalova,16 A. Balaceanu,17 F. Barbato,18,12 R. J. Barreira Luz,2K. H. Becker,19 J. A. Bellido,5C. Berat,20 M. E. Bertaina,10,4 X. Bertou,6P. L. Biermann,21 T. Bister,13 J. Biteau,22 J. Blazek,16 C. Bleve,20 M. Boháčová,16 D. Boncioli,23,24 C. Bonifazi,25 L. Bonneau Arbeletche,26 N. Borodai,27 A. M. Botti,7J. Brack,28 T. Bretz,13 F. L. Briechle,13 P. Buchholz,29 A. Bueno,30 S. Buitink,31 M. Buscemi,32,33 K. S. Caballero-Mora,34 L. Caccianiga,35,36 A. Cancio,8,7 F. Canfora,1,37 I. Caracas,19 J. M. Carceller,30 R. Caruso,32,33 A. Castellina,3,4 F. Catalani,38 G. Cataldi,39 L. Cazon,2M. Cerda,40 J. A. Chinellato,41 K. Choi,42 J. Chudoba,16 L. Chytka,43 R. W. Clay,5A. C. Cobos Cerutti,44 R. Colalillo,18,12 A. Coleman,45 M. R. Coluccia,46,39 R. Conceição,2A. Condorelli,47,24 G. Consolati,36,48 F. Contreras,14 F. Convenga,46,39 C. E. Covault,49,†S. Dasso,50,51 K. Daumiller,52 B. R. Dawson,5J. A. Day,5R. M. de Almeida,53 J. de Jesús,7,52 S. J. de Jong,1,37 G. De Mauro,1,37 J. R. T. de Mello Neto,25,54 I. De Mitri,47,24 J. de Oliveira,53 D. de Oliveira Franco,41 V. de Souza,55 E. De Vito,46,39 J. Debatin,56 M. del Río,14 O. Deligny,57 H. Dembinski,52 N. Dhital,27 A. Di Matteo,4C. Dobrigkeit,41 J. C. D’Olivo,58 R. C. dos Anjos,59 M. T. Dova,60 J. Ebr,16 R. Engel,56,52 I. Epicoco,46,39 M. Erdmann,13 C. O. Escobar,61 A. Etchegoyen,7,8 H. Falcke,1,62,37 J. Farmer,63 G. Farrar,64 A. C. Fauth,41 N. Fazzini,65 F. Feldbusch,66 F. Fenu,10,4 B. Fick,67 J. M. Figueira,7A. Filipčič,68,69 T. Fodran,1M. M. Freire,70 T. Fujii,63,‡ A. Fuster,7,8 C. Galea,1C. Galelli,35,36 B. García,44 A. L. Garcia Vegas,13 H. Gemmeke,66 F. Gesualdi,7,52 A. Gherghel-Lascu,17 P. L. Ghia,57 U. Giaccari,1M. Giammarchi,36 M. Giller,71 J. Glombitza,13 F. Gobbi,40 F. Gollan,7 G. Golup,6M. Gómez Berisso,6P. F. Gómez Vitale,14 J. P. Gongora,14 N. González,7I. Goos,6,52 D. Góra,27 A. Gorgi,3,4 M. Gottowik,19 T. D. Grubb,5F. Guarino,18,12 G. P. Guedes,72 E. Guido,4,10 S. Hahn,52,7 R. Halliday,49 M. R. Hampel,7 P. Hansen,60 D. Harari,6V. M. Harvey,5A. Haungs,52 T. Hebbeker,13 D. Heck,52 G. C. Hill,5C. Hojvat,65 J. R. Hörandel,1,37 P. Horvath,43 M. Hrabovský,43 T. Huege,52,31 J. Hulsman,7,52 A. Insolia,32,33 P. G. Isar,73 J. A. Johnsen,74 J. Jurysek,16 A. Kääpä,19 K. H. Kampert,19 B. Keilhauer,52 J. Kemp,13 H. O. Klages,52 M. Kleifges,66 J. Kleinfeller,40 M. Köpke,56 G. Kukec Mezek,69 B. L. Lago,75 D. LaHurd,49 R. G. Lang,55 N. Langner,13 M. A. Leigui de Oliveira,76 V. Lenok,52 A. Letessier-Selvon,77 I. Lhenry-Yvon,57 D. Lo Presti,32,33 L. Lopes,2R. López,78 R. Lorek,49 Q. Luce,56 A. Lucero,7 J. P. Lundquist,69 A. Machado Payeras,41 G. Mancarella,46,39 D. Mandat,16 B. C. Manning,5J. Manshanden,79 P. Mantsch,65 S. Marafico,57 A. G. Mariazzi,60 I. C. Mariş,42 G. Marsella,46,39 D. Martello,46,39 H. Martinez,55 O. Martínez Bravo,78 M. Mastrodicasa,23,24 H. J. Mathes,52 J. Matthews,80 G. Matthiae,81,82 E. Mayotte,19 P. O. Mazur,65 G. Medina-Tanco,58 D. Melo,7A. Menshikov,66 K.-D. Merenda,74 S. Michal,43 M. I. Micheletti,70 L. Miramonti,35,36 S. Mollerach,6 F. Montanet,20 C. Morello,3,4 M. Mostafá,83 A. L. Müller,7,52 M. A. Muller,41,25,§ K. Mulrey,31 R. Mussa,4M. Muzio,64 W. M. Namasaka,19 L. Nellen,58 M. Niculescu-Oglinzanu,17 M. Niechciol,29 D. Nitz,67,∥D. Nosek,84 V. Novotny,84 L. Nožka,43 A. Nucita,46,39 L. A. Núñez,85 M. Palatka,16 J. Pallotta,86 P. Papenbreer,19 G. Parente,9A. Parra,78 M. Pech,16 F. Pedreira,9J. P¸ekala,27 R. Pelayo,87 J. Peña-Rodriguez,85 J. Perez Armand,26 M. Perlin,7,52 L. Perrone,46,39 S. Petrera,47,24 T. Pierog,52 M. Pimenta,2V. Pirronello,32,33 M. Platino,7B. Pont,1M. Pothast,37,1 P. Privitera,63 M. Prouza,16 A. Puyleart,67 S. Querchfeld,19 J. Rautenberg,19 D. Ravignani,7M. Reininghaus,52,7 J. Ridky,16 F. Riehn ,2M. Risse,29 P. Ristori,86 V. Rizi,23,24 W. Rodrigues de Carvalho,26 J. Rodriguez Rojo,14 M. J. Roncoroni,7M. Roth,52 E. Roulet,6A. C. Rovero,50 P. Ruehl,29 S. J. Saffi,5A. Saftoiu,17 F. Salamida,23,24 H. Salazar,78 G. Salina,82 J. D. Sanabria Gomez,85 F. Sánchez,7 E. M. Santos,26 E. Santos,16 F. Sarazin,74 R. Sarmento,2C. Sarmiento-Cano,7R. Sato,14 P. Savina,46,39,57 C. M. Schäfer,52 V. Scherini,39 H. Schieler,52 M. Schimassek,56,7 M. Schimp,19 F. Schlüter,52,7 D. Schmidt,56 O. Scholten,88,31 P. Schovánek,16 F. G. Schröder,45,52 S. Schröder,19 J. Schulte,13 S. J. Sciutto,60 M. Scornavacche,7,52 R. C. Shellard,89 G. Sigl,79 G. Silli,7,52 O. Sima,17,¶ R. Šmída,63 P. Sommers,83 J. F. Soriano,11 J. Souchard,20 R. Squartini,40 M. Stadelmaier,52,7 D. Stanca,17 S. Stanič,69 J. Stasielak,27 P. Stassi,20 A. Streich,56,7 M. Suárez-Durán,85 T. Sudholz,5 T. Suomijärvi,22 A. D. Supanitsky,7J. Šupík,43 Z. Szadkowski,90 A. Taboada,56 A. Tapia,91 C. Timmermans,37,1 O. Tkachenko,52 P. Tobiska,16 C. J. Todero Peixoto,38 B. Tom´e,2G. Torralba Elipe,9A. Travaini,40 P. Travnicek,16 C. Trimarelli,23,24 M. Trini,69 M. Tueros,60 R. Ulrich,52 M. Unger,52 L. Vaclavek,43 M. Vacula,43 J. F. Vald´es Galicia,58 I. Valiño,47,24 L. Valore,18,12 E. Varela,78 V. Varma K. C.,7,52 A. Vásquez-Ramírez,85 D. Veberič,52 C. Ventura,54 I. D. Vergara Quispe,60 V. Verzi,82 J. Vicha,16 J. Vink,92 S. Vorobiov,69 H. Wahlberg,60 A. A. Watson,93 M. Weber,66 PHYSICAL REVIEW LETTERS 126, 152002 (2021) 0031-9007=21=126(15)=152002(11) 152002-1 Published by the American Physical Society A. Weindl,52 L. Wiencke,74 H. Wilczyński,27 T. Winchen,31 M. Wirtz,13 D. Wittkowski,19 B. Wundheiler,7A. Yushkov,16 O. Zapparrata,42 E. Zas,9D. Zavrtanik,69,68 M. Zavrtanik,68,69 L. Zehrer,69 and A. Zepeda94 (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, South Australia, 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 9Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de Santiago de Compostela, Santiago de Compostela, Spain 10Dipartimento di Fisica, Universit`a Torino, Torino, Italy 11Department of Physics and Astronomy, Lehman College, City University of New York, Bronx, New York, USA 12INFN, Sezione di Napoli, Napoli, Italy 13III. Physikalisches Institut A, RWTH Aachen University, Aachen, Germany 14Observatorio Pierre Auger and Comisión Nacional de Energía Atómica, Malargüe, Argentina 15University Politehnica of Bucharest, Bucharest, Romania 16Institute of Physics of the Czech Academy of Sciences, Prague, Czech Republic 17“Horia Hulubei”National Institute for Physics and Nuclear Engineering, Bucharest-Magurele, Romania 18Dipartimento di Fisica “Ettore Pancini,”Universit`a di Napoli “Federico II,”Napoli, Italy 19Department of Physics, Bergische Universität Wuppertal, Wuppertal, Germany 20Universit´e Grenoble Alpes, CNRS, Grenoble Institute of Engineering Universit´e Grenoble Alpes, LPSC-IN2P3, 38000 Grenoble, France 21Max-Planck-Institut für Radioastronomie, Bonn, Germany 22Universit´e Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France 23Dipartimento di Scienze Fisiche e Chimiche, Universit`a dell’Aquila, L’Aquila, Italy 24INFN Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy 25Instituto de Física, Universidade Federal do Rio de Janeiro, Rio de Janeiro, Brazil 26Instituto de Física, Universidade de São Paulo, São Paulo, Brazil 27Institute of Nuclear Physics PAN, Krakow, Poland 28Colorado State University, Fort Collins, Colorado, USA 29Universität Siegen, Fachbereich 7 Physik—Experimentelle Teilchenphysik, Siegen, Germany 30Universidad de Granada and C.A.F.P.E., Granada, Spain 31Vrije Universiteit Brussels, Brussels, Belgium 32Dipartimento di Fisica e Astronomia, Universit`a di Catania, Catania, Italy 33INFN, Sezione di Catania, Catania, Italy 34Universidad Autónoma de Chiapas, Tuxtla Guti´errez, Chiapas, M´exico 35Dipartimento di Fisica, Universit`a di Milano, Milano, Italy 36INFN, Sezione di Milano, Milano, Italy 37Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF), Science Park, Amsterdam, The Netherlands 38Universidade de São Paulo, Escola de Engenharia de Lorena, Lorena, São Paulo, Brazil 39INFN, Sezione di Lecce, Lecce, Italy 40Observatorio Pierre Auger, Malargüe, Argentina 41Universidade Estadual de Campinas, IFGW, Campinas, São Paulo, Brazil 42Universit´e Libre de Bruxelles (ULB), Brussels, Belgium 43Palacky University, RCPTM, Olomouc, Czech Republic 44Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), and Universidad Tecnológica Nacional—Facultad Regional Mendoza (CONICET/CNEA), Mendoza, Argentina 45Department of Physics and Astronomy, Bartol Research Institute, University of Delaware, Newark, Delaware, USA 46Dipartimento di Matematica e Fisica “E. De Giorgi,”Universit`a del Salento, Lecce, Italy 47Gran Sasso Science Institute, L’Aquila, Italy 48Dipartimento di Scienze e Tecnologie Aerospaziali, Politecnico di Milano, Milano, Italy 49Case Western Reserve University, Cleveland, Ohio, USA 50Instituto de Astronomía y Física del Espacio (IAFE, CONICET-UBA), Buenos Aires, Argentina PHYSICAL REVIEW LETTERS 126, 152002 (2021) 152002-2 51Departamento 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 52Karlsruhe Institute of Technology, Institut für Kernphysik, Karlsruhe, Germany 53Universidade Federal Fluminense, EEIMVR, Volta Redonda, Rio de Janeiro, Brazil 54Universidade Federal do Rio de Janeiro (UFRJ), Observatório do Valongo, Rio de Janeiro, Brazil 55Instituto de Física de São Carlos, Universidade de São Paulo, São Carlos, São Paulo, Brazil 56Karlsruhe Institute of Technology, Institute for Experimental Particle Physics (ETP), Karlsruhe, Germany 57CNRS/IN2P3, IJCLab, Universit´e Paris-Saclay, Orsay, France 58Universidad Nacional Autónoma de M´exico, M´exico, Distrito Federal, M´exico 59Universidade Federal do Paraná, Setor Palotina, Palotina, Brazil 60IFLP, Universidad Nacional de La Plata and CONICET, La Plata, Argentina 61Fermi National Accelerator Laboratory, Batavia, Illinois, USA 62Stichting Astronomisch Onderzoek in Nederland (ASTRON), Dwingeloo, The Netherlands 63University of Chicago, Enrico Fermi Institute, Chicago, Illinois, USA 64New York University, New York, New York, USA 65Fermi National Accelerator Laboratory, Fermilab, Batavia, Illinois, USA 66Karlsruhe Institute of Technology, Institut für Prozessdatenverarbeitung und Elektronik, Karlsruhe, Germany 67Michigan Technological University, Houghton, Michigan, 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 Universidad Nacional de Rosario, Rosario, Argentina 71Faculty of Astrophysics, University of Łód´z, Łód´z, Poland 72Universidade Estadual de Feira de Santana, Feira de Santana, Brazil 73Institute of Space Science, Bucharest-Magurele, Romania 74Colorado School of Mines, Golden, Colorado, USA 75Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Nova Friburgo, Brazil 76Universidade Federal do ABC, Santo Andr´e, São Paulo, Brazil 77Laboratoire de Physique Nucl´eaire et de Hautes Energies (LPNHE), Sorbonne Universit´e, Universit´e de Paris, CNRS-IN2P3, Paris, France 78Benem´erita Universidad Autónoma de Puebla, Puebla, M´exico 79Universität Hamburg, II. Institut für Theoretische Physik, Hamburg, Germany 80Louisiana State University, Baton Rouge, Louisiana, USA 81Dipartimento di Fisica, Universit`a di Roma “Tor Vergata,”Roma, Italy 82INFN, Sezione di Roma “Tor Vergata,”Roma, Italy 83Pennsylvania State University, University Park, Pennsylvania, USA 84Charles University, Faculty of Mathematics and Physics, Institute of Particle and Nuclear Physics, Prague, Czech Republic 85Universidad Industrial de Santander, Bucaramanga, Colombia 86Centro de Investigaciones en Láseres y Aplicaciones, CITEDEF and CONICET, Villa Martelli, Argentina 87Unidad Profesional Interdisciplinaria en Ingeniería y Tecnologías Avanzadas del Instituto Polit´ecnico Nacional (UPIITA-IPN), M´exico, Distrito Federal, M´exico 88KVI—Center for Advanced Radiation Technology, University of Groningen, Groningen, The Netherlands 89Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, Brazil 90Faculty of High-Energy Astrophysics, University of Łód´z, Łód´z, Poland 91Universidad de Medellín, Medellín, Colombia 92Faculty of Science, Universiteit van Amsterdam, Amsterdam, The Netherlands 93School of Physics and Astronomy, University of Leeds, Leeds, United Kingdom 94Centro de Investigación y de Estudios Avanzados del IPN (CINVESTAV), M´exico, Distrito Federal, M´exico (Received 4 October 2020; revised 28 January 2021; accepted 12 February 2021; published 16 April 2021) We present the first measurement of the fluctuations in the number of muons in extensive air showers produced by ultrahigh energy cosmic rays. We find that the measured fluctuations are in good agreement with predictions from air shower simulations. This observation provides new insights into the origin of the previously reported deficit of muons in air shower simulations and constrains models of hadronic 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. Funded by SCOAP3. PHYSICAL REVIEW LETTERS 126, 152002 (2021) 152002-3 interactions at ultrahigh energies. Our measurement is compatible with the muon deficit originating from small deviations in the predictions from hadronic interaction models of particle production that accumulate as the showers develop. DOI: 10.1103/PhysRevLett.126.152002 Introduction.—Ultrahigh energy cosmic rays (UHECRs) are particles coming from outer space, with energies exceeding 1018 eV. They provide the only experimental opportunity to explore particle physics beyond energies reachable by Earth-based accelerators, which go up to cosmic ray energies of 9×1016 eV. The Pierre Auger Observatory [1] detects extensive air showers that are initiated by the UHECRs colliding with the nuclei in the atmosphere. Information about UHECRs is extracted using simulations based on hadronic interaction models which rely on extrapolations of accelerator measurementsto unexplored regions of phase space, most notably the forward and highest-energy region. In addition, accelerator experiments at the highest energies either probe the interactions between protons or of protons with heavy nuclei, while most interactions within air showers are between pions and light nuclei. A further challenge is that the UHECR mass has to be measured despite not being yet completely decoupled from the hadronic uncertainties. The observable with the least dependence on hadronic interactions is the atmospheric depth at which the longitudinal development of the electromagnetic (EM) component of the shower reaches the maximum number of particles, namely, Xmax [2]. In hadronic cascades, the energy of each interacting particle is distributed among the secondaries, mostly pions. Neutral pions rapidly decay into two photons, feeding a practically decoupled electromagnetic cascade (other resonances decaying into π0’s, electrons, and/or photons also contribute). Charged pions (and other long-lived mesons like kaons) tend to further interact until their individual energies are below a critical value, below which they are more likely to decay. Muons, which are products of hadronic decays, are thus predominantly produced in the final shower stages. In sufficiently inclined showers, the pure EM component is absorbed in the atmosphere and the particles that reach the ground (muons and muon decay products) directly sample the muon content [3,4], reflecting the hadronic component of the shower. Air showers are mainly detected at the Pierre Auger Observatory by the surface detector (SD), an array of water-Cherenkov detector stations, and the fluorescence detector (FD), consisting of 24 fluorescence telescopes. By selecting the subsample of events reconstructed with both the SD and FD, and with zenith angles exceeding 62°, both the muon content and the energy of the shower are simultaneously measured. The results obtained indicate that all the simulations underestimate the number of muons in the showers [5,6]. These analyses come with the caveat that they cannot distinguish a muon rescaling from a shift in the absolute energy scale of the FD measurement. However, muon content and energy scale were disentangled in a complementary technique based on showers with zenith angles below 60°. Using the longitudinal profile of the shower in the atmosphere obtained with the FD and the signals at the ground measured with the SD, it was shown that the muonic component still has to be scaled up to match observed data, while no rescaling of the EM component and the FD energy is required [7]. The measurements with the FD also show that both the position of the shower maximum in the atmosphere (Xmax) and the entire shape of the EM shower are well described by the simulations [8,9].Atlower energies, down to ∼1017.3eV, in a measurement using the subarray of buried scintillators of the Pierre Auger Observatory, a direct count of the muons independent of EM contamination was obtained, which also shows that simulations produce too few muons [10]. There is much evidence that all the simulations underpredict the average number of muons in the showers: a comprehensive study of muon number measurements made with different experiments has shown that the muon deficit in simulations starts around ∼1016 eV and steadily increases with energy. Depending on model and experiment, the deficit at ∼1020 eV ranges between tens of percent up to a factor of 2 [11]. The increased statistics obtained at the Pierre Auger Observatory allow us to now take a further step and explore fluctuations in the number of muons between showers, hereinafter referred to as “physical fluctuations." The ratio of the physical fluctuations to the average number of muons (relative fluctuations) has been shown to be mostly dominated by the first interaction, rather than the lower energy interactions deeper in the shower development [12,13]. Here, we exploit the sensitivity of fluctuations to the first interaction to explore hadronic interactions well above the energies achievable in accelerator experiments. Methodology.—Our analysis here is based on the set of inclined air showers (62°<θ<80°) that are reconstructed both with the SD and FD between January 1, 2004 and December 31, 2017. For each event, we obtain independent measurements of the muon content (with the SD) and the calorimetric energy (with the FD). To ensure the showers can be reconstructed with small uncertainties, we select only events with at least four triggered stations in the SD array and we further require that all the stations surrounding the impact point of the shower on the ground are operational at the time of the event. Only events with good atmospheric conditions (few clouds and a low aerosol content) are accepted in order PHYSICAL REVIEW LETTERS 126, 152002 (2021) 152002-4 to guarantee a good energy reconstruction with the FD. In addition, it is required that the entire shower profile and, in particular, Xmax is within the field of view of our telescopes. Since heavy primaries penetrate the atmosphere less than light ones, the acceptance with this selection would be mass dependent. To avoid this bias, we constrain the field of view to the region where all values of Xmax are accepted. Further details are given in [5,14]. These selection criteria result in a total number of events of 786. In addition, only events with energy larger than 4×1018 eV, which ensures full trigger efficiency of the SD [3], are used to extract the fluctuations (281 events). The number of muons is reconstructed by fitting a 2D model of the lateral profile of the muon density at the ground to the observed signals in the SD array. The free parameters of the fit are the zenith and azimuth angles of the shower, the impact point of the shower on the ground (shower core position), and a normalization factor with respect to a reference muon density profile in simulated proton showers at 1019 eV [3]. There exists a residual pure EM component in showers with low zenith angles and stations very close to the shower core position (at 400m and 64° it is ∼6%), which has been subtracted using a parametrization [4]. The dimensionless normalization factor we obtain from the fit is then transformed to the dimensionless quantity Rμ, which is given by the integrated number of muons at the ground divided by a reference given by the average number of muons in simulated proton showers at 1019 eV and the given zenith angle. At 1019 eV and an inclination of 60°, Rμ¼1corresponds to 2.148 ×107 muons. For more details, see [5]. In the following, we refer to Rμas the number of muons for short. The calorimetric energy of the air showers Ecal is reconstructed by integrating the longitudinal shower profiles observed with the FD [9,15]. The total energy of the shower is then obtained by adding the average energy carried away by muons and neutrinos, the so-called invisible energy E¼Ecal þEinv.At1019 eV, Einv accounts for 14% of the total energy in air showers [16–20]. In Fig. 1the muon number Rμis shown as a function of the measured energy. Markers on the top of the frame define the bins in energy for which we will extract the fluctuations, with the number of events in each bin shown above. The bins are chosen such that the number of events in each is similar. Based on models of air shower development and given the gradual change of the composition in this energy range (single logarithmic dependence on energy) [8,21–23], the number of muons is related to the primary energy by a single power law hRμiðEÞ¼a½E=ð1019 eVÞb;ð1Þ which can be fitted following a procedure described in the text below. The best-fit parameters are given at the beginning of the next section. The scattering in the data has three sources: experimental uncertainties in the energy sEand in the muon number sμfrom event reconstruction (both represented by the error bars), and the physical fluctuations in the muon number denoted as σ. Given Eq. (1), the variance of the muon number is σ2þb2hsEi2þhsμi2. In this Letter, we adopt a method based on maximizing the likelihood of a probability distribution function (PDF). The PDF incorporates the various contributions to the fluctuations, treating each energy bin independently while also accounting for the effect of the migration of events between bins [5,24]. The model assumes that measurements of Eand Rμfollow Gaussian distributions centered at the true value, with widths given by the detector resolution sEand sμ, which are the uncertainties obtained in each individual event reconstruction [3,25]. Physical fluctuations are also assumed to follow a Gaussian distribution of width σ. Simulations have shown this is an acceptable approximation given the event number in each bin. The total PDF is obtained through the convolution of the detector response and the physical fluctuations with the probability distribution of the hybrid events measured at the Pierre Auger Observatory. The log-likelihood function is then given by ln Lða; b; ˆσ1;…;ˆσ6Þ¼X i ln X 6 k¼0ZEk Ek−1 dEhðEÞCðEÞ × exp − 1 2ðEi−EÞ2 s2 E × exp − 1 2ðRμ;i −hRμiðEÞÞ2 s2 μþðˆσk·hRμiðEÞÞ2: ð2Þ FIG. 1. Number of muons as a function of the measured energy. The black line is the fitted hRμi¼a½E=ð1019 eVÞb. Markers on the top of the frame define the bins in which the fluctuations are evaluated. The numbers givethe events in each bin.The effectof the uncertainty of the absolute energy scale is indicated by σsysðEÞ. The best-fit values for parameters a; b and the deviance per degree of freedom (n.d.f.) in the fit are shown on the lower right. PHYSICAL REVIEW LETTERS 126, 152002 (2021) 152002-5 The probability of hybrid events hðEÞ(product of the energy spectrum of cosmic rays and the efficiency of detection) can be obtained from the data, as explained in and [10,24,26]. The rhs of Eq. (2) depends on the parameters aand bvia Eq. (1). To obtain the energy dependence of the fluctuations, we parametrize σby six independent values such that σðEÞ¼ˆσk·hRμiðEÞ, where the constants ˆσkare the relative fluctuations in the kth energy bin with limits ½Ek−1;E k, where kruns from one to six. In Eq. (2),k¼0corresponds to the contributions from the interval ½0;E thr, where the SD is not fully efficient. The fluctuations here are assumed to take the value of the first fitted bin ˆσ0≡ˆσ1. The sum over the index iin Eq. (2) (the usual sum over the log-likelihoods of events) includes only events above the energy threshold of 4×1018 eV. The function CðEÞis the normalization factor from the double Gaussian. The result of the fit for the parameters aand bare shown in Fig. 1. The fluctuations are shown in Fig. 2. The distribution of the number of muons and the PDF in the individual energy bins can be found in the Supplemental Material [17]. The dominant systematic uncertainties of σcome from the uncertainties in the resolutions sEand sμ.Forsμwe estimate the uncertainty using simulations and data. In simulations, the uncertainty was estimated by the spread in a sample of simulated showers, where each shower is reconstructed multiple times, each time changing only the impact point at the ground. For data, we reconstruct the same event multiple times, leaving out the signals from one of the detector stations. The average relative resolution hsμ=Rμiand its systematic uncertainty is thus ð10 3Þ% at 1019 eV. We verified the values of sEby studying the difference in the energy reconstruction of events measured independently by two or more FD stations. The width of the distribution of these energy differences is found to be compatible with sE. We therefore take the statistical 1-σuncertainties of this cross check as a conservative upper limit of the systematic uncertainty of sE[27]. The average relative energy resolution hsE=Eiis about ð8.42.9Þ%at 1019 eV. We have further confirmed that there are no significant contributions to the fluctuations from differences between the individual FD stations, neither related to the longtime performance evolution of the SD and FD detectors. Any residual electromagnetic component in the signal would affect the lower zenith angles more. We therefore split the event sample at the median zenith angle (66°) and compare the resulting fluctuations. We find no significant difference between the more and the less inclined sample. In another test, we do find a small modulation of hRμi with the azimuth angle (<1%), which we correct for. This modulation is related to the approximations used in the reconstruction, which deal with the azimuthal asymmetry of the muon densities at the ground due to the Earth’s magnetic field [3]. Finally, we have run an end-to-end validation of the whole analysis method described in this Letter on samples of simulated proton, helium, oxygen, and iron showers. Because of the almost linear relation between Rμand E, the systematic uncertainty on σdue to the uncertainty of the absolute energy scale of 14% [25] practically cancels out in the relative fluctuations. The systematic uncertainty in the absolute scale of Rμof 11% [5] drops out for the same reason. The systematic effects for the bin around 1019 eV are summarized in Table I. Over all energies, the systematic uncertainties are below 8%. Results and discussion.—The best-fit value for the average relative number of muons at 1019 eV (parameter a) is hRμið1019eVÞ¼1.860.02ðstatÞþ0.36 −0.31 ðsystÞ. For the slope (parameter b) we find dhln Rμi=d ln E¼0.99  0.02ðstatÞþ0.03 −0.03 ðsystÞ. These values are consistent with the values previously reported [5,17]. FIG. 2. Measured relative fluctuations in the number of muons as a function of the energy and the predictions from three interaction models for proton (red) and iron (blue) showers. The gray band represents the expectations from the measured mass composition interpreted with the interaction models. The statistical uncertainty in the measurement is represented by the error bars. The total systematic uncertainty is indicated by the square brackets. TABLE I. Contributions to the systematic uncertainty in the relative fluctuations around 1019 eV (1018.97–1019.15 eV). The central value is σ=hRμi¼0.102 0.029ðstatÞ0.007ðsystÞ. Source of uncertainty Uncertainty (%) Eabsolute scale hEi<0.1 Eresolution sE4.6 Rμabsolute scale hRμi0.5 Rμresolution sμ5.2 Rμazimuthal modulation hRμiðϕÞ0.5 Total systematics 7.0 PHYSICAL REVIEW LETTERS 126, 152002 (2021) 152002-6 The measured relative fluctuations as a function of the energy are shown in Fig. 2. We note that the measurement falls within the range that is expected from current hadronic interaction models for pure proton and pure iron primaries [28–36]. To estimate the effect of a mixed composition, we take the fractions of the four mass components (proton, helium, nitrogen, and iron) derived from the Xmax measurements [8,37,38] and, using the simulations of the pure primaries, calculate the corresponding fluctuations in the number of muons. The gray band in Fig. 2encompasses the predicted σ=hRμiof the three interaction models QGSJET II-04, EPOS-LHC, and Sibyll 2.3d given the inferred composition mix for each [17]. In Fig. 3, the effects of different composition scenarios on both the fluctuations and the average number of muons can be shown by drawing, at a fixed primary energy of 1019 eV, the relative fluctuations σ=hRμiagainst the average number of muons hRμi. Given any one of the interaction models, any particular mixture of the four components p, He, N, and Fe falls somewhere within one of the areas enclosed by the corresponding colored lines. The points of pure composition in this contour are labeled accordingly. For each model, the expected values for σ=hRμiand hRμigiven the composition mixture obtained from the Xmax measurements [8] is indicated within each contour by the correspondingly colored star marker. The shaded areas surrounding the star markers indicate the statistical and systematic uncertainties inherited from the Xmax measurements [39]. Finally, our measurement with statistical and systematic uncertainty is shown by the black marker. Within the uncertainty, none of the predictions from the interaction models and the Xmax composition (star markers) are consistent with our measurement. The predictions from the interaction models QGSJET II-04, EPOS-LHC, and Sibyll 2.3d can be reconciled with our measurement by an increase in the average number of muons of 43%, 35%, and 26%, respectively. For the fluctuations, no rescaling is necessary for any model. Taken together, the average value and fluctuations of the muon flux constrain the way hadronic interaction models should be changed to agree with air shower data. To see this, we briefly discuss the origin of the fluctuations. The average number of muons in a proton shower of energy Ehas been shown in simulations to scale as hN μi¼CEβ, where β≃0.9[12,13,22,23]. If we assume all the secondaries from the first interaction produce muons following the same relation as given for protons above, we obtain the number of muons in the shower as Nμ¼X m j¼1 CEβ j¼hN μiX m j¼1 xβ j¼hN μiα1;ð3Þ where index jruns over msecondary particles which reinteract hadronically and xj¼Ej=E is the fraction of energy fed to the hadronic shower by each [41]. In this expression, the fluctuations in Nμare induced by α1in the first generation, which fluctuates because the multiplicity m and the energies xjof the secondaries fluctuate [13]. We can continue this reasoning for the subsequent generations to obtain Nμ hN μi¼α1α2αiαn;ð4Þ here the subindex iruns over ngenerations, until the cascade stops. We note that, for the calculation of α2, in the second generation, there are mparticles contributing. Assuming the distributions of the α’s for each one are similar, when adding up the muons produced by each, the fluctuations produced by one are statistically likely to be compensated by another. In other words, the α2distribution is narrower by a factor ∼1=ffiffiffiffi m p. The deeper the generation, the sharper the corresponding αiis expected to be. As a result, the dominant part of the fluctuations comes from the first interaction. This has also been observed with simulations. The model can be generalized for primary nuclei with mass Ausing the superposition model and fixing the number of participants to Aprotons, which reduces the different contributions to the fluctuations by a factor ∼1=ffiffiffiffi A p. FIG. 3. Data (black, with error bars) compared to models for the fluctuations and the average number of muons for showers with a primary energy of 1019 eV. Fluctuations are evaluated in the energy range from 1018.97 to 1019.15 eV. The statistical uncertainty is represented by the error bars. The total systematic uncertainty is indicated by the square brackets. The expectation from the interaction models for any mixture of the four components p, He, N, Fe is illustrated by the colored contours. The values preferred by the mixture derived from the Xmax measurements are indicated by the star symbols. The shaded areas show the regions allowed by the statistical and systematic uncertainties of the Xmax measurement [39]. PHYSICAL REVIEW LETTERS 126, 152002 (2021) 152002-7 There are two options to increase the average number of muons in air showers. One is to increase αin a specific generation, notably the first where the energy is the highest and exotic phenomena could conceivably play a role, i.e., α1→α1þδα1. Note that, if only the first generation is modified (implying some sort of threshold effect for new physics), the increase in Nμis linear with the modification. There are several examples in the literature where this approach has been used assuming different mechanisms [43–47]. For the fluctuations, the change depends on the model. Alternatively, the number of muons can be increased by introducing small deviations in the hadronic energy fraction δα in all generations. Accumulated along a number nof generations, these small deviations build up as Nμ∝ðαþδαÞn. For instance, a 5% deviation per generation converts into ∼30% deviation after six generations [48]. On the other hand, a change of 5% in the fluctuations of αis not amplified in the muon fluctuations because of the suppression in later generations. This approach characterizes the increase in the number of muons in the current hadronic interaction models with regard to previous models [32,50–54]. It is also compatible with the increase of the discrepancy in the average number of muons across a wide range of energies reported in [11]. The present analysis finding that fluctuations are consistent with model predictions means that the increase in muon number may be a small effect accumulating over many generations or a very particular modification of the first interaction that changes Nμwithout changing the fluctuations [17]. Summary.—We have presented for the first time a measurement of the fluctuations in the number of muons in inclined air showers, as a function of the UHECR primary energy. Within the current uncertainties, the relative fluctuations show no discrepancy with respect to the expectation from current high-energy hadronic interaction models and the composition taken from Xmax measurements. This agreement between models and data for the fluctuations, combined with the significant deficit in the predicted total number of muons, points to the origin of the models’muon deficit being a small deficit at every stage of the shower that accumulates along the shower development, rather than a discrepancy in the first interaction. Adjustments to models to address the current muon deficit must therefore not alter the predicted relative fluctuations. The Pierre Auger Observatory is currently undergoing an upgrade that includes the deployment of scintillators on top of the SD stations [55] to help disentangle the muonic and electromagnetic content of the showers, as well as an array of radio antennas [56]. It has been shown that radio arrays can provide an estimate of the calorimetric energy [57], and therefore, it will soon be possible to perform an analysis similar to the one presented here with much larger statistics using hybrid events measured by the high-duty-cycle radio and surface detector arrays [56]. 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´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); Ministry of Education, Youth and Sports of the Czech Republic— Grants No. MSMT CR LTT18004, No. LM2015038, No. LM2018102, No. CZ.02.1.01/0.0/0.0/16_013/ 0001402, No. CZ.02.1.01/0.0/0.0/18_046/0016010, and No. 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. ANR11-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) Grant No. 167733, Universidad Nacional Autónoma de M´exico (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 Science and Higher Education, Grant No. DIR/WK/2018/11, National Science Centre, Grants No. 2013/08/M/ST9/00322, PHYSICAL REVIEW LETTERS 126, 152002 (2021) 152002-8 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. I00033, 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 de Astropartículas (FPA2015-68783-REDT), and María de Maeztu Unit of Excellence (MDM-2016-0692); U.S.— Department of Energy, Awards No. DE-AC0207CH11359, No. DE-FR02-04ER41300, No. DE-FG0299ER41107, 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 †Also at Radboud Universtiy Nijmegen, Nijmegen, The Netherlands. ‡Present address: Hakubi Center for Advanced Research and Graduate School of Science, Kyoto University, Kyoto, Japan. §Also at Universidade Federal de Alfenas, Poços de Caldas, Brazil. ∥Also at Karlsruhe Institute of Technology, Karlsruhe, Germany. ¶Also at University of Bucharest, Physics Department, Bucharest, Romania. [1] A. Aab et al. (Pierre Auger Collaboration), The Pierre Auger Cosmic Ray Observatory, Nucl. Instrum. Methods Phys. Res., Sect. A 798, 172 (2015). [2] J. Linsley, Structure of large showers at depth 834 g=cm2, applications, in Proceedings of the 15th Int. Cosmic Ray Conf., Plovdiv, Bulgaria (Bulgarian Academy of Sciences, Bulgaria, 1977), Vol. 12, p. 89. [3] A. Aab et al. (Pierre Auger Collaboration), Reconstruction of inclined air showers detected with the Pierre Auger Observatory, J. Cosmol. Astropart. Phys. 08 (2014) 019. [4] I. Valiño, J. Alvarez-Muñiz, M. Roth, R. A. Vazquez, and E. Zas, Characterisation of the electromagnetic component in ultra-high energy inclined air showers, Astropart. Phys. 32, 304 (2010). [5] A. Aab et al. (Pierre Auger Collaboration), Muons in air showers at the Pierre Auger Observatory: Mean number in highly inclined events, Phys. Rev. D 91, 032003 (2015); Erratum, Phys. Rev. D91, 059901 (2015). [6] S. J. Sciutto, Air showers, hadronic models, and muon production, EPJ Web Conf. 210, 02007 (2019). [7] A. Aab et al. (Pierre Auger Collaboration), Testing Hadronic Interactions at Ultrahigh Energies with Air Showers Measured by the Pierre Auger Observatory, Phys. Rev. Lett. 117, 192001 (2016). [8] A. Aab et al. (Pierre Auger Collaboration), Depth of maximum of air-shower profiles at the Pierre Auger Observatory. I. Measurements at energies above 1017.8eV, Phys. Rev. D 90, 122005 (2014). [9] A. Aab et al. (Pierre Auger Collaboration), Measurement of the average shape of longitudinal profiles of cosmic-ray air showers at the Pierre Auger Observatory, J. Cosmol. Astropart. Phys. 03 (2019) 018. [10] A. Aab et al. (Pierre Auger Collaboration), Direct measurement of the muonic content of extensive air showers between 2×1017 and 2×1018 eV at the Pierre Auger Observatory, Eur. Phys. J. C 80, 751 (2020). [11] L. Cazon et al. (EAS-MSU1, IceCube, KASCADE-Grande, NEVOD-DECOR, Pierre Auger, SUGAR, Telescope Array, and Yakutsk EASArray Collaborations), Working group report on the combined analysis of muon density measurements from eight air shower experiments, Proc. Sci., ICRC2019 (2019) 214. [12] S. Fukui, H. Hasegawa, T. Matano, I. Miura, M. Oda, K. Suga, G. Tanahashi, and Y. Tanaka, A study on the structure of the extensive air shower, Prog. Theor. Phys. Suppl. 16,1 (1960). [13] L. Cazon, R. Conceição, and F. Riehn, Probing the energy spectrum of hadrons in proton air interactions at ultrahigh energies through the fluctuations of the muon content of extensive air showers, Phys. Lett. B 784, 68 (2018). [14] P. Abreu, M. Aglietta, E. J. Ahn, D. Allard, I. Allekotte, J. Allen, J. Alvarez Castillo, J. Alvarez-Muñiz, M. Ambrosio, A. Aminaei et al., The exposure of the hybrid detector of the Pierre Auger Observatory, Astropart. Phys. 34, 368 (2011). [15] J. Abraham et al. (Pierre Auger Collaboration), The fluorescence detector of the Pierre Auger Observatory, Nucl. Instrum. Methods Phys. Res., Sect. A 620, 227 (2010). [16] A. Aab et al. (Pierre Auger Collaboration), Data-driven estimation of the invisible energy of cosmic ray showers with the Pierre Auger Observatory, Phys. Rev. D 100, 082003 (2019). [17] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.126.152002 for additional information on the energy measurement, the comparison of the full model distribution for Rμin the individual energy bins, updated versions of select figures from [5], the predictions by the individual models, and additional information on the average number of muons and the fluctuations. [18] H. M. J. Barbosa, F. Catalani, J. A. Chinellato, and C. Dobrigkeit, Determination of the calorimetric energy in extensive air showers, Astropart. Phys. 22, 159 (2004). [19] M. Risse and D. Heck, Energy release in air showers, Astropart. Phys. 20, 661 (2004). [20] T. Pierog, M. Alekseeva, T. Bergmann, V. Chernatkin, R. Engel, D. Heck, N. N. Kalmykov, S. Ostapchenko, M. Unger, and K. Werner, Dependence of the longitudinal shower profile on the characteristics of hadronic multiparticle production, Int. Cosmic Ray Conf. 7, 103 (2005), PHYSICAL REVIEW LETTERS 126, 152002 (2021) 152002-9