Measurement of the fluctuations in the number of muons in extensive air showers with the Pierre Auger Observatory
Pierre Auger Collaboration, A. Aab, P. Abreu, M. Aglietta, J.M. Albury, I. Allekotte, A. Almela, J. Alvarez-Muñiz, R. Alves Batista, G.A. Anastasi, L. Anchordoqui, B. Andrada, S. Andringa, C. Aramo, P.R. Araújo Ferreira, H. Asorey, P. Assis, G. Avila, A.M. Badescu, A. Bakalova, A. Balaceanu, F. Barbato, R.J. Barreira Luz, K.H. Becker, J.A. Bellido, C. Berat, M.E. Bertaina, X. Bertou, P.L. Biermann, T. Bister, J. Biteau, J. Blazek, C. Bleve, M. Bohá?ová, D. Boncioli, C. Bonifazi, L. Bonneau Arbeletche, N. Borodai, A.M. Botti, J. Brack, T. Bretz, F.L. Briechle, P. Buchholz, A. Bueno, S. Buitink, M. Buscemi, K.S. Caballero-Mora, L. Caccianiga, A. Cancio, F. Canfora, I. Caracas, J.M. Carceller, R. Caruso, A. Castellina, F. Catalani, G. Cataldi, L. Cazon, M. Cerda, J.A. Chinellato, K. Choi, J. Chudoba, L. Chytka, R.W. Clay, A.C. Cobos Cerutti, R. Colalillo, A. Coleman, M.R. Coluccia, R. Conceição, A. Condorelli, G. Consolati, F. Contreras, F. Convenga, C.E. Covault, S. Dasso, K. Daumiller, B.R. Dawson, J.A. Day, R.M. de Almeida, J. de Jesús, S.J. de Jong, G. De Mauro, J.R.T. de Mello Neto, I. De Mitri, J. de Oliveira, D. de Oliveira Franco, V. de Souza, E. De Vito, J. Debatin, M. del Río, O. Deligny, H. Dembinski, N. Dhital, A. Di Matteo, C. Dobrigkeit, J.C. D'Olivo, R.C. dos Anjos, M.T. Dova, J. Ebr, R. Engel, I. Epicoco, et al. (268 additional authors not shown)
MMeasurement of the fluctuations in the number of muons in extensive air showerswith the Pierre Auger Observatory
A. Aab, P. Abreu, M. Aglietta,
3, 4
J.M. Albury, I. Allekotte, A. Almela,
7, 8
J. Alvarez-Mu˜niz, R. Alves Batista, G.A. Anastasi,
10, 4
L. Anchordoqui, B. Andrada, S. Andringa, C. Aramo, P.R. Ara´ujo Ferreira, H. Asorey, P. Assis, G. Avila, A.M. Badescu, A. Bakalova, A. Balaceanu, F. Barbato,
18, 12
R.J. Barreira Luz, K.H. Becker, J.A. Bellido, C. Berat, M.E. Bertaina,
10, 4
X. Bertou, P.L. Biermann, T. Bister, J. Biteau, J. Blazek, C. Bleve, M. Boh´aˇcov´a, D. Boncioli,
23, 24
C. Bonifazi, L. Bonneau Arbeletche, N. Borodai, A.M. Botti, J. Brack, T. Bretz, F.L. Briechle, P. Buchholz, A. Bueno, S. Buitink, M. Buscemi,
32, 33
K.S. Caballero-Mora, L. Caccianiga,
35, 36
A. Cancio,
8, 7
F. Canfora,
1, 37
I. Caracas, J.M. Carceller, R. Caruso,
32, 33
A. Castellina,
3, 4
F. Catalani, G. Cataldi, L. Cazon, M. Cerda, J.A. Chinellato, K. Choi, J. Chudoba, L. Chytka, R.W. Clay, A.C. Cobos Cerutti, R. Colalillo,
18, 12
A. Coleman, M.R. Coluccia,
46, 39
R. Concei¸c˜ao, A. Condorelli,
47, 24
G. Consolati,
36, 48
F. Contreras, F. Convenga,
46, 39
C.E. Covault,
49, 50
S. Dasso,
51, 52
K. Daumiller, B.R. Dawson, J.A. Day, R.M. de Almeida, J. de Jes´us,
7, 53
S.J. de Jong,
1, 37
G. De Mauro,
1, 37
J.R.T. de Mello Neto,
25, 55
I. De Mitri,
47, 24
J. de Oliveira, D. de Oliveira Franco, V. deSouza, E. De Vito,
46, 39
J. Debatin, M. del R´ıo, O. Deligny, N. Dhital, A. Di Matteo, C. Dobrigkeit, J.C. D’Olivo, R.C. dos Anjos, M.T. Dova, J. Ebr, R. Engel,
57, 53
I. Epicoco,
46, 39
M. Erdmann, C.O. Escobar, A. Etchegoyen,
7, 8
H. Falcke,
1, 63, 37
J. Farmer, G. Farrar, A.C. Fauth, N. Fazzini, F. Feldbusch, F. Fenu,
10, 4
B. Fick, J.M. Figueira, A. Filipˇciˇc,
69, 70
T. Fodran, M.M. Freire, T. Fujii,
64, 72
A. Fuster,
7, 8
C. Galea, C. Galelli,
35, 36
B. Garc´ıa, A.L. Garcia Vegas, H. Gemmeke, F. Gesualdi,
7, 53
A. Gherghel-Lascu, P.L. Ghia, U. Giaccari, M. Giammarchi, M. Giller, J. Glombitza, F. Gobbi, F. Gollan, G. Golup, M. G´omez Berisso, P.F. G´omez Vitale, J.P. Gongora, N. Gonz´alez, I. Goos,
6, 53
D. G´ora, A. Gorgi,
3, 4
M. Gottowik, T.D. Grubb, F. Guarino,
18, 12
G.P. Guedes, E. Guido,
4, 10
S. Hahn,
53, 7
R. Halliday, M.R. Hampel, P. Hansen, D. Harari, V.M. Harvey, A. Haungs, T. Hebbeker, D. Heck, G.C. Hill, C. Hojvat, J.R. H¨orandel,
1, 37
P. Horvath, M. Hrabovsk´y, T. Huege,
53, 31
J. Hulsman,
7, 53
A. Insolia,
32, 33
P.G. Isar, J.A. Johnsen, J. Jurysek, A. K¨a¨ap¨a, K.H. Kampert, B. Keilhauer, J. Kemp, H.O. Klages, M. Kleifges, J. Kleinfeller, M. K¨opke, G. Kukec Mezek, B.L. Lago, D. LaHurd, R.G. Lang, N. Langner, M.A. Leigui de Oliveira, V. Lenok, A. Letessier-Selvon, I. Lhenry-Yvon, D. Lo Presti,
32, 33
L. Lopes, R. L´opez, R. Lorek, Q. Luce, A. Lucero, J.P. Lundquist, A. Machado Payeras, G. Mancarella,
46, 39
D. Mandat, B.C. Manning, J. Manshanden, P. Mantsch, S. Marafico, A.G. Mariazzi, I.C. Mari¸s, G. Marsella,
46, 39
D. Martello,
46, 39
H. Martinez, O. Mart´ınezBravo, M. Mastrodicasa,
23, 24
H.J. Mathes, J. Matthews, G. Matthiae,
83, 84
E. Mayotte, P.O. Mazur, G. Medina-Tanco, D. Melo, A. Menshikov, K.-D. Merenda, S. Michal, M.I. Micheletti, L. Miramonti,
35, 36
S. Mollerach, F. Montanet, C. Morello,
3, 4
M. Mostaf´a, A.L. M¨uller,
7, 53
M.A. Muller,
41, 86, 25
K. Mulrey, R. Mussa, M. Muzio, W.M. Namasaka, L. Nellen, M. Niculescu-Oglinzanu, M. Niechciol, D. Nitz,
68, 87
D. Nosek, V. Novotny, L. Noˇzka, A Nucita,
46, 39
L.A. N´u˜nez, M. Palatka, J. Pallotta, P. Papenbreer, G. Parente, A. Parra, M. Pech, F. Pedreira, J. P¸ekala, R. Pelayo, J. Pe˜na-Rodriguez, J. Perez Armand, M. Perlin,
7, 53
L. Perrone,
46, 39
S. Petrera,
47, 24
T. Pierog, M. Pimenta, V. Pirronello,
32, 33
M. Platino, B. Pont, M. Pothast,
37, 1
P. Privitera, M. Prouza, A. Puyleart, S. Querchfeld, J. Rautenberg, D. Ravignani, M. Reininghaus,
53, 7
J. Ridky, F. Riehn, M. Risse, P. Ristori, V. Rizi,
23, 24
W. Rodrigues deCarvalho, J. Rodriguez Rojo, M.J. Roncoroni, M. Roth, E. Roulet, A.C. Rovero, P. Ruehl, S.J. Saffi, A. Saftoiu, F. Salamida,
23, 24
H. Salazar, G. Salina, J.D. Sanabria Gomez, F. S´anchez, E.M. Santos, E. Santos, F. Sarazin, R. Sarmento, C. Sarmiento-Cano, R. Sato, P. Savina,
46, 39, 58
C.M. Sch¨afer, V. Scherini, H. Schieler, M. Schimassek,
57, 7
M. Schimp, F. Schl¨uter,
53, 7
D. Schmidt, O. Scholten,
92, 31
P. Schov´anek, F.G. Schr¨oder,
45, 53
S. Schr¨oder, J. Schulte, S.J. Sciutto, M. Scornavacche,
7, 53
R.C. Shellard, G. Sigl, G. Silli,
7, 53
O. Sima,
17, 94
R. ˇSm´ıda, P. Sommers, J.F. Soriano, J. Souchard, R. Squartini, M. Stadelmaier,
53, 7
D. Stanca, S. Staniˇc, J. Stasielak, P. Stassi, A. Streich,
57, 7
M. Su´arez-Dur´an, T. Sudholz, T. Suomij¨arvi, A.D. Supanitsky, J. ˇSup´ık, Z. Szadkowski, A. Taboada, A. Tapia, C. Timmermans,
37, 1
O. Tkachenko, P. Tobiska, C.J. Todero Peixoto, B. Tom´e, A. Travaini, P. Travnicek, C. Trimarelli,
23, 24
M. Trini, M. Tueros, R. Ulrich, M. Unger, L. Vaclavek, M. Vacula, J.F. Vald´es Galicia, L. Valore,
18, 12
E. Varela, V. Varma K.C.,
7, 53
A. V´asquez-Ram´ırez, a r X i v : . [ h e p - e x ] F e b D. Veberiˇc, C. Ventura, I.D. Vergara Quispe, V. Verzi, J. Vicha, J. Vink, S. Vorobiov, H. Wahlberg, A.A. Watson, M. Weber, A. Weindl, L. Wiencke, H. Wilczy´nski, T. Winchen, M. Wirtz, D. Wittkowski, B. Wundheiler, A. Yushkov, O. Zapparrata, E. Zas, D. Zavrtanik,
70, 69
M. Zavrtanik,
69, 70
L. Zehrer, A. Zepeda, H. Dembinski, G. Torralba Elipe, and I. Vali˜no
47, 24 (The Pierre Auger Collaboration) ∗ IMAPP, Radboud University Nijmegen, Nijmegen, The Netherlands Laborat´orio de Instrumenta¸c˜ao e F´ısica Experimental de Part´ıculas – LIP andInstituto Superior T´ecnico – IST, Universidade de Lisboa – UL, Lisboa, Portugal Osservatorio Astrofisico di Torino (INAF), Torino, Italy INFN, Sezione di Torino, Torino, Italy University of Adelaide, Adelaide, S.A., Australia Centro At´omico Bariloche and Instituto Balseiro (CNEA-UNCuyo-CONICET), San Carlos de Bariloche, Argentina Instituto de Tecnolog´ıas en Detecci´on y Astropart´ıculas (CNEA, CONICET, UNSAM), Buenos Aires, Argentina Universidad Tecnol´ogica Nacional – Facultad Regional Buenos Aires, Buenos Aires, Argentina Instituto Galego de F´ısica de Altas Enerx´ıas (IGFAE),Universidade de Santiago de Compostela, Santiago de Compostela, Spain Universit`a Torino, Dipartimento di Fisica, Torino, Italy Department of Physics and Astronomy, Lehman College, City University of New York, Bronx, NY, USA INFN, Sezione di Napoli, Napoli, Italy RWTH Aachen University, III. Physikalisches Institut A, Aachen, Germany Observatorio Pierre Auger and Comisi´on Nacional de Energ´ıa At´omica, Malarg¨ue, Argentina University Politehnica of Bucharest, Bucharest, Romania Institute of Physics of the Czech Academy of Sciences, Prague, Czech Republic “Horia Hulubei” National Institute for Physics and Nuclear Engineering, Bucharest-Magurele, Romania Universit`a di Napoli “Federico II”, Dipartimento di Fisica “Ettore Pancini”, Napoli, Italy Bergische Universit¨at Wuppertal, Department of Physics, Wuppertal, Germany Univ. Grenoble Alpes, CNRS, Grenoble Institute of EngineeringUniv. Grenoble Alpes, LPSC-IN2P3, 38000 Grenoble, France, France Max-Planck-Institut f¨ur Radioastronomie, Bonn, Germany Universit´e Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France, France Universit`a dell’Aquila, Dipartimento di Scienze Fisiche e Chimiche, L’Aquila, Italy INFN Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy Universidade Federal do Rio de Janeiro, Instituto de F´ısica, Rio de Janeiro, RJ, Brazil Universidade de S˜ao Paulo, Instituto de F´ısica, S˜ao Paulo, SP, Brazil Institute of Nuclear Physics PAN, Krakow, Poland Colorado State University, Fort Collins, CO, USA Universit¨at Siegen, Fachbereich 7 Physik – Experimentelle Teilchenphysik, Siegen, Germany Universidad de Granada and C.A.F.P.E., Granada, Spain Vrije Universiteit Brussels, Brussels, Belgium Universit`a di Catania, Dipartimento di Fisica e Astronomia, Catania, Italy INFN, Sezione di Catania, Catania, Italy Universidad Aut´onoma de Chiapas, Tuxtla Guti´errez, Chiapas, M´exico Universit`a di Milano, Dipartimento di Fisica, Milano, Italy INFN, Sezione di Milano, Milano, Italy Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF), Science Park, Amsterdam, The Netherlands Universidade de S˜ao Paulo, Escola de Engenharia de Lorena, Lorena, SP, Brazil INFN, Sezione di Lecce, Lecce, Italy Observatorio Pierre Auger, Malarg¨ue, Argentina Universidade Estadual de Campinas, IFGW, Campinas, SP, Brazil Universit´e Libre de Bruxelles (ULB), Brussels, Belgium Palacky University, RCPTM, Olomouc, Czech Republic Instituto de Tecnolog´ıas en Detecci´on y Astropart´ıculas (CNEA, CONICET, UNSAM),and Universidad Tecnol´ogica Nacional – Facultad Regional Mendoza (CONICET/CNEA), Mendoza, Argentina University of Delaware, Department of Physics and Astronomy, Bartol Research Institute, Newark, DE, USA Universit`a del Salento, Dipartimento di Matematica e Fisica “E. De Giorgi”, Lecce, Italy Gran Sasso Science Institute, L’Aquila, Italy Politecnico di Milano, Dipartimento di Scienze e Tecnologie Aerospaziali , Milano, Italy Case Western Reserve University, Cleveland, OH, USA also at Radboud Universtiy Nijmegen, Nijmegen, The Netherlands Instituto de Astronom´ıa y F´ısica del Espacio (IAFE, CONICET-UBA), Buenos Aires, Argentina Departamento de F´ısica and Departamento de Ciencias de la Atm´osfera y los Oc´eanos,FCEyN, Universidad de Buenos Aires and CONICET, Buenos Aires, Argentina Karlsruhe Institute of Technology, Institut f¨ur Kernphysik, Karlsruhe, Germany Universidade Federal Fluminense, EEIMVR, Volta Redonda, RJ, Brazil Universidade Federal do Rio de Janeiro (UFRJ),Observat´orio do Valongo, Rio de Janeiro, RJ, Brazil Universidade de S˜ao Paulo, Instituto de F´ısica de S˜ao Carlos, S˜ao Carlos, SP, Brazil Karlsruhe Institute of Technology, Institute for Experimental Particle Physics (ETP), Karlsruhe, Germany CNRS/IN2P3, IJCLab, Universit´e Paris-Saclay, Orsay, France Universidad Nacional Aut´onoma de M´exico, M´exico, D.F., M´exico Universidade Federal do Paran´a, Setor Palotina, Palotina, Brazil IFLP, Universidad Nacional de La Plata and CONICET, La Plata, Argentina Fermi National Accelerator Laboratory, USA Stichting Astronomisch Onderzoek in Nederland (ASTRON), Dwingeloo, The Netherlands University of Chicago, Enrico Fermi Institute, Chicago, IL, USA New York University, New York, NY, USA Fermi National Accelerator Laboratory, Fermilab, Batavia, IL, USA Karlsruhe Institute of Technology, Institut f¨ur Prozessdatenverarbeitung und Elektronik, Karlsruhe, Germany Michigan Technological University, Houghton, MI, USA Experimental Particle Physics Department, J. Stefan Institute, Ljubljana, Slovenia Center for Astrophysics and Cosmology (CAC),University of Nova Gorica, Nova Gorica, Slovenia Instituto de F´ısica de Rosario (IFIR) – CONICET/U.N.R. and Facultadde Ciencias Bioqu´ımicas y Farmac´euticas U.N.R., Rosario, Argentina now at Hakubi Center for Advanced Research and Graduate School of Science, Kyoto University, Kyoto, Japan University of (cid:32)L´od´z, Faculty of Astrophysics, (cid:32)L´od´z, Poland Universidade Estadual de Feira de Santana, Feira de Santana, Brazil Institute of Space Science, Bucharest-Magurele, Romania Colorado School of Mines, Golden, CO, USA Centro Federal de Educa¸c˜ao Tecnol´ogica Celso Suckow da Fonseca, Nova Friburgo, Brazil Universidade Federal do ABC, Santo Andr´e, SP, Brazil Laboratoire de Physique Nucl´eaire et de Hautes Energies (LPNHE),Sorbonne Universit´e, Universit´e de Paris, CNRS-IN2P3, Paris, France Benem´erita Universidad Aut´onoma de Puebla, Puebla, M´exico Universit¨at Hamburg, II. Institut f¨ur Theoretische Physik, Hamburg, Germany Louisiana State University, Baton Rouge, LA, USA Universit`a di Roma “Tor Vergata”, Dipartimento di Fisica, Roma, Italy INFN, Sezione di Roma “Tor Vergata”, Roma, Italy Pennsylvania State University, University Park, PA, USA also at Universidade Federal de Alfenas, Po¸cos de Caldas, Brazil also at Karlsruhe Institute of Technology, Karlsruhe, Germany Charles University, Faculty of Mathematics and Physics,Institute of Particle and Nuclear Physics, Prague, Czech Republic Universidad Industrial de Santander, Bucaramanga, Colombia Centro de Investigaciones en L´aseres y Aplicaciones,CITEDEF and CONICET, Villa Martelli, Argentina Unidad Profesional Interdisciplinaria en Ingenier´ıa y Tecnolog´ıas Avanzadasdel Instituto Polit´ecnico Nacional (UPIITA-IPN), M´exico, D.F., M´exico KVI – Center for Advanced Radiation Technology,University of Groningen, Groningen, The Netherlands Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, RJ, Brazil also at University of Bucharest, Physics Department, Bucharest, Romania University of (cid:32)L´od´z, Faculty of High-Energy Astrophysics,(cid:32)L´od´z, Poland Universidad de Medell´ın, Medell´ın, Colombia Universiteit van Amsterdam, Faculty of Science, Amsterdam, The Netherlands School of Physics and Astronomy, University of Leeds, Leeds, United Kingdom Centro de Investigaci´on y de Estudios Avanzados del IPN (CINVESTAV), M´exico, D.F., M´exico (Dated: February 17, 2021)We present the first measurement of the fluctuations in the number of muons in extensive airshowers produced by ultra-high energy cosmic rays. We find that the measured fluctuations arein good agreement with predictions from air shower simulations. This observation provides newinsights into the origin of the previously reported deficit of muons in air shower simulations andconstrains models of hadronic interactions at ultra-high energies. Our measurement is compatiblewith the muon deficit originating from small deviations in the predictions from hadronic interactionmodels of particle production that accumulate as the showers develop.
INTRODUCTION
Ultra High Energy Cosmic Rays (UHECRs) are par-ticles coming from outer space, with energies exceeding10 eV. They provide the only experimental opportunityto explore particle physics beyond energies reachable byEarth-based accelerators, which go up to cosmic ray en-ergies of 9 × eV.The Pierre Auger Observatory [1] detects extensive airshowers that are initiated by the UHECRs colliding withthe nuclei in the atmosphere. Information about UHE-CRs is extracted using simulations based on hadronicinteraction models which rely on extrapolations of ac-celerator measurements to unexplored regions of phasespace, most notably the forward and highest-energy re-gion. In addition, accelerator experiments at the highestenergies either probe the interactions between protonsor of protons with heavy nuclei, while most interactionswithin air showers are between pions and light nuclei.A further challenge is that the UHECR mass has tobe measured despite not being yet completely decou-pled from the hadronic uncertainties. The observablewith the least dependence on hadronic interactions isthe atmospheric depth at which the longitudinal de-velopment of the electromagnetic (EM) component ofthe shower reaches the maximum number of particles,namely X max [2].In hadronic cascades the energy of each interactingparticle is distributed among the secondaries, mostly pi-ons. Neutral pions rapidly decay into two photons feedinga practically decoupled electromagnetic cascade (otherresonances decaying into π ’s, electrons and or photonsalso contribute). Charged pions (and other long livedmesons like kaons) tend to further interact until their in-dividual energies are below a critical value, below whichthey are more likely to decay. Muons, which are productsof hadronic decays, are thus predominantly produced inthe final shower stages. In sufficiently inclined showers,the pure EM component is absorbed in the atmosphereand the particles that reach the ground (muons and muondecay products) directly sample the muon content [3, 4],reflecting the hadronic component of the shower.Air showers are mainly detected at the Pierre AugerObservatory by the Surface Detector (SD), an array ofwater-Cherenkov detector stations, and the FluorescenceDetector (FD) consisting of 24 fluorescence telescopes.By selecting the sub-sample of events reconstructed withboth the SD and the FD and with zenith angles exceeding62 ◦ , both the muon content and the energy of the showerare simultaneously measured.The results obtained indicate that all the simulationsunderestimate the number of muons in the showers [5, 6].These analyses come with the caveat that they cannotdistinguish a muon rescaling from a shift in the absoluteenergy scale of the FD measurement. However, muon content and energy scale were disentangled in a comple-mentary technique based on showers with zenith anglesbelow 60 ◦ . Using the longitudinal profile of the showerin the atmosphere obtained with the FD and the sig-nals at the ground measured with the SD, it was shownthat the muonic component still has to be scaled up tomatch observed data, while no rescaling of the EM com-ponent and the FD energy is required [7]. The measure-ments with the FD also show that both the position ofthe shower maximum in the atmosphere ( X max ) and theentire shape of the EM shower are well described by thesimulations [8, 9]. At lower energies, down to ∼ . eV,in a measurement using the sub-array of buried scintil-lators of the Pierre Auger Observatory, a direct countof the muons independent of EM contamination was ob-tained, which also shows that simulations produce toofew muons [10]. There is much evidence that all thesimulations underpredict the average number of muonsin the showers: a comprehensive study of muon num-ber measurements made with different experiments hasshown that the muon deficit in simulations starts around ∼ eV and steadily increases with energy. Dependingon model and experiment, the deficit at ∼ eV rangesbetween tens of percent up to a factor of two [11].The increased statistics obtained at the Pierre AugerObservatory allows us to now take a further step and ex-plore fluctuations in the number of muons between show-ers, hereinafter referred to as physical fluctuations . Theratio of the physical fluctuations to the average numberof muons (relative fluctuations) has been shown to bemostly dominated by the first interaction, rather thanthe lower energy interactions deeper in the shower de-velopment [12, 13]. Here, we exploit the sensitivity offluctuations to the first-interaction to explore hadronicinteractions well above the energies achievable in accel-erator experiments. METHODOLOGY
Our analysis here is based on the set of inclined airshowers (62 ◦ < θ < ◦ ) that are reconstructed bothwith the SD and with the FD between 1 January 2004and 31 December 2017. For each event we obtain in-dependent measurements of the muon content (with theSD) and the calorimetric energy (with the FD). To ensurethe showers can be reconstructed with small uncertain-ties, we select only events with at least four triggered sta-tions in the SD array and we further require that all thestations surrounding the impact point of the shower onthe ground are operational at the time of the event. Onlyevents with good atmospheric conditions, (few clouds anda low aerosol content) are accepted in order to guaranteea good energy reconstruction with the FD. In addition itis required that the entire shower profile and in particular X max is within the field of view of our telescopes. Sinceheavy primaries penetrate the atmosphere less than lightones, the acceptance with this selection would be massdependent. To avoid this bias we constrain the field ofview to the region where all values of X max are accepted.Further details are given in [5, 14]. These selection crite-ria result in a total number of events of 786.In addition only events with energy larger than 4 × 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 2-Dmodel of the lateral profile of the muon density at theground to the observed signals in the SD array. The freeparameters of the fit are the zenith and azimuth anglesof the shower, the impact point of the shower on theground (shower core position) and a normalization fac-tor with respect to a reference muon density profile insimulated proton showers at 10 eV [3]. There exists aresidual pure EM component in showers with low zenithangles and stations very close to the shower core position(at 400 m and 64 ◦ it is ∼ R µ , which is given by theintegrated number of muons at the ground divided by areference given by the average number of muons in sim-ulated proton showers at 10 eV and the given zenithangle. At 10 eV and an inclination of 60 ◦ , R µ = 1 cor-responds to 2 . × 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, E cal , is re-constructed by integrating the longitudinal shower pro-files observed with the FD [9, 15]. The total energy ofthe shower is then obtained by adding the average energycarried away by muons and neutrinos, the so-called invis-ible energy, E = E cal + E inv . At 10 eV E inv accountsfor 14% of the total energy in air showers [16–20].In Fig. 1 the muon number, R µ , is shown as a func-tion of the measured energy. Markers on the top of theframe define the bins in energy for which we will extractthe fluctuations, with the numbers of events in each binshown above. The bins are chosen such that the num-ber of events in each is similar. Based on models of airshower development and given the gradual change of thecomposition in this energy range (single logarithmic de-pendence on energy) [8, 21–23], the number of muons isrelated to the primary energy by a single power-law (cid:104) R µ (cid:105) ( E ) = a ( E/ (10 eV)) b , (1)which can be fitted following a procedure described inthe text below. The best-fit parameters are given atthe beginning of the next section. The scattering in thedata has three sources: experimental uncertainties in theenergy s E and in the muon number s µ from event re-construction (both represented by the error bars), andthe physical fluctuations in the muon number denoted as E / eV10 R µ
44 46 48 47 47 49 σ sys ( E ) D/ n . d . f . = 289 . / a = 1 . ± . b = 0 . ± . data ± σ res (281 events) h R µ i ( E ) FIG. 1. Number of muons as a function of the measuredenergy. The black line is the fitted (cid:104) R µ (cid:105) = a ( E/ (10 eV)) b .Markers on the top of the frame define the bins in which thefluctuations are evaluated. The numbers give the events ineach bin. The effect of the uncertainty of the absolute energyscale is indicated by σ sys ( E ). σ . Given Eq. (1) the variance of the muon number is σ + b (cid:104) s E (cid:105) + (cid:104) s µ (cid:105) .In this work, we adopt a method based on maximiz-ing the likelihood of a probability distribution function(PDF). The PDF incorporates the various contributionsto the fluctuations, treating each energy bin indepen-dently while also accounting for the effect of the migra-tion of events between bins [5, 24].The model assumes that measurements of E and R µ follow Gaussian distributions centered at the true value,with widths given by the detector resolution s E and s µ ,which are the uncertainties obtained in each individualevent reconstruction [3, 25].Physical fluctuations are also assumed to follow aGaussian distribution of width σ . Simulations haveshown this is an acceptable approximation given theevent number in each bin.The total PDF is obtained through the convolutionof the detector response and the physical fluctuationswith the probability distribution of the hybrid eventsmeasured at the Pierre Auger Observatory. The log-likelihood function is then given byln L ( a, b, ˆ σ , .., ˆ σ ) = (cid:88) i ln (cid:88) k =0 E k (cid:90) E k − d E h ( E ) C ( E ) × exp (cid:18) −
12 ( E i − E ) s E (cid:19) × exp (cid:18) −
12 ( R µ,i − (cid:104) R µ (cid:105) ( E ) ) s µ + (ˆ σ k · (cid:104) R µ (cid:105) ( E )) (cid:19)(cid:21) . (2)The probability of hybrid events, h ( E ), (product of theenergy spectrum of CRs and the efficiency of detection),can be obtained from the data itself as explained in [24]and [10, 26]. The RHS of Eq. (2) depends on the param-eters a, b via Eq. (1). To obtain the energy dependenceof the fluctuations we parameterize σ by six independentvalues such that σ ( E ) = ˆ σ k · (cid:104) R µ (cid:105) ( E ) where the con-stants ˆ σ k are the relative fluctuations in the k -th energybin with limits [ E k − , E k ] where k runs from one to six.In Eq. (2), k = 0 corresponds to the contributions fromthe interval [0 , E thr ] where the SD is not fully efficient.The fluctuations here are assumed to take the value ofthe first fitted bin, ˆ σ ≡ ˆ σ .The sum over the index i in Eq. (2) (the usual sum overthe log-likelihoods of events), includes only events abovethe energy threshold of 4 × eV. The function C ( E ) isthe normalization factor from the double Gaussian. Theresult of the fit for the parameters a and b are shownin Fig. 1. The fluctuations are shown in Fig. 2. Thedistribution of the number of muons and the PDF in theindividual energy bins can be found in the SupplementalMaterial [27].The dominant systematic uncertainties of σ come fromthe uncertainties in the resolutions s E and s µ . For s µ weestimate the uncertainty using simulations and data. Insimulations the uncertainty was estimated by the spreadin a sample of simulated showers where each shower is re-constructed multiple times, each time changing only theimpact point at the ground. For data we reconstruct thesame event multiple times, leaving out the signals fromone of the detector stations. The average relative reso-lution, (cid:104) s µ /R µ (cid:105) , and its systematic uncertainty is thus(10 ±
3) % at 10 eV.We verified the values of s E by studying the differencein the energy reconstruction of events measured inde-pendently by two or more FD stations. The width ofthe distribution of these energy-differences is found tobe compatible with s E . We therefore take the statisti-cal one- σ uncertainties of this cross check as a conserva-tive upper limit of the systematic uncertainty of s E [28].The average relative energy resolution, (cid:104) s E /E (cid:105) , is about(8 . ± .
9) % at 10 eV.We have further confirmed that there are no significantcontributions to the fluctuations from differences between TABLE I. Contributions to the systematic uncertainty in therelative fluctuations around 10 eV (10 . eV to 10 . eV).The central value is σ/ (cid:104) R µ (cid:105) = 0 . ± .
029 (stat . ) ± .
007 (syst . ).Source of uncertainty Uncertainty E absolute scale (cid:104) E (cid:105) < . E resolution s E R µ absolute scale (cid:104) R µ (cid:105) R µ resolution s µ R µ azimuthal modulation (cid:104) R µ (cid:105) ( φ ) 0.5 %Total systematics 7.0 % the individual FD stations, neither related to the longtime performance evolution of the SD and FD detectors.Any residual electromagnetic component in the signalwould affect the lower zenith angles more. We thereforesplit the event sample at the median zenith angle (66 ◦ )and compare the resulting fluctuations. We find no sig-nificant difference between the more and the less inclinedsample.In another test we do find a small modulation of (cid:104) R µ (cid:105) with the azimuth angle ( < R µ and E ,the systematic uncertainty on σ due to the uncertainty ofthe absolute energy scale of 14 % [25] practically cancelsout in the relative fluctuations. The systematic uncer-tainty in the absolute scale of R µ of 11 % [5] drops outfor the same reason. The systematic effects for the binaround 10 eV are summarized in Tab. I. Over all ener-gies the systematic uncertainties are below 8%. RESULTS AND DISCUSSION
The best-fit value for the average relative numberof muons at 10 eV (parameter a ) is (cid:104) R µ (cid:105) (10 eV) =1 . ± .
02 (stat . ) +0 . − . (syst . ). For the slope (pa-rameter b ) we find d (cid:104) ln R µ (cid:105) / d ln E = 0 . ± .
02 (stat . ) +0 . − . (syst . ). These values are consistent withthe values previously reported [5, 29].The measured relative fluctuations as a function of theenergy are shown in Fig. 2. We note that the measure-ment falls within the range that is expected from cur-rent hadronic interaction models for pure proton andpure iron primaries [30–38]. To estimate the effect ofa mixed composition, we take the fractions of the fourmass components (proton, helium, nitrogen and iron) de-rived from the X max measurements [8, 39, 40] and, using E / eV0 . . . . . . . σ / h R µ i pFeEPOS-LHCQGSJetII-04SIBYLL-2.3d data4-mass- X max -fit+models10080 90 200 300 400Proton-proton equivalent CM energy √ s / TeV
FIG. 2. Measured relative fluctuations in the number ofmuons as a function of the energy and the predictions fromthree interaction models for proton (red) and iron (blue)showers. The gray band represents the expectations fromthe measured mass composition interpreted with the interac-tion models. The statistical uncertainty in the measurementis represented by the error bars. The total systematic uncer-tainty is indicated by the square brackets. the simulations of the pure primaries, calculate the cor-responding fluctuations in the number of muons. Thegray band in Fig. 2 encompasses the predicted σ/ (cid:104) R µ (cid:105) of the three interaction models Qgsjet
II-04,
Epos-lhc and
Sibyll eV, the relative fluctuations σ/ (cid:104) R µ (cid:105) againstthe average number of muons (cid:104) R µ (cid:105) . Given any one ofthe interaction models, any particular mixture of thefour components p, He, N and Fe falls somewhere withinone of the areas enclosed by the corresponding coloredlines. The points of pure composition in this contour arelabeled accordingly. For each model the expected val-ues for σ/ (cid:104) R µ (cid:105) and (cid:104) R µ (cid:105) given the composition mixtureobtained from the X max measurements [8] is indicatedwithin each contour by the correspondingly colored starmarker. The shaded areas surrounding the star markersindicate the statistical and systematic uncertainties in-herited from the X max measurements [42]. Finally ourmeasurement with statistical and systematic uncertaintyis shown by the black marker.Within the uncertainty none of the predictions fromthe interaction models and the X max -composition (starmarkers) are consistent with our measurement. Thepredictions from the interaction models Qgsjet
II-04,
Epos-lhc and
Sibyll . . . . . . h R µ i . . . . . . . . . σ / h R µ i p He N Fe E = 10 eV4-mass- X max -fit+modelEPOS-LHC QGSJetII-04 SIBYLL-2.3d data FIG. 3. Data (black, with error bars) compared to mod-els for the fluctuations and the average number of muons forshowers with a primary energy of 10 eV. Fluctuations areevaluated in the energy range from 10 . eV and 10 . eV.The statistical uncertainty is represented by the error bars.The total systematic uncertainty is indicated by the squarebrackets. The expectation from the interaction models forany mixture of the four components p, He, N, Fe is illus-trated by the colored contours. The values preferred by themixture derived from the X max measurements are indicatedby the star symbols. The shaded areas show the regions al-lowed by the statistical and systematic uncertainties of the X max measurement [42]. muons of 43%, 35%, and 26%, respectively. For the fluc-tuations, no rescaling is necessary for any model.Taken together, the average value and fluctuations ofthe muon flux constrain the way hadronic interactionmodels should be changed to agree with air shower data.To see this we briefly discuss the origin of the fluctua-tions.The average number of muons in a proton showerof energy E has been shown in simulations to scale as (cid:104) N ∗ µ (cid:105) = C E β where β (cid:39) . N µ = m (cid:88) j =1 C E βj = (cid:104) N ∗ µ (cid:105) m (cid:88) j =1 x βj = (cid:104) N ∗ µ (cid:105) α , (3)where index j runs over m secondary particles which rein-teract hadronically and x j = E j /E is the fraction of en-ergy fed to the hadronic shower by each [43]. In this ex-pression the fluctuations in N µ are induced by α in thefirst generation which fluctuates because the multiplicity m and the energies x j of the secondaries fluctuate [13].We can continue this reasoning for the subsequent gen-erations to obtain N µ (cid:104) N ∗ µ (cid:105) = α · α · . . . α i . . . α n , (4)here the subindex i runs over n generations, until thecascade stops. We note that for the calculation of α , inthe second generation, there are m particles contribut-ing. Assuming the distributions of the α ’s for each oneare similar, when adding up the muons produced by each,the fluctuations produced by one are statistically likely tobe compensated by another. In other words the α distri-bution is narrower by a factor ∼ / √ m . The deeper thegeneration, the sharper the corresponding α i is expectedto be. As a result, the dominant part of the fluctuationscomes from the first interaction. This has also been ob-served with simulations. The model can be generalisedfor primary nuclei with mass A using the superpositionmodel, and fixing the number of participants to A pro-tons, which reduce the different contributions to the fluc-tuations by a factor ∼ / √ A .There are two options to increase the average numberof muons in air showers. One is to increase α in a specificgeneration, notably the first where the energy is the high-est and exotic phenomena could conceivably play a role,i.e. α → α + δα . Note that if only the first generationis modified (implying some sort of threshold effect fornew physics) the increase in N µ is linear with the mod-ification. There are several examples in the literaturewhere this approach has been used assuming differentmechanisms [44–48]. For the fluctuations the change de-pends on the model. Alternatively, the number of muonscan be increased by introducing small deviations in thehadronic energy fraction, δα , in all generations. Accu-mulated along a number n of generations, these smalldeviations build up as N µ ∝ ( α + δα ) n . For instance,a 5% deviation per generation converts into ∼
30% de-viation after six generations [49]. On the other hand, achange of 5% in the fluctuations of α is not amplified inthe muon fluctuations because of the suppression in latergenerations. This approach characterises the increase inthe number of muons in the current hadronic interac-tion models with regard to previous models. It is alsocompatible with the increase of the discrepancy in theaverage number of muons across a wide range of energiesreported in [11].The present analysis finding that fluctuations are con-sistent with model predictions means that the increase inmuon number may be a small effect accumulating overmany generations, or a very particular modification ofthe first interaction that changes N µ without changingthe fluctuations [50]. SUMMARY
We have presented for the first time a measurementof the fluctuations in the number of muons in inclinedair showers, as a function of the UHECR primary en-ergy. Within the current uncertainties, the relative fluc-tuations show no discrepancy with respect to the ex-pectation from current high-energy hadronic interactionmodels and the composition taken from X max measure-ments. This agreement between models and data for thefluctuations, combined with the significant deficit in thepredicted total number of muons, points to the origin ofthe models’ muon deficit being a small deficit at everystage of the shower which accumulates along the showerdevelopment, rather than a discrepancy in the first in-teraction. Adjustments to models to address the currentmuon deficit, must therefore not alter the predicted rel-ative fluctuations.The Pierre Auger Observatory is currently undergoingan upgrade which includes the deployment of scintilla-tors on top of the SD stations [51] to help disentanglethe muonic and electromagnetic content of the showers,as well as an array of radio antennas [52]. It has beenshown that radio arrays can provide an estimate of thecalorimetric energy [53] and therefore it will soon be pos-sible to perform an analysis similar to the one presentedhere with much larger statistics using hybrid events mea-sured by the high-duty-cycle radio and surface detectorarrays [52]. ACKNOWLEDGMENTS
The successful installation, commissioning, and oper-ation of the Pierre Auger Observatory would not havebeen possible without the strong commitment and effortfrom the technical and administrative staff in Malarg¨ue.We are very grateful to the following agencies and orga-nizations for financial support:Argentina – Comisi´on Nacional de Energ´ıa At´omica;Agencia Nacional de Promoci´on Cient´ıfica y Tecnol´ogica(ANPCyT); Consejo Nacional de InvestigacionesCient´ıficas y T´ecnicas (CONICET); Gobierno de laProvincia de Mendoza; Municipalidad de Malarg¨ue;NDM Holdings and Valle Las Le˜nas; in gratitude fortheir continuing cooperation over land access; Australia– the Australian Research Council; Brazil – ConselhoNacional de Desenvolvimento Cient´ıfico e Tecnol´ogico(CNPq); Financiadora de Estudos e Projetos (FINEP);Funda¸c˜ao de Amparo `a Pesquisa do Estado de Rio deJaneiro (FAPERJ); S˜ao Paulo Research Foundation(FAPESP) Grants No. 2019/10151-2, No. 2010/07359-6and No. 1999/05404-3; Minist´erio da Ciˆencia, Tecnologia,Inova¸c˜oes e Comunica¸c˜oes (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 andCZ.02.1.01/0.0/0.0/17 049/0008422; France – Cen-tre de Calcul IN2P3/CNRS; Centre National de laRecherche Scientifique (CNRS); Conseil R´egionalIle-de-France; D´epartement Physique Nucl´eaire et Cor-pusculaire (PNC-IN2P3/CNRS); D´epartement Sciencesde l’Univers (SDU-INSU/CNRS); Institut Lagrangede Paris (ILP) Grant No. LABEX ANR-10-LABX-63 within the Investissements d’Avenir ProgrammeGrant No. ANR-11-IDEX-0004-02; Germany – Bun-desministerium f¨ur Bildung und Forschung (BMBF);Deutsche Forschungsgemeinschaft (DFG); Finanzmin-isterium Baden-W¨urttemberg; Helmholtz Alliance forAstroparticle Physics (HAP); Helmholtz-GemeinschaftDeutscher Forschungszentren (HGF); Ministerium f¨urInnovation, Wissenschaft und Forschung des LandesNordrhein-Westfalen; Ministerium f¨ur Wissenschaft,Forschung und Kunst des Landes Baden-W¨urttemberg;Italy – Istituto Nazionale di Fisica Nucleare (INFN);Istituto Nazionale di Astrofisica (INAF); Ministerodell’Istruzione, dell’Universit´a e della Ricerca (MIUR);CETEMPS Center of Excellence; Ministero degliAffari Esteri (MAE); M´exico – Consejo Nacional deCiencia y Tecnolog´ıa (CONACYT) No. 167733; Uni-versidad Nacional Aut´onoma de M´exico (UNAM);PAPIIT DGAPA-UNAM; The Netherlands – Min-istry of Education, Culture and Science; NetherlandsOrganisation for Scientific Research (NWO); Dutchnational e-infrastructure with the support of SURFCooperative; Poland -Ministry of Science and HigherEducation, grant No. DIR/WK/2018/11; NationalScience 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 fundswithin Programa Operacional Factores de Competitivi-dade through Funda¸c˜ao para a Ciˆencia e a Tecnologia(COMPETE); Romania – Romanian Ministry of Edu-cation and Research, the Program Nucleu within MCI(PN19150201/16N/2019 and PN19060102) and projectPN-III-P1-1.2-PCCDI-2017-0839/19PCCDI/2018 withinPNCDI III; Slovenia – Slovenian Research Agency,grants P1-0031, P1-0385, I0-0033, N1-0111; Spain –Ministerio de Econom´ıa, Industria y Competitividad(FPA2017-85114-P and FPA2017-85197-P), Xuntade Galicia (ED431C 2017/07), Junta de Andaluc´ıa(SOMM17/6104/UGR), Feder Funds, RENATA RedNacional Tem´atica de Astropart´ıculas (FPA2015-68783-REDT) and Mar´ıa de Maeztu Unit of Excellence (MDM-2016-0692); USA – Department of Energy, ContractsNo. DE-AC02-07CH11359, No. DE-FR02-04ER41300,No. DE-FG02-99ER41107 and No. DE-SC0011689;National Science Foundation, Grant No. 0450696; TheGrainger Foundation; Marie Curie-IRSES/EPLANET;European Particle Physics Latin American Network;and UNESCO. ∗ et al. (Pierre Auger Collaboration), “ThePierre Auger Cosmic Ray Observatory,” Nucl. Instrum.Meth. A , 172–213 (2015), arXiv:1502.01323 [astro-ph.IM].[2] J. Linsley, “Structure of large showers at depth834 g / cm , applications,” (1977), in Proceedings of the15th Int. Cosmic Ray Conf., Plovdiv, Bulgaria, vol. 12,p. 89.[3] Alexander Aab et al. (Pierre Auger Collaboration), “Re-construction of inclined air showers detected with thePierre Auger Observatory,” JCAP , 019 (2014),arXiv:1407.3214 [astro-ph.HE].[4] I. Valino, J. Alvarez-Muniz, M. Roth, R.A. Vazquez,and E. Zas, “Characterisation of the electromagneticcomponent in ultra-high energy inclined air showers,”Astropart. Phys. , 304–317 (2010), arXiv:0910.2873[astro-ph.HE].[5] Alexander Aab et al. (Pierre Auger Collaboration),“Muons in air showers at the Pierre Auger Observa-tory: Mean number in highly inclined events,” Phys.Rev. D , 032003 (2015), [Erratum: Phys. Rev. D91,no.5,059901(2015)], arXiv:1408.1421 [astro-ph.HE].[6] S. J. Sciutto, “Air showers, hadronic models, and muonproduction,” Proceedings, Ultra High Energy CosmicRays (UHECR 2018): Paris, France, October 8-12, 2018 ,EPJ Web Conf. , 02007 (2019), arXiv:1904.12056[astro-ph.HE].[7] Alexander Aab et al. (Pierre Auger Collaboration),“Testing Hadronic Interactions at Ultrahigh Energieswith Air Showers Measured by the Pierre AugerObservatory,” Phys. Rev. Lett. , 192001 (2016),arXiv:1610.08509 [hep-ex].[8] Alexander Aab et al. (Pierre Auger Collaboration),“Depth of maximum of air-shower profiles at thePierre Auger Observatory. I. Measurements at energiesabove 10 . eV,” Phys. Rev. D , 122005 (2014),arXiv:1409.4809 [astro-ph.HE].[9] Alexander Aab et al. (Pierre Auger Collaboration),“Measurement of the average shape of longitudinal pro-files of cosmic-ray air showers at the Pierre Auger Obser-vatory,” JCAP , 018 (2019), arXiv:1811.04660 [astro-ph.HE].[10] A. Aab et al. (Pierre Auger Collaboration), “Direct mea-surement of the muonic content of extensive air showersbetween × and × eV at the Pierre AugerObservatory,” Eur. Phys. J. C , 751 (2020).[11] Lorenzo Cazon et al. (EAS-MSU1, IceCube, KASCADE-Grande, NEVOD-DECOR, Pierre Auger, SUGAR, Tele-scope Array, and Yakutsk EASArray Collaborations),“Working Group Report on the Combined Analysis ofMuon Density Measurements from Eight Air ShowerEx-periments,” Proceedings for the 36th International Cos-mic Ray Conference (ICRC 2019) , PoS
ICRC2019 , 214(2019).[12] S. Fukui, H. Hasegawa, T. Matano, I. Miura, M. Oda,K. Suga, G. Tanahashi, and Y. Tanaka, “A Study onthe Structure of the Extensive Air Shower,” Progress ofTheoretical Physics Supplement , 1–53 (1960).[13] Lorenzo Cazon, Ruben Concei¸c˜ao, and Felix Riehn,“Probing the energy spectrum of hadrons in proton air in- teractions at ultrahigh energies through the fluctuationsof the muon content of extensive air showers,” Phys. Lett.B , 68–76 (2018), arXiv:1803.05699 [hep-ph].[14] P. Abreu, M. Aglietta, E.J. Ahn, D. Allard, I. Allekotte,J. Allen, J. Alvarez Castillo, J. Alvarez-Mu˜niz, M. Am-brosio, A. Aminaei, and et al., “The exposure of thehybrid detector of the Pierre Auger Observatory,” As-troparticle Physics , 368–381 (2011).[15] J. Abraham et al. (Pierre Auger Collaboration), “TheFluorescence Detector of the Pierre Auger Observa-tory,” Nucl. Instrum. Meth. A , 227–251 (2010),arXiv:0907.4282 [astro-ph.IM].[16] A. Aab et al. (Pierre Auger Collaboration), “Data-drivenestimation of the invisible energy of cosmic ray showerswith the Pierre Auger Observatory,” Phys. Rev. D ,082003 (2019), arXiv:1901.08040 [astro-ph.IM].[17] (), See Supplemental Material at [URL will be insertedby publisher] for additional information on the energymeasurement.[18] Henrique M.J. Barbosa, F. Catalani, J.A. Chinellato,and Carola Dobrigkeit, “Determination of the calorimet-ric energy in extensive air showers,” Astropart. Phys. ,159–166 (2004), arXiv:astro-ph/0310234.[19] Markus Risse and Dieter Heck, “Energy release in airshowers,” Astropart. Phys. , 661 (2004), arXiv:astro-ph/0308158.[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 longitudi-nal shower profile on the characteristics of hadronic mul-tiparticle production,” , International CosmicRay Conference, , 103 (2005).[21] Karl-Heinz Kampert and Michael Unger, “Measure-ments of the Cosmic Ray Composition with Air ShowerExperiments,” Astropart. Phys. , 660–678 (2012),arXiv:1201.0018 [astro-ph.HE].[22] J. Matthews, “A heitler model of extensive air showers,”Astropart. Phys. , 387–397 (2005).[23] Ralph Engel, Dieter Heck, and Tanguy Pierog, “Exten-sive air showers and hadronic interactions at high en-ergy,” Ann. Rev. Nucl. Part. Sci. , 467–489 (2011).[24] Hans Peter Dembinski, Balazs K´egl, Ioana C. Mari¸s,Markus Roth, and Darko Veberiˇc, “A likelihood methodto cross-calibrate air-shower detectors,” Astropart. Phys. , 44–51 (2016), arXiv:1503.09027 [astro-ph.IM].[25] Bruce Dawson et al. (Pierre Auger Collaboration), “TheEnergy Scale of the Pierre Auger Observatory,” ThePierre Auger Observatory: Contributions to the 36th In-ternational Cosmic Ray Conference (ICRC 2019) , PoS
ICRC2019 , 231 (2019).[26] Alexander Aab et al. (Pierre Auger Collaboration),“Measurement of the cosmic-ray energy spectrum above2 . × eV using the Pierre Auger Observatory,” Phys.Rev. D , 062005 (2020), arXiv:2008.06486 [astro-ph.HE].[27] (), See Supplemental Material at [URL will be insertedby publisher] for the comparison of the full model distri-bution for R µ in the individual energy bins.[28] The resolution systematics estimated in [25] are smallerby about a factor three, but were derived for differentquality cuts than the ones applied here.[29] (), See Supplemental Material at [URL will be insertedby publisher] for updated versions of select figures from the previous publication.[30] T. Pierog and K. Werner, “Epos model and ultra highenergy cosmic rays,” Nucl. Phys. Proc. Suppl. , 102–105 (2009).[31] T. Pierog, Iu. Karpenko, J. M. Katzy, E. Yatsenko, andK. Werner, “EPOS LHC: Test of collective hadronizationwith data measured at the CERN Large Hadron Col-lider,” Phys. Rev. C , 034906 (2015), arXiv:1306.0121[hep-ph].[32] Sergey Ostapchenko, “Monte Carlo treatment ofhadronic interactions in enhanced Pomeron scheme: I.QGSJET-II model,” Phys. Rev. D , 014018 (2011),arXiv:1010.1869 [hep-ph].[33] Eun-Joo Ahn, Ralph Engel, Thomas K. Gaisser, PaoloLipari, and Todor Stanev, “Cosmic ray interaction eventgenerator SIBYLL 2.1,” Phys. Rev. D , 094003 (2009).[34] Felix Riehn, Ralph Engel, Anatoli Fedynitch, Thomas K.Gaisser, and Todor Stanev, “Hadronic interaction modelSIBYLL 2.3d and extensive air showers,” Phys. Rev. D , 063002 (2020).[35] T. Bergmann et al. , “One-dimensional hybrid approachto extensive air shower simulation,” Astropart. Phys. ,420–432 (2007).[36] T. Pierog et al. , “First results of fast one-dimensionalhybrid simulation of eas using conex,” Nucl. Phys. Proc.Suppl. , 159–162 (2006).[37] Alfredo Ferrari et al. , “FLUKA: A multi-particle trans-port code,” CERN-2005-010, SLAC-R-773, INFN-TC-05-11 (2005).[38] T.T. B˘ohlen et al. , “The FLUKA Code: Developmentsand Challenges for High Energy and Medical Applica-tions,” Nuclear Data Sheets , 211–214 (2014).[39] Jose Bellido (Pierre Auger Collaboration), “Depth ofmaximum of air-shower profiles at the Pierre Auger Ob-servatory: Measurements above 10 . eV and Com-position Implications,” The Pierre Auger Observatory:Contributions to the 35th International Cosmic RayConference (ICRC 2017) , PoS
ICRC2017 , 506 (2018),[,40(2017)].[40] The mass fractions that are used for Sibyll here werederived for Sibyll 2.3. Given the small difference of5 g / cm in (cid:104) X max (cid:105) and σ ( X max ) between Sibyll 2.3 andSibyll 2.3d, the mass fractions for Sibyll 2.3d are expectedto be very similar.[41] (), See Supplemental Material at [URL will be inserted bypublisher] for the predictions by the individual models.[42] While the 68% contours of Epos and Sibyll cover a rangeof different mixtures, the QGSjet contour is elongatedalong the p-He line. The latter is an artifact originatingfrom the poor description of the measured X max distri-bution with the QGSjet model [54].[43] The energy fed to the electromagnetic shower in the firstinteraction is E − (cid:80) mj =1 E j , and it rapidly decreases fornext generations [55].[44] Roberto Aloisio, Denise Boncioli, Armando di Matteo,Piera L. Ghia, Aurelio F. Grillo, Sergio Petrera, andFrancesco Salamida, “Are Cosmic Rays still a valuableprobe of Lorentz Invariance Violations in the Auger era?” Proceedings, Vulcano Workshop 2014: Frontier Objectsin Astrophysics and Particle Physics: Vulcano, Italy,May 18-24, 2014 , Frascati Phys. Ser. , 274 (2014),arXiv:1408.5213 [astro-ph.HE].[45] J. Alvarez-Muniz, L. Cazon, R. Concei¸c˜ao, J. Dias de Deus, C. Pajares, and M. Pimenta, “Muon produc-tion and string percolation effects in cosmic rays at thehighest energies,” (2012), arXiv:1209.6474 [hep-ph].[46] Luis A. Anchordoqui, Haim Goldberg, and Thomas J.Weiler, “Strange fireball as an explanation of the muonexcess in Auger data,” Phys. Rev. D , 063005 (2017),arXiv:1612.07328 [hep-ph].[47] Glennys R. Farrar, “Particle Physics at Ultrahigh Ener-gies,” in Proceedings, 18th International Symposium onVery High Energy Cosmic Ray Interactions (ISVHECRI2014): Geneva, Switzerland, August 18-22, 2014 (2019)arXiv:1902.11271 [hep-ph].[48] Glennys R. Farrar and Jeffrey D. Allen, “A new phys-ical phenomenon in ultra-high energy collisions,”
Pro-ceedings, International Symposium on Future Directionsin UHECR Physics (UHECR2012): CERN, Geneva,Switzerland, February 13-16, 2012 , EPJ Web Conf. ,07007 (2013), arXiv:1307.2322 [hep-ph].[49] The number of generations is difficult to define as it de-pends on the details of the measurement and the showers,like the zenith angle, the muon energy threshold, distancefrom shower axis. Six is a reasonable minimum of genera-tions before the muons reach the Auger detectors [22, 56].[50] See Supplemental Material at [URL will be inserted bypublisher] for additional information on the average num- ber of muons and the fluctuations.[51] Alexander Aab et al. (Pierre Auger Collaboration), “ThePierre Auger Observatory Upgrade - Preliminary DesignReport,” (2016), arXiv:1604.03637 [astro-ph.IM].[52] Bjarni Pont et al. (Pierre Auger Collaboration), “A LargeRadio Detector at the Pierre AugerObservatory – Mea-suring the Properties of CosmicRays up to the HighestEnergies,” The Pierre Auger Observatory: Contributionsto the 36th International Cosmic Ray Conference (ICRC2019) , PoS
ICRC2019 , 395 (2019).[53] A. Aab et al. (Pierre Auger Collaboration), “Measure-ment of the radiation energy in the radio signal of ex-tensive air showers as a universal estimator of cosmic-rayenergy,” Phys. Rev. Lett. , 241101 (2016).[54] A. Aab et al. (Pierre Auger Collaboration), “Depth ofmaximum of air-shower profiles at the Pierre Auger Ob-servatory. II. Composition implications,” Phys. Rev. D , 122006 (2014), arXiv:1409.5083 [astro-ph.HE].[55] Lorenzo Cazon, “Probing High-Energy Hadronic Interac-tions with Extensive Air Showers,” PoS ICRC2019 , 005(2019), arXiv:1909.02962 [hep-ex].[56] Jaime Alvarez-Muniz, Ralph Engel, T.K. Gaisser, Jefer-son A. Ortiz, and Todor Stanev, “Hybrid simulations ofextensive air showers,” Phys. Rev. D66