Atmospheric muons measured with the KM3NeT detectors in comparison with updated numeric predictions

Eur. Phys. J. C (2024) 84:696 https://doi.org/10.1140/epjc/s10052-024-13018-8 Regular Article - Experimental Physics Atmospheric muons measured with the KM3NeT detectors in comparison with updated numeric predictions KM3NeT Collaboration,a S. Aiello1, A. Albert2,50, M. Alshamsi3,4, S. Alves Garre5, A. Ambrosone6,7,F.Ameli 8, M. Andre9, E. Androutsou10, M. Anguita11, L. Aphecetche12, M. Ardid13, S. Ardid13, A. Arsenic14, H. Atmani15, J. Aublin16, F. Badaracco17, L. Bailly-Salins18, Z. Bardaˇcová19,20, B. Baret16, A. Bariego-Quintana5, S. Basegmez du Pree21, Y. Becherini16, M. Bendahman15,16, F. Benfenati22,23, M. Benhassi6,24, D. M. Benoit25, E. Berbee21, V. Bertin3, S. Biagi26, M. Boettcher27, D. Bonanno26, J. Boumaaza15, M. Bouta28, M. Bouwhuis21, C. Bozza6,29, R. M. Bozza6,7, H. Brânza¸s30, F. Bretaudeau12, M. Breuhaus3, R. Bruijn21,31, J. Brunner3, R. Bruno1,E.Buis 21,32, R. Buompane6,24,J.Busto 3,B.Caiffi 17,D.Calvo 5, S. Campion8,33, A. Capone8,33, F. Carenini22,23, V. Carretero21,31, T. Cartraud16, P. Castaldi22,34, V. Cecchini5, S. Celli8,33, L. Cerisy3, M. Chabab35, M. Chadolias14, A. Chen36, S. Cherubini26,37, T. Chiarusi22, M. Circella38, R. Cocimano26, J. A. B. Coelho16, A. Coleiro16, A. Condorelli6,7, R. Coniglione26, P. Coyle3, A. Creusot16, G. Cuttone26, R. Dallier12, Y. Darras14, A. De Benedittis6, B. De Martino3, V. Decoene12,R.DelBurgo 6, I. Del Rosso22,23,L.S.DiMauro 26,I.DiPalma 8,33,A.F.Díaz 11,C.Diaz 11, D. Diego-Tortosa26, C. Distefano26,A.Domi 14, C. Donzaud16, D. Dornic3, M. Dörr51, E. Drakopoulou10, D. Drouhin2,50, J.-G. Ducoin3, R. Dvornický20, T. Eberl14, E. Eckerová19,20, A. Eddymaoui15, T. van Eeden21, M. Eff16,D.vanEijk 21, I. El Bojaddaini28, S. El Hedri16, A. Enzenhöfer3, G. Ferrara26, M. D. Filipovi´c39, F. Filippini22,23, D. Franciotti26, L. A. Fusco6,29, S. Gagliardini8,33,T.Gal 14, J. García Méndez13, A. Garcia Soto5, C. Gatius Oliver21, N. Geißelbrecht14, H. Ghaddari28, L. Gialanella6,24,B.K.Gibson 25, E. Giorgio26, I. Goos16, P. Goswami16, S. R. Gozzini5, R. Gracia14, K. Graf14, C. Guidi17,40, B. Guillon18, M. Gutiérrez41, C. Haack14, H. van Haren42, A. Heijboer21, A. Hekalo51, L. Hennig14, J. J. Hernández-Rey5, W. Idrissi Ibnsalih6, G. Illuminati22,23,D.Joly 3, M. de Jong21,43, P. de Jong21,31,B.J.Jung 21, P. Kalaczy´nski52, O. Kalekin14, U. F. Katz14, A. Khatun20, G. Kistauri44,45, C. Kopper14, A. Kouchner16,46, V. Kueviakoe21, V. Kulikovskiy17,b, R. Kvatadze45, M. Labalme18, R. Lahmann14, G. Larosa26, C. Lastoria3, A. Lazo5,S.LeStum 3, G. Lehaut18, E. Leonora1, N. Lessing5,G.Levi 22,2, F. Longhitano1, F. Magnani3, J. Majumdar21, L. Malerba17, F. Mamedov19, J. Ma´nczak5, A. Manfreda6, M. Marconi17,40, A. Margiotta22,23, A. Marinelli6,7, C. Markou10, L. Martin12, F. Marzaioli6,24, M. Mastrodicasa8,33, S. Mastroianni6, S. Miccichè26, G. Miele6,7, P. Migliozzi6, E. Migneco26, M. L. Mitsou6, C. M. Mollo6, L. Morales-Gallegos6,24, G. Moretti8,33, A. Moussa28, I. Mozun Mateo18, R. Muller22, M. R. Musone6,24, M. Musumeci26,S.Navas 41, A. Nayerhoda38, C. A. Nicolau8,B.Nkosi 36, B. Ó. Fearraigh17, V. Oliviero6,7, A. Orlando26, E. Oukacha16, D. Paesani26, J. Palacios González5, G. Papalashvili38,44, V. Parisi17,40, E. J. Pastor Gomez5,A.M.P˘aun30,G.E.P˘av˘ala¸s30, I. Pelegris10, S. Peña Martínez16, M. Perrin-Terrin3, J. Perronnel18, V. Pestel18,R.Pestes 16, P. Piattelli26, C. Poirè6,29, V. Popa30, T. Pradier2, J. Prado5, S. Pulvirenti26, C. A. Quiroz-Rangel13, U. Rahaman5, N. Randazzo1, S. Razzaque47,I.C.Rea 6, D. Real5, G. Riccobene26, J. Robinson27, A. Romanov17,18,40,c, A. Šaina5, F. Salesa Greus5, D. F. E. Samtleben21,43, A. Sánchez Losa5, S. Sanfilippo26, M. Sanguineti17,40, C. Santonastaso6,24, D. Santonocito26, P. Sapienza26, J. Schnabel14, J. Schumann14, H. M. Schutte27, J. Seneca21, N. Sennan28, B. Setter14, I. Sgura38, R. Shanidze44, A. Sharma16, Y. Shitov19,F.Šimkovic 20, A. Simonelli6, A. Sinopoulou1,M.V.Smirnov 14, B. Spisso6, M. Spurio22,23, D. Stavropoulos10, I. Štekl19, M. Taiuti17,40, Y. Tayalati15, H. Thiersen27, I. Tosta e Melo1,37, E. Tragia10, B. Trocmé16, V. Tsourapis10, A. Tudorache8,33, E. Tzamariudaki10, A. Vacheret18, A. Valer Melchor21, V. Valsecchi26, V. Van Elewyck16,46, G. Vannoye3, G. Vasileiadis48, F. Vazquez de Sola21, A. Veutro8,33, S. Viola26, D. Vivolo6,24, J. Wilms49,E.deWolf 21,31, H. Yepes-Ramirez13,I.Yvon 16, G. Zarpapis10, S. Zavatarelli17, A. Zegarelli8,33,D.Zito 26, J. D. Zornoza5, J. Zúñiga5, N. Zywucka27 1INFN, Sezione di Catania, (INFN-CT), Via Santa Sofia 64, 95123 Catania, Italy 2Université de Strasbourg, CNRS, IPHC UMR 7178, 67000 Strasbourg, France 3Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France 4University of Sharjah, Sharjah Academy for Astronomy, Space Sciences, and Technology, University Campus, POB 27272, Sharjah, 123 696 Page 2 of 19 Eur. Phys. J. C (2024) 84 :696 United Arab Emirates 5IFIC-Instituto de Física Corpuscular (CSIC-Universitat de València), c/Catedrático José Beltrán, 2, 46980 Paterna, Valencia, Spain 6INFN, Sezione di Napoli, Complesso Universitario di Monte S. Angelo, Via Cintia ed. G, 80126 Naples, Italy 7Dip. Scienze Fisiche “E. Pancini”, Complesso Universitario di Monte S. Angelo, Università di Napoli “Federico II”, Via Cintia ed. G, 80126 Naples, Italy 8INFN, Sezione di Roma, Piazzale Aldo Moro 2, 00185 Rome, Italy 9Laboratori d’Aplicacions Bioacústiques, Centre Tecnològic de Vilanova i la Geltrú, Universitat Politècnica de Catalunya, Avda. Rambla Exposició, s/n, 08800 Vilanova i la Geltrú, Spain 10 NCSR Demokritos, Institute of Nuclear and Particle Physics, Ag. Paraskevi Attikis, 15310 Athens, Greece 11 Dept. of Computer Architecture and Technology/CITIC, University of Granada, 18071 Granada, Spain 12 Subatech, IMT Atlantique, IN2P3-CNRS, Nantes Université, 4 rue Alfred Kastler-La Chantrerie, BP 20722 44307, Nantes, France 13 Instituto de Investigación para la Gestión Integrada de las Zonas Costeras, Universitat Politècnica de València, C/ Paranimf, 1, 46730 Gandia, Spain 14 Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Straße 2, 91058 Erlangen, Germany 15 University Mohammed V in Rabat, Faculty of Sciences, 4 av. Ibn Battouta, B.P. 1014, R.P. 10000 Rabat, Morocco 16 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France 17 INFN, Sezione di Genova, Via Dodecaneso 33, 16146 Genoa, Italy 18 LPC CAEN, Normandie Univ, ENSICAEN, UNICAEN, CNRS/IN2P3, 6 boulevard Maréchal Juin, 14050 Caen, France 19 Institute of Experimental and Applied Physics, Czech Technical University in Prague, Husova 240/5, Prague 110 00, Czech Republic 20 Department of Nuclear Physics and Biophysics, Comenius University in Bratislava, Mlynska dolina F1, 842 48 Bratislava, Slovak Republic 21 Nikhef, National Institute for Subatomic Physics, PO Box 41882, 1009 DB Amsterdam, Netherlands 22 INFN, Sezione di Bologna, v.le C. Berti-Pichat, 6/2, 40127 Bologna, Italy 23 Dipartimento di Fisica e Astronomia, Università di Bologna, v.le C. Berti-Pichat, 6/2, 40127 Bologna, Italy 24 Dipartimento di Matematica e Fisica, Università degli Studi della Campania “Luigi Vanvitelli”, viale Lincoln 5, 81100 Caserta, Italy 25 E. A. Milne Centre for Astrophysics, University of Hull, Hull HU6 7RX, UK 26 INFN, Laboratori Nazionali del Sud, (LNS) Via S. Sofia 62, 95123 Catania, Italy 27 Centre for Space Research, North-West University, Private Bag X6001, Potchefstroom 2520, South Africa 28 Faculty of Sciences, BV Mohammed VI, University Mohammed I, B.P. 717, R.P. 60000 Oujda, Morocco 29 Dipartimento di Fisica, Università di Salerno e INFN Gruppo Collegato di Salerno, Via Giovanni Paolo II 132, 84084 Fisciano, Italy 30 ISS, Atomistilor 409, 077125 M˘agurele, Romania 31 Institute of Physics/IHEF, University of Amsterdam, PO Box 94216, 1090 GE Amsterdam, Netherlands 32 TNO, Technical Sciences, PO Box 155, 2600 AD Delft, Netherlands 33 Dipartimento di Fisica, Università La Sapienza, Piazzale Aldo Moro 2, 00185 Rome, Italy 34 Dipartimento di Ingegneria dell’Energia Elettrica e dell’Informazione ”Guglielmo Marconi”, Università di Bologna, Via dell’Università 50, Cesena 47521, Italy 35 Physics Department, Faculty of Science Semlalia, Cadi Ayyad University, Av. My Abdellah, P.O.B. 2390, 40000 Marrakech, Morocco 36 School of Physics, University of the Witwatersrand, Wits, Private Bag 3, Johannesburg 2050, South Africa 37 Dipartimento di Fisica e Astronomia “Ettore Majorana”, Università di Catania (INFN-CT), Via Santa Sofia 64, 95123 Catania, Italy 38 INFN, Sezione di Bari, via Orabona, 4, 70125 Bari, Italy 39 School of Computing, Engineering and Mathematics, Western Sydney University, Locked Bag 1797, Penrith, NSW 2751, Australia 40 Università di Genova, Via Dodecaneso 33, 16146 Genoa, Italy 41 Dpto. de Física Teórica y del Cosmos & C.A.F.P.E., University of Granada, 18071 Granada, Spain 42 NIOZ (Royal Netherlands Institute for Sea Research), PO Box 59, Den Burg, 1790 AB Texel, The Netherlands 43 Leiden Institute of Physics, Leiden University, PO Box 9504, 2300 RA Leiden, Netherlands 44 Department of Physics, Tbilisi State University, 3, Chavchavadze Ave., 0179 Tbilisi, Georgia 45 Institute of Physics, The University of Georgia, Kostava str. 77, 0171 Tbilisi, Georgia 46 Institut Universitaire de France, 1 rue Descartes, 75005 Paris, France 47 Department Physics, University of Johannesburg, PO Box 524, Auckland Park 2006, South Africa 48 Laboratoire Univers et Particules de Montpellier, Place Eugène Bataillon-CC 72, 34095 Montpellier Cédex 05, France 49 Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Remeis Sternwarte, Sternwartstraße 7, 96049 Bamberg, Germany 50 Université de Haute Alsace, rue des Frères Lumière, 68093 Mulhouse Cedex, France 51 Fakultät für Physik und Astronomie, Institut für Theoretische Physik und Astrophysik, Lehrstuhl für Astronomie, Julius-Maximilians-Universität Würzburg, Emil-Fischer-Straße 31, 97074 Würzburg, Germany 52 AstroCeNT, Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Rektorska 4, 00-614 Warsaw, Poland Received: 18 March 2024 / Accepted: 11 June 2024 / Published online: 15 July 2024 © The Author(s) 2024 123 Eur. Phys. J. C (2024) 84 :696 Page 3 of 19 696 Abstract The measurement of the flux of muons produced in cosmic ray air showers is essential for the study of primary cosmic rays. Such measurements are important in extensive air shower detectors to assess the energy spectrum and the chemical composition of the cosmic ray flux, complementary to the information provided by fluorescence detectors. Detailed simulations of the cosmic ray air showers are carried out, using codes such as CORSIKA, to estimate the muon flux at sea level. These simulations are based on the choice of hadronic interaction models, for which improvements have been implemented in the post-LHC era. In this work, a deficit in simulations that use state-of-the-art QCD models with respect to the measurement deep underwater with the KM3NeT neutrino detectors is reported. The KM3NeT/ARCA and KM3NeT/ORCA neutrino telescopes are sensitive to TeV muons originating mostly from primary cosmic rays with energies around 10 TeV. The predictions of state-of-the-art QCD models show that the deficit with respect to the data is constant in zenith angle; no dependency on the water overburden is observed. The observed deficit at a depth of several kilometres is compatible with the deficit seen in the comparison of the simulations and measurements at sea level. 1 Introduction Primary Cosmic Rays (CRs) are ionized nuclei, mainly protons,thatwerediscoveredbyVictorHessandothersinanumberofballoonflights[1].Despiteextensivestudiessincethen, the origin, propagation, and interaction of high-energy CRs in the atmosphere are not yet fully understood [2]. This work focuses on the latter point demonstrating a deficit of deeply penetratingmuonsinthesimulationsofCRinteractionsinthe atmosphere using the up-to-date hadronic interaction models. Although the KM3NeT data are in better agreement with older parametric simulations [3] (Figure 2), those are based on underground and underwater muon measurements and the hadronic interaction models used are inconsistent with the LHC results. The measurement of the number of muons in CR events with the Pierre Auger Observatory has revealed a deficit of GeV muons in CR simulations with respect to the data [4,5]. This result triggered similar analyses by several experiments, such as the EAS-MSU Array [6], the IceCube Neutrino Observatory [7], the KASCADE-Grande experiment [8], the NEVOD-DECOR detector [9], the SUGAR array [10], the Telescope Array [11], and the Yakutsk array [12]. ae-mail: [email protected] be-mail: vladimir.kuliko[email protected] (corresponding author) ce-mail: [email protected] (corresponding author) For several of these experiments, the muon deficit becomes more and more pronounced as the energy increases, starting from 20 to 40 PeV up to the highest energies of tens of EeV. The issue connected with this deficit is known as the Muon Puzzle [13]. In this work, the zenith distribution of the rate of highenergy muons arising from Extensive Air Showers (EAS) is probed at a depth of several kilometres with the KM3NeT underwatertelescopes.Theminimal muon energyatsea level required to reach the KM3NeT detectors is around 500 GeV [14] with the majority of muons having energies in the TeV range as shown below in Sect.3. A verification that the production of TeV muons is described properly serves also as a check that the neutrino and gamma ray fluxes from cosmic sources are correctly predicted.Sinceneutrinosandmuonsaremostlyproducedinthe weak decays of mesons, one may expect that the increase of muons is due to an increase of mesons decaying predominantly in muons, and compensated by a decrease in mesons decaying electromagnetically with a production of gamma rays (see, for example, [15,16]). The increased neutrino production and decreased gamma ray production means that the gamma ray to neutrino conversion predictions for the cosmic sources used in studies of the neutrino telescope collaborations could gain a boost for the neutrino fluxes [17–19]. Several other large-volume neutrino telescopes reported results on atmospheric muon measurements, namely AMANDA [20], ANTARES [21], Baikal-GVD [22], IceCube [23], NEMO [24], and NESTOR [25]. The quoted resultswerecomparedwiththepre-LHChadronicinteraction models in contrast to this work where the post-LHC models are used. This paper is organised as follows. The description of the KM3NeT detectors is given in Sect.2. Details on the models used in the simulation are provided in Sect.3. The motivation and procedure of the tuning of a fast Monte Carlo (MC) generator based on the detailed MC simulations are described in Sect.4. The muon reconstruction performance and the estimation of systematic uncertainties are given in Sect.5and Sect.6, respectively. The main results of the comparison of KM3NeT data with MC simulations are presented and discussed in Sect.7. Finally, the conclusions are summarised in Sect.8. 2 The KM3NeT neutrino telescopes TheKM3NeTresearchinfrastructurecomprisestwoneutrino telescopes at the bottom of the Mediterranean Sea [26]. The detection technology is the same for both detectors. In the KM3NeT telescopes, Cherenkov photons induced by highlyrelativistic chargedparticles are detected byanarray of threeinch PhotoMultiplier Tubes (PMTs). The PMTs are placed in 123 696 Page 4 of 19 Eur. Phys. J. C (2024) 84 :696 pressure-resistant 17-inch diameter glass spheres and there are thirty one PMTs in each sphere. This covers most of the light arrival directions which is suitable for the downwardgoing atmospheric muons detection. Such sphere, together with the associated readout electronics inside it [27], comprises the Digital Optical Module (DOM) [28]. The KM3NeT DOMs are attached to a pair of vertical ropes, forming a Detection Unit (DU) consisting of eighteen DOMs. Each DU is secured to the sea floor with an anchor and has a buoy at the top to keep the DU in an almost vertical position. The detector geometry varies slowly with deepwater currents. The position and orientation of each DOM are continuously measured with an acoustic positioning system, and with tilt and compass sensors [29]. The dynamic positioning calibration of the detectors can be further corroborated and enhanced using muons passing through the detector, by means of a likelihood maximisation procedure [30]. Whenever the signal pulse read-out of a PMT exceeds a predefined threshold voltage, the signal is recorded. In particular, the start time of the signal and the pulse duration are saved. Together with the PMT identification number these data form a “hit”. A set of causally-connected hits triggers an “event”. The data acquisition is based on the all-data-toshore concept [26]. The event triggering and reconstruction are performed onshore. The data are taken continuously in consecutive runs of a few hours. Neutrino telescopes can detect different event topologies, comingfrom all-flavour neutrino interactionsand from atmospheric muons. Track-like events correspond to muons crossing the instrumented volume of the detector. This topology is mostly caused by atmospheric muons or by muons from charged-current neutrino interactions. Shower-like (or cascade-like) events are caused by interactions of neutrinos where only electromagnetic and hadronic showers are produced. In this work, events are reconstructed under the assumption of a track-like topology [31]. The overwhelming majorityofsuchreconstructedeventsisduetothedownwardgoing atmospheric muons. The KM3NeT/ARCA telescope is under construction offshore the coast of Sicily, Italy, at a depth of about 3.5 km. In its final configuration, the KM3NeT/ARCA detector will consist of two “building blocks”, each comprising 115 DUs. The horizontal spacing between the DUs is about 90m and the vertical distance between the DOMs on each DU is 36m, optimised for the detection of high-energy cosmic neutrino signals. The KM3NeT/ORCA detector is located about 40km offshore the coast of Toulon, France, at about 2.5 km depth. KM3NeT/ORCA will comprise one building block of 115 DUs. The distance between the DUs is around 20m and the DOMs on each DU are on average 9m apart, in order to measure neutrino oscillations of GeV atmospheric neutrinos. The overburden difference between the KM3NeT detectors is almost one kilometre. This, together with different detector energy thresholds and geometries, allows for measurements of deeply-penetrating atmospheric muons in complementary energy ranges. The analysis presented in this work is performed with data taken by both the KM3NeT/ARCA and KM3NeT/ORCA telescopes in a configuration of six DUs at each site. In the following, these configurations are denoted as ARCA6 and ORCA6 for the two detectors, correspondingly. 3 Simulation of high-energy muons in KM3NeT The simulation of muons for the KM3NeT telescopes is divided into four steps: 1. Simulation of the interactions of the primary CRs with the airnucleiandthesubsequentpropagation,interaction,and decay in the atmosphere of the secondary particles. This is done with a detailed MC simulation software. By detailed simulations, it is implied that all the relevant interaction, decay, and energy loss processes are taken into account and treated stochastically using step-by-step Monte Carlo sampling methods. 2. Propagation of muons in water from sea level down to the KM3NeT detectors. The propagation routine uses MC methods for proper treatment of the stochastic muon energy losses and the muon direction changes. 3. Generation of the Cherenkov light from muons and its detection by the PMTs. Fast sampling of the PMT photoelectrons for nominal quantum efficiency is done from pre-generated tables. Conversion of the simulated photoelectrons to hits is performed using the measured efficiencies, PMT timing characteristics, and the optical background. 4. Reconstruction of muon tracks. 3.1 Simulations of primary cosmic-ray interactions For the first step, CORSIKA version 7.7410 [32] has been used for this work. Since the minimal muon energy required toreachthe depthofthe KM3NeTdetectorsisabout500GeV at sea level [14], the lower limit on the primary CR energy is conservatively set to 1 TeV per nucleon in the simulation. For each nucleus, the energy range in which simulations are carried out is divided into 3 sub-ranges (named “TeV low”, “TeV high”, and “PeV”) to obtain sufficient statistics for all energies of interest. A summary of the production energy ranges and the number of generated CR showers in each production is reported in Table 1. Five nuclei are used as primaries in the simulation: hydrogen (1 1H), helium (4 2He), carbon (12 6C), oxygen (16 8O), and 123 Eur. Phys. J. C (2024) 84 :696 Page 5 of 19 696 Table 1 Five nuclei and three energy ranges have been used for the simulation using the CORSIKA code. The Table reports for each nucleus and for the three energy bands the interval of energy and, in brackets, the number of simulated events Nucleus Energy range per nucleus [TeV] (Number of generated showers) TeV low TeV high PeV 1 1H 1–6(3×109)6–1.1×103(3 ×109)1.1×103–9×104(4 ×107) 4 2He 4 – 10 (2 ×109)10–1.1×103(2 ×109)1.1×103–9×104(4 ×107) 12 6C12–30(1×109)30–1.1×103(7.5 ×108)1.1×103–9×104(2.5 ×107) 16 8O16–30(1×109)30–1.1×103(7.5 ×108)1.1×103–9×104(2.5 ×107) 56 26Fe 56 – 100 (5 ×108) 100–1.1×103(2 ×108)1.1×103–9×104(2 ×107) iron (56 26Fe). These simulation samples are used together to describe the expected muon flux at sea level. Each sample must be weighted in order to reproduce the input primary spectrum and its composition. Other primaries are taken into account by increasing the flux weights of C, O, and Fe, according to the flux of nuclei that have not been simulated. The Global Spline Fit (GSF) is chosen as the model for the CR flux and mass-composition [33]. This model is a data-driven parameterisation of the CR flux, obtained using only the assumption that the flux is smoothly varying with energy. Since the GSF, in contrast with other mass composition models such as H3a [34] and GST [35], does not rely on the assumption that the primary CR components should be described by power-law energy spectra with rigiditydependent cutoffs, it describes data without assumptions on the physical origin of the measured CRs. Additionally, the GSF incorporates the data-driven uncertainties on the total flux and on the flux of individual nuclear species, which are used in this work for the systematic uncertainties studies described in Sect.6. The UrQMD 1.3 model [37] is used to simulate the elastic andinelastic interactionsof hadronsbelow80 GeVinair.The simulation of high-energy hadronic interactions is performed with the post-LHC model Sibyll 2.3d [38]. The flux of particles produced in the CR air shower depends on the density profile of the atmosphere. NRLMSIS 2.0 [36] is used to obtain the atmosphere density profile. Since the atmospheric conditions vary depending on time and location, a fit is performed on the model predictions averaged over a 3 year period (2019–2021) and over the KM3NeT/ORCA and KM3NeT/ARCA sites, 42◦48N 06◦02E and 36◦16N16 ◦06E, respectively. In CORSIKA the atmosphere is divided into five layers with the fifth layer at the top of the atmosphere. The atmospheric overburden as a function of the altitude above sea level, T(h)(g/cm2), is parameterised in each layer [32]. For the first four layers, T(h)is parameterised as an exponential function: T(h)=ai+bi·e−h ci,1≤i≤4,(1) while in the fifth layer the T(h)dependence is linear: T(h)=a5−b5·h c5.(2) Here, iis the layer number, and ai,bi, and ciare the fitting parameters. The values of fitted parameters are reported in Table 2. 3.2 Propagation of muons in seawater To obtain muon kinematic properties in seawater at a cylindricalsurfacesurroundingtheKM3NeTdetector,the can,the gSeaGen code [39] is used. The propagation of atmospheric muons resulting from CORSIKA is added to the gSeaGen code for this work since originally gSeaGen was developed as a GENIE-based [40] application for the neutrino simulation. The muon propagation routine used in gSeaGen is the PROPOSAL software [41], specifically version 6 of the code. The can height and radius correspond to those of the instrumented volume enlarged by 4 times the maximum of the absorption lengths in seawater, i.e., 70m ×4, to account for the light emitted along the path of muons passing outside the can, but close enough to reach the instrumented volume. The following muon energy loss processes are taken into account in the PROPOSAL calculations: ionisation (Groom, Mokhov, Striganov parameterisation [42] with the density effect correction calculated using Sternheimer parameterisation [43]), bremsstrahlung (Kelner et al. [44]), e+/e−pair production(KokoulinandPetrukhin[45]),and photo-nuclear interactions (Abramowicz and Levy [46]). The LandauPomeranchuk-Migdal effect [47] is also accounted for in PROPOSAL, although the effect becomes noticeable only at muon energies greater than 10 PeV; it is negligible for atmospheric muons detected in KM3NeT. The uncertainty on the underwater flux of muons caused by the inaccuracies in the models of the energy loss processes is about a few per cent at most [48]. The seawater composition and mass density used in the simulations correspond to the “ANTARESWater” model in PROPOSAL. The composition is taken from [49] and corrected for the salinity at the ANTARES site (38.480 /00 [50]). The ANTARES location was near the KM3NeT/ORCA 123 696 Page 6 of 19 Eur. Phys. J. C (2024) 84 :696 Table 2 Values of the parameters used in Eqs.(1)and(2) resulting from the T(h)fit Layer Height (km) Fitted aivalue (g/cm2) Fitted bivalue (g/cm2) Fitted civalue (cm) 1 0–17.5 −58.45 10.71 ×1028.65 ×105 2 17.5–45 0.50 13.81 ×1026.23 ×105 3 45–73 −0.015 4.50 ×1027.90 ×105 4 73–101.3 −4.29 ×10−55.71 ×1035.99 ×105 5 101.3–125 3.0 ×10−338.86 1.41 ×1011 Table 3 Relative content of the seawater elements used in the “ANTARESWater” model and in the updated model obtained in this work. The content of hydrogen atoms is predefined to be equal to 2, the content of other elements is calculated relative to hydrogen Model H O Na K Mg Ca Cl S “ANTARESWater” 2.0 1.0088 0.00943 0.000209 0.001087 0.000209 0.01106 0.00582 Updated 2.0 1.0024 0.00962 0.000209 0.001083 0.000211 0.01119 0.000579 detector; the salinity at the KM3NeT/ARCA site differs from theANTARESone atalevelbelow1%[51]. With thiscorrection, the seawater content of all elements except for hydrogen and oxygen is multiplied by a factor corresponding to the relative difference in salinity (38.48/35); the correction for the oxygen coming from the sea salt sulphate was neglected. In this work, the updated seawater composition [52] (Table 4), corrected to account also for oxygen in sulphates, is tested. First, the abundance of each sea salt molecule is re-scaled accounting for higher salinity at the ANTARES site with respecttothe salinityof the“ReferenceSeawater”,i.e., 350 /00. Then, molar masses of seawater molecules are calculated as the ratio between the relative content of molecules constituting the sea salt and the atomic weights of elements as listed in PROPOSAL, which also accounts for different isotope abundances.Afterthat, the hydrogen atomiccontentis normalised such that there are two atoms of hydrogen in the “solution” which is used in PROPOSAL. The contribution of other elements is then re-scaled using the proportion obtained in the hydrogen normalisation. The updated seawater composition and the composition from the “ANTARESWater” model are given in Table 3. A test is performed in order to evaluate the influence of the new water composition on the muon energy losses. The differences in the energy losses between the two composition models for TeV muons travelling 2 and 3km in water are found at a level below 0.5%. The original “ANTARESWater” model is used for the following in this work. 3.3 Detector response simulation The light generation, detector response simulation, and track reconstruction are performed with internal KM3NeT software. The detector response simulation takes into account the calibrated position and orientation of the PMTs, values of PMT quantum efficiencies and the environmental optical background. This background includes the contribution of photons induced in the water by the products of the 40K β-decays, bioluminescent light [53], and other radioactive decays of nuclei in the DOM glass and water. 3.4 Characteristics of the reconstructed events The same reconstruction algorithms are applied to both data and simulations assuming the track-like topology of events. The events produced by the optical background are filtered by applying quality cuts on the output of the track reconstruction algorithm. Cuts are applied on the number of hits used by the algorithm, nhits >50, and on the likelihood of the reconstruction, >20. Hereafter, those events that passed such quality cuts are referred to as “reconstructed events”. The rate of reconstructed events in the simulations as a function of the nucleus energy of primary CRs is shown in Fig.1. The highlighted area indicates the 90% fraction of the total number of events counting from the maximum of the distribution, i.e., the highest density interval. This energy range spans from 3 to 320 TeV for ORCA6 and from 4 TeV to 1 PeV for ARCA6; the median values of the distributions are 26 TeV and 60 TeV, respectively. Most of the events originate from proton and helium nuclei. Muons from the same CR interaction arrive in packed bundles at the KM3NeT detectors and they are reconstructed as a single track. The expected rate of muons at sea level from the bundles that are reconstructed with the ORCA6 and ARCA6 detectors as a function of the individual muon energy is shown in Fig.2for different ranges of the cosines of the zenith angle, cos θ, where cos θ=1 indicates vertically down-going muons. Each muon from the reconstructed bundles that reached the can is used to fill the distributions. 123 Eur. Phys. J. C (2024) 84 :696 Page 7 of 19 696 Fig. 1 The rate of simulated atmospheric muon events reconstructed with the ORCA6 (left) and ARCA6 (right) detectors as a function of the primary CR energy. The unshaded area corresponds to the 90% fraction of events (the highest density interval). Statistical uncertainties are shown as vertical error bars More inclined muons have higher energies at sea level since a larger amount of the water overburden is travelled. The pseudorapidity variable characterises the direction of the outgoing muon with respect to the primary CR. It is defined as η=−ln[tan(α/2)], where αis the angle between the primary particle and the secondary muon at sea level. The pseudorapidity distributions of reconstructed muons reaching the ORCA6 and ARCA6 detectors are plotted in Fig.3. The distributions are filled using each muon that reaches the can. The peak of the distributions is located at η≃9. Hence, muons detected by the KM3NeT telescopes originate from hadronic interactions in the very forward region of pseudorapidities. This region is not fully covered by accelerator experiments [13]. 4 Tuning of the MUPAGE parameters on the CORSIKA simulation The CORSIKA software provides detailed simulations of EAS. Given the steeply falling primary cosmic ray spectrum, high-energy events are much rarer than low-energy events, and the properties of air showers are highly stochastic. Thus, in order to get representative statistics of atmospheric muons deep underwater, the demand for CPU time for simulations is high. The KM3NeT simulations are based on the Run-byRun (RbR) approach similar to the ANTARES one [54], i.e., simulations are subdivided into batches corresponding to the actual data-taking periods to reproduce the time variability of the detector conditions. In particular, environmental noise is sampled from the data, the calibrated position and orientation of PMTs and the values of the PMT quantum efficiencies are considered individually for each run in the simulation. A separate muon sample is generated for each run to prevent bias associated with using the same sample. Performing RbR simulations with CORSIKA is not presently feasible due to the extensive CPU time required. Currently,thesimulationofatmosphericmuonsinKM3NeT is based on the fast MC generator MUPAGE [55,56]. It generates the muon bundle kinematic properties for a certain water depth and zenith angle at a plane perpendicular to the bundle axis. This generation is based on parametric formulas describing the flux of single and multiple muons in the bundle, the differential energy spectrum, and the distance from the bundle axis. The default MUPAGE parameters were obtained from a detailed EAS simulation performed with the HEMAS package [57] and then fitting the results to the MACRO measurements [58]. A framework has been developed to tune the MUPAGE parameters using the CORSIKA simulation software and the most recent available models to describe high-energy hadronic interactions, Sibyll 2.3d [38], and the CR flux, GSF [33]. The goal of this tuning is to combine the advantage of a quick parameterised simulation with the features coming from a detailed MC simulation. The MUPAGE code and the parametric formulas themselves are kept unmodified. The numerical values of some parameters have been changed according to the method described below. In this tuning procedure, atmospheric muons at sea level from CORSIKA simulations are propagated down through thewatertoaplaneperpendiculartotheCRshoweraxisusing thePROPOSALsoftware[41].The propagationis performed for seven different depths, from 2km down to 3.5 km in steps of 250m. Muon bundle kinematic properties are fitted at each depth following the approach described in the MUPAGE paper [55]. Five fits are performed to obtain the new values of the MUPAGE parameters that aim to describe the following muon bundle kinematic properties: 123 696 Page 8 of 19 Eur. Phys. J. C (2024) 84 :696 Fig. 2 Sea level rate of generated muons from the reconstructed bundles as a function of the muon energy. Different colours indicate different zenith angle ranges; vertical muons have cosθ=1. Statistical uncertainties are shown as vertical error bars Fig. 3 Pseudorapidity distribution of muons reaching the ORCA6 (left plot) and ARCA6 (right plot) cans. Only reconstructed events are used in the distribution. The unshaded area shows the pseudorapidity range that includes the 90% fraction of events (the highest density interval). Statistical uncertainties are shown as vertical error bars – flux of single muons in the bundles as a function of the zenith angle, – flux of multiple muons in the bundles as a function of the bundle multiplicity, – normalised energy spectrum of single muons in the bundles, – lateral spread of multiple muons in the bundles, – normalised energy spectrum of multiple muons in the bundles. The new values of the MUPAGE parameters obtained using CORSIKA simulations are given in Tables 4and 5. Comparison between the MUPAGE tuned on CORSIKA, original MUPAGE, and CORSIKA for the five aforementioned kinematic properties of muon bundles at the can are plotted in Fig.4. A comparison of the tuned and unmodified (nominal) MUPAGE withCORSIKAforthesinglemuon flux as a funcTable 4 Best fit values of the parameters [55] that describe the zenith dependence of the flux of single muons in the bundles and the number of muons in the bundles obtained in this work Parameter name Best fit value Parameter name Best fit value K0a7.98 ×10−3ν0a−6.48 ×10−2 K0b−1.896 ν0b0.433 K1a−0.606 ν0c2.475 K1a−0.110 ν1a4.27 ×10−2 ν1b0.476 tion of the cosine of the simulated zenith angle is given in Fig.4a. The MUPAGE functions with the tuned parameters well describe the CORSIKA distributions for muons with a direction ranging from vertical to very inclined. The statistical uncertainties are underestimated for the very inclined CORSIKA muons where the proton contribution is missing due to the lack of statistics in the “TeV high” production (see 123 Eur. Phys. J. C (2024) 84 :696 Page 9 of 19 696 Table 5 Best fit values of the parameters [55] that describe the energy dependence of single muon bundles and the muon lateral spread obtained in this work Parameter name Best fit value Parameter name Best fit value γ0−0.343 ρ0a−2.126 γ13.991 ρ0b27.23 0a4.41 ×10−4ρ1−1.025 0b1.288 θ010.0 1a−2.36 ×10−2α0a−1.104 1b0.765 α0b7.493 α1a7.94 ×10−2 α1b9.86 ×10−2 Table 1for the production label and Fig.1for the number of simulated proton showers). The muon bundle fluxes are compared in Fig.4b. A good agreement between the tuned MUPAGE and CORSIKA is present up to a multiplicity of five. The difference for larger multiplicities can be ascribed to the radial distribution of muons in bundles: muons in the bundles with multiplicity four and above are sampled from the same radial distribution in the MUPAGE code, which is dominated by bundles consisting of four muons. However, in the CORSIKA simulations the radial distribution becomes more compact as the multiplicity grows, hence, more muons are crossing the can with respect to MUPAGE at high multiplicities. A possible fix for the issue is the change of the MUPAGE formulas which may be performed in future works. Here, the discrepancy in the bundle multiplicity does not affect the results as shown in Fig.6since the inclusive muon rate underwater is dominated by bundles with low multiplicity. With inclusive atmospheric muon flux all muons originating from the whole CR energy spectrum are indicated, in contrast to the muons measured in EAS detectors from Ultra-High Energy (UHE) CR interactions only. The energy spectra of single muons are compared in Fig.4c. Both the nominal and the tuned MUPAGE follow the spectrum obtained with CORSIKA. The fit of the muon lateral spread distribution does not converge in the presented version of the tuning procedure. In particular, for the lateral distribution of inclined bundles, the MUPAGE parametric formulas do not approximate well the distribution as simulated with CORSIKA. The comparison between the CORSIKA and MUPAGE distributions for the lateral spread of two-muon bundles reaching the can is shown in Fig.4d. This distribution is dominated by vertical muons for which the MUPAGE functions with the adjusted parameters provide a good description. Several MUPAGE functions that aim to describe the multiple muon bundle energy dependence do not approximate the distributions resulting from CORSIKA with Sibyll 2.3d and GSF on the full parameter space, and thus, the parameter values are kept nominal and are not reported in the two tables above. The rate of events as a function of the bundle energy is plotted in Fig.4e. Although the MUPAGE parameters for the energy distribution of multi-muon bundles were not changed, the agreement between the tuned MUPAGE and CORSIKA has improved, thanks to the changes of parameter values connected with the multiplicity distributions. 5 Muon reconstruction performance To estimate the performance of the reconstruction algorithms, a comparison between the distributions of the cosine of the zenith angle from the MC reconstructed and generated events separately for the ORCA6 and ARCA6 detectors is done as shown in the top part of Fig.5. A discrepancy between the two distributions starts to emerge for events with cos θ<0.5 for ORCA6 and with cos θ<0.6 for ARCA6. Even though the KM3NeT angular resolution is at the subdegree level [26], the very steep dependence of the muon flux on the true cos θcauses the few well-reconstructed inclined events to be hidden by a small fraction of mis-reconstructed more vertical muons. The scatter plot for the true generated cos θversus the reconstructed one is given at the bottom of Fig.5. It can be noted also that all muons are assumed to be collinear with the bundle axis in MUPAGE. The median angular spread between muons in a bundle and the bundle axis is evaluated using CORSIKA. This value is around 0.14◦, with 90% of muons deviating from the bundle axis by less than 0.38◦. Hence, the parallel muon approximation should not affect the zenith distribution. The performance of the energy reconstruction for atmospheric muons is also evaluated. This analysis is carried out with data collected in a six-DU configuration for both KM3NeT/ORCA and KM3NeT/ARCA, this corresponds to 5% and 2.5% of the final design detector size, respectively. With the limited instrumented volume of the detectors, only a fraction of the muon energy is deposited within the sensitive volume of the telescopes. Therefore, the energy reconstruction algorithms fail to properly estimate the muon energy and the reconstructed spectrum does not match with the true MC one. The reconstruction performance is expected to improve for larger detector configurations. It is worth noticing that for muons with energies below 100 GeV radiative losses are negligible and their energy is reconstructed from the visible track length. The track length of 100 GeV muons is about 400m, which is longer than the sensitive size of the KM3NeT/ORCA detector, especially considering slightly-inclined muons that do not pass through the whole vertical length of the detector (∼300m). This limits the performance of the energy reconstruction for KM3NeT/ORCA. For the KM3NeT/ARCA detector muons 123 696 Page 16 of 19 Eur. Phys. J. C (2024) 84 :696 Fig. 9 Vertical muon spectrum at sea level from the CORSIKA simulation used in this work compared to two analytical models, Gaisser (orange line) and Bugaev (green line), and to data from ground-based experiments. The ratio of models and data to CORSIKA is shown on the bottom plot. Statistical and systematic uncertainties are shown as vertical error bars. The sea-level energy intervals for 90% of vertical muons (0.975 <cos θ<1) from reconstructed events for ARCA6 andORCA6areshownas coloured bands the Sibyll 2.1 (pre-LHC) hadronic interaction model [73] and two CR flux models, H3a [34] and GST [35]. Even though the models are different with respect to the ones used in this work, there is also an underestimation of the muon flux in the simulations with respect to the IceCube data at the level of 20% for vertical muons. The TeV muon energy spectrum at sea level predicted by the Sibyll 2.1 model is about 10% higher with respect to the one predicted by Sibyll 2.3d and the difference is almost energy independent; the H3a primary proton flux used in IceCube is higher by <10% than the one from GSF used in this work. Therefore, the IceCube result is compatiblewith theonefound inthisworkfor verticalmuons once the differences between the models are accounted for. In the IceCube analysis, the toy model simulation was fitted to data in order to get the power law index of the CR spectrum which provides the flat data/MC ratio. The spectral indices obtained for the trigger and high-quality data were −2.715 and −2.855, respectively. Both values differ from the one that provides the flat data/MC ratio for ORCA6 in this work, which is −2.607. Differences in the KM3NeT and IceCube results for the spectral index may be explained by the systematics of the two detectors and by the different hadronic models used. The ANTARES measurement of the zenith distribution of the atmospheric muon flux has been published in [21]. The detailed MC simulations with CORSIKA and pre-LHC hadronic interaction models underestimate the ANTARES data, although, the difference with respect to the MC simulation is within the systematic uncertainties for ANTARES. The underwater muon rate as a function of the zenith angle hasbeen also measuredby theBaikalCollaboration [22].The results show good agreement between the data and simulations. The hadronic interaction model used in these simulations is pre-LHC and the CR composition does not include the most recent direct measurements. 7.4 Comparison with the inclusive muon flux at sea level The measured discrepancy between data and simulations in KM3NeTcan be investigatedconsideringalso sea-level measurements of the muon flux. The high-energy muon flux can be measured at sea level up to several TeV. This threshold depends on the maximum detectable momentum of a muon, given by the relative momentum resolution of the spectrometer used to carry out such a measurement [74]. Several experiments were able to directly measure the absolute muon spectrum at sea level in the energy range relevant for KM3NeT: the L3+cosmic experiment [74], the Nottingham CR spectrometer [75,76], and the Kiel spectrographs [77]. TheverticalmuonfluxatsealevelresultingfromtheCORSIKA simulations used in this work is compared to the data from the experiments mentioned above, and to two analytical models, namely the Gaisser [2] and Bugaev [78] models. This comparison is shown in Fig.9. The CORSIKA simulations using Sibyll 2.3d and GSF underestimate the sea-level muon flux predicted by the analytical models at a level of ∼30%. The first L3+cosmic data point, the results from the Nottingham spectrometer obtained in 1984 [76], and the Kiel spectrograph data point are also above the CORSIKA predictions by ∼20–40%. This is compatible with the discrepancy seen by ORCA6 within the quoted uncertainties. The Nottingham data of 1968 [75] and the second L3+cosmic point are in agreement with the CORSIKA predictions and in disagreement with the two models considered here. 123 Eur. Phys. J. C (2024) 84 :696 Page 17 of 19 696 A similar discrepancy between the data for the inclusive flux of lower energy (GeV) muons at sea level and the MC simulations with Sibyll 2.3d and GSF is reported in [48]. This discrepancy is at a level of ∼20–30% for muons in the energy range from 1 GeV up to several hundred GeV. Muons contributing to the inclusive flux originate mainly from the first interactions in the EAS. Thus, they are mainly coming from CRs whose energies are higher by one order of magnitude with respect to the muon energy [71]. The inclusive TeV muon flux detected by the KM3NeT telescopes is also mainly driven by the first interactions. Therefore, it is worth highlighting that a similar muon deficit is observed for GeV muons coming from GeV CRs [48] and for TeV muons coming from TeV primaries as shown in this work. On the other hand, the deficit described in the Muon Puzzle is also for GeV muons but those muons are decay products of pions born after several stages of the shower cascade development. 7.5 On the atmospheric neutrino flux Atmospheric muons with energies of a few TeV are produced mainly in the decays of charged pions [71] (Figure 5) which also give rise to neutrinos. The neutrino energy is lower than that of the muon of the same pion decay by a factor ∼4[2] (Equation (6.14)). Hence, the TeV muon deficit observed in this work would correspond to the muon neutrino deficit in the ∼250 GeV energy domain. However, the inclusive neutrino flux for energies above 200 GeV is dominated by neutrinos coming from kaon decays [71] (Figure 5). Therefore, the deficit of atmospheric muons does not directly lead to the deficit of atmospheric neutrinos. Comparisons with atmospheric neutrino data can provide a complementary measurement for hadronic interaction studies. 7.6 On the first interactions probed by the Pierre Auger Observatory The Pierre Auger Observatory has reported measurements of the fluctuations in the number of muons in EAS produced by UHECRs [79]. These measurements agree with the MC simulationswiththerecenthigh-energyhadronicinteractions models, in particular with Sibyll 2.3d. The fluctuations in the muon number are believed to be mainly determined by the firstinteractionof theCRprimaries withtheatmosphere[80]. Since muons detected by the KM3NeT telescopes originate mainly from the first interactions in EAS, these detectors provide direct probes of such interactions. As demonstrated inthiswork,thereisadiscrepancybetweentheKM3NeTdata and the Sibyll 2.3d hadronic model. Therefore, it is important to notice that even though the hadronic models are able to describe the fluctuations of the muon number for UHECRs, they fail in the description of the absolute number of TeV muons originating from the first interactions for lower CR energies (TeV–few PeV range) as reported here. 8 Conclusions The most recent direct measurements of the primary CR flux describe the per-nucleus spectrum up to several hundreds of TeV, mainly driven by the CREAM and AMS-02 detectors [33]. Above these energies, the CR flux can only be measured indirectly through air shower observations, whose results are strongly dependent on various systematic uncertainties, most importantly those concerning the nature of hadronic interactions. Even though the knowledge on hadronic interaction has been significantly improved by the LHC measurements, the lack of accelerator data in the forward interaction region does not allow to describe CR air showers accurately. This leaves space for discrepancies, such as those observed recently in the EAS shower measurements. In this work, the high-energy muon contribution from CR showers is studied with the KM3NeT underwater neutrino detectors. A detailed simulation using the CORSIKA MC software has been used in this work to extract updated parameters for the fast underwater muon flux generator MUPAGE, namely using recent input for the hadronic interaction model (Sybil 2.3d) and the CR flux (the GSF model). This parameterisation allows for a comparison of simulations with data from the ORCA6 and ARCA6 detectors, and reveals a deficit in the simulations with respect to the data at the ∼40% level for TeV-energy atmospheric muons. This deficit is weakly dependent on the muon inclinations and thus on the muon pathlengthormuonenergy atsea level.Thedeficit iscompatible with the sea-level measurement and models of TeV-scale muon flux. The observed deficit of TeV muons indicates that the neutrino production in cosmic sources may be underestimated with respect to the flux of the accelerated nuclei and with respect to the gamma ray flux. Once the new hadronic interaction model is obtained as a solution to the observed atmospheric muon deficit, the gamma ray and neutrino production in cosmic ray accelerators could be revisited. An overview of the discrepancies in different muon energiescomingfromdifferentprimaryCRenergiesisreportedin thiswork.Thisprovidesadditionalinputsto the Muon Puzzle observed in the measurement of GeV muons from ultra-high energy CR air showers. The recent attempts to solve the discrepanciesintheMuonPuzzleshouldbeextendedtodescribe a broader phase space region where the issue is observed. Other muon kinematic properties important for understanding the discrepancy origin include muon bundle multiplicity, lateral spread, and underwater energy spectrum. Proper reconstruction of these observables is expected with the completed KM3NeT detectors. 123 696 Page 18 of 19 Eur. Phys. J. C (2024) 84 :696 Acknowledgements We would like to thank Rémi Adam for the discussions on the implications for the gamma ray emission mechanism. Funding The authors acknowledge the financial support of the funding agencies: Czech Science Foundation (GAˇ CR 24-12702S); Agence Nationale de la Recherche (contract ANR-15-CE31-0020), Centre NationaldelaRechercheScientifique(CNRS),CommissionEuropéenne (FEDER fund and Marie Curie Program), LabEx UnivEarthS (ANR10-LABX-0023andANR-18-IDEX-0001),ParisÎle-de-FranceRegion, France; Shota Rustaveli National Science Foundation of Georgia (SRNSFG,FR-22-13708),Georgia;TheGeneralSecretariatofResearch and Innovation (GSRI), Greece; Istituto Nazionale di Fisica Nucleare (INFN) and Ministero dell’Università e della Ricerca (MUR), through PRIN 2022 program (Grant PANTHEON 2022E2J4RK, Next Generation EU, Grant ALICA 2022A7ZC3K) and PON R&I program (Avviso n. 424 del 28 febbraio 2018, Progetto PACK-PIR01 00021), Italy; A. De Benedittis, R. Del Burgo, W. Idrissi Ibnsalih, A. Nayerhoda, G. Papalashvili, I. C. Rea, S. Santanastaso, A. Simonelli have been supported by the Italian Ministero dell’Università e della Ricerca (MUR), Progetto CIR01 00021 (Avviso n. 2595 del 24 dicembre 2019); Ministry of Higher Education, Scientific Research and Innovation, Morocco, and the Arab Fund for Economic and Social Development, Kuwait; Nederlandse organisatie voor Wetenschappelijk Onderzoek (NWO), the Netherlands; The National Science Centre, Poland (2021/41/N/ST2/01177); The grant “AstroCeNT: Particle Astrophysics Science and Technology Centre”, carried out within the International Research Agendas programme of the Foundation for Polish Science financed by the European Union under the European Regional Development Fund; National Authority for Scientific Research (ANCS), Romania; MCIN for PID2021-124591NB-C41, -C42, -C43, funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, for ASFAE/2022/014, ASFAE/2022 /023, with funding from the EU NextGenerationEU (PRTR-C17.I01), Generalitat Valenciana, and for CSIC-INFRA23013, Generalitat Valenciana for PROMETEO/2020/019, for Grant AST22_6.2 with funding from Consejería de Universidad, Investigación e Innovación and Gobierno de España and European Union - NextGenerationEU, for CIDEGENT/2018/034, /2019/043, /2020/049, /2021/23 and for GRISOLIAP/2021/192 and EU for MSC/101025085, Spain; The European Union’s Horizon 2020 Research and Innovation Programme (ChETEC-INFRA - Project no. 101008324). Data Availability Statement Data will be made available on reasonable request. [Author’s comment: The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.] CodeAvailabilityStatement Code/software will be made available on reasonable request. [Author’s comment: The code/software generated during and/or analysed during the current study is available from the corresponding author on reasonable request.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, whichpermitsuse,sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP3. References 1. V. F. Hess, Phys. Zeits. 13, (1912). https://inspirehep.net/literature/ 1623161 2. T.K. Gaisser, R. Engel, E. Resconi (Cambridge University Press, 2016). https://doi.org/10.1017/CBO9781139192194 3. B. Ó Fearraigh for the KM3NeT Collaboration, PoS(ICRC2021)1176. https://pos.sissa.it/395/1176/pdf 4. A. Aab et al., (the Pierre Auger Collaboration). Phys. Rev. D 91, 032003 (2015). https://doi.org/10.1103/PhysRevD.91.032003 5. A. Aab et al., (the Pierre Auger Collaboration). Eur. Phys. J. C 80, 751 (2020). https://doi.org/10.1140/epjc/s10052-020-8055-y 6. Y.A. Fomin et al., (the EAS-MSU Collaboration). Astropart. Phys. 92, 1 (2017). https://doi.org/10.1016/j.astropartphys.2017.04.001 7. J. G. Gonzalez for the IceCube Collaboration, EPJ Web Conf. 208, 03003 (2019). https://doi.org/10.1051/epjconf/201920803003 8. W.D. Apel et al., (the KASCADE-Grande Collaboration). Astropart. Phys. 95, 25 (2017). https://doi.org/10.1016/j. astropartphys.2017.07.001 9. A.G. Bogdanov et al., (the NEVOD-DECOR Collaboration). Astropart. Phys. 98, 13 (2018). https://doi.org/10.1016/j. astropartphys.2018.01.003 10. J.A. Bellido et al., (the SUGAR Collaboration). Phys. Rev. D 98(2), 023014 (2018). https://doi.org/10.1103/PhysRevD.98.023014 11. R.U. Abbasi et al., (the Telescope Array Collaboration). Phys. Rev. D 98(2), 022002 (2018). https://doi.org/10.1103/PhysRevD. 98.022002 12. A.V. Glushkov, A.V. Saburov (the Yakutsk EAS Collaboration), JETP Lett. 109(9), 559 (2019). https://doi.org/10.1134/ S0021364019090091 13. J. Albrecht et al., Astrophys. Space Sci. 367, 27 (2022). https://doi. org/10.1007/s10509-022-04054-5 14. S. Aiello et al., the KM3NeT Collaboration, Eur. Phys. J. C 83, 344 (2023). https://doi.org/10.1140/epjc/s10052-023-11401-5 15. S. Baur, H. Dembinski, M. Perlin, T. Pierog, R. Ulrich, K. Werner, Phys. Rev. D 107, 094031 (2023). https://doi.org/10.1103/ PhysRevD.107.094031 16. F. Riehn, R. Engel, A. Fedynitch, PoS(ICRC2023)429. https://doi. org/10.22323/1.444.0429 17. F.L. Villante, F. Vissani, Phys. Rev. D 78, 103007 (2008). https:// doi.org/10.1103/PhysRevD.78.103007 18. C. Mascaretti, F. Vissani, JCAP 08, 004 (2019). https://doi.org/10. 1088/1475-7516/2019/08/004 19. S.R. Kelner, F.A. Aharonian, V.V. Bugayov, Phys. Rev. D 74, 034018 (2006). https://doi.org/10.1103/PhysRevD.74.034018 20. E. Andres et al., (the AMANDA Collaboration), Astropart. Phys. 13, 1 (2000). https://doi.org/10.1016/S0927-6505(99)00092-4 21. J. A. Aguilar et al. (the ANTARES Collaboration), Astropart. Phys. 34, 179 (2010). https://doi.org/10.1016/j.astropartphys.2010.07. 001 22. V.A. Allakhverdyan et al., (the Baikal-GVD Collaboration), Eur. Phys. J. C 81, 1025 (2021). https://doi.org/10.1140/epjc/ s10052-021-09825-y 23. M.G. Aartsen et al., (the IceCube Collaboration), Astropart. Phys. 78, 1 (2016). https://doi.org/10.1016/j.astropartphys.2016.01.006 24. S. Aiello et al., (the NEMO Collaboration), Astropart. Phys. 33, 263 (2010). https://doi.org/10.1016/j.astropartphys.2010.02.009 25. G. Aggouras et al., (the NESTOR Collaboration), Astropart. Phys. 23, 377 (2005). https://doi.org/10.1016/j.astropartphys.2005.02. 001 26. S. Adrián-Martínez et al., (the KM3NeT Collaboration), J. Phys. GNucl.Part.Phys.43, 084001 (2016). https://doi.org/10.1088/ 0954-3899/43/8/084001 123 Eur. Phys. J. C (2024) 84 :696 Page 19 of 19 696 27. S. Aiello et al., (the KM3NeT Collaboration), Comput. Phys. Commun. 296, 109036 (2024). https://doi.org/10.1016/j.cpc.2023. 109036 28. S. Aiello et al., (the KM3NeT Collaboration), JINST 17, P07038 (2022). https://doi.org/10.1088/1748-0221/17/07/P07038 29. G. Riccobene for the KM3NeT Collaboration, EPJ Web of Conferences 207, 07005 (2019). https://doi.org/10.1051/epjconf/ 201920707005 30. L. Bailly-Salins for the KM3NeT Collaboration, PoS(ICRC2023)218. https://pos.sissa.it/444/218/pdf 31. K. Melis, A. Heijboer, M. de Jong for the KM3NeT Collaboration, PoS(ICRC2017)950. https://doi.org/10.22323/1.301.0950 32. D. Heck et al., FZKA-6019, FZKA Karlsruhe (1998). https:// inspirehep.net/literature/469835 33. H.P. Dembinski et al., PoS(ICRC2017)533. https://doi.org/10. 22323/1.301.0533 34. T.K. Gaisser, Astropart. Phys. 35, 12 (2012). https://doi.org/10. 1016/j.astropartphys.2012.02.010 35. T.K. Gaisser, T. Stanev, S. Tilav, Front. Phys. 8, 748–758 (2013). https://doi.org/10.1007/s11467-013-0319-7 36. J. T. Emmert et al., Earth Space Sci. 8, e2020EA001321 (2021). https://doi.org/10.1029/2020EA001321 37. T. Lang et al., Phys. Rev. C 93(1), 014901 (2016). https://doi.org/ 10.1103/PhysRevC.93.014901 38. F. Riehn, R. Engel, A. Fedynitch, T.K. Gaisser, T. Stanev, Phys. Rev. D 102, 063002 (2020). https://doi.org/10.1103/PhysRevD. 102.063002 39. S. Aiello et al., (the KM3NeT Collaboration), Comput. Phys. Commun. 256, 107477 (2020). https://doi.org/10.1016/j.cpc.2020. 107477 40. C. Andreopoulos et al., Nucl. Instrum. Methods Phys. Res. Sect. A 614, 87 (2010). https://doi.org/10.1016/j.nima.2009.12.009 41. J.H. Koehne, K. Frantzen, M. Schmitz, T. Fuchs, W. Rhode, D. Chirkin, J.B. Tjus, Comput. Phys. Commun. 184, 2070 (2013). https://doi.org/10.1016/j.cpc.2013.04.001 42. D.E. Groom, N.V. Mokhov, S.I. Striganov, At. Data Nucl. Data Tables 78, 183 (2001). https://doi.org/10.1006/adnd.2001.0861 43. R.M. Sternheimer, Phys. Rev. 88, 851 (1952). https://doi.org/10. 1103/PhysRev.89.1309.2 44. S. Kelner et al., MEphI Report No. MEPhI 024-95 (1995). https:// inspirehep.net/literature/400124 45. R.P. Kokoulin, A.A. Petrukhin, Proceedings of the XII Int. Conf. on Cosmic Rays 6(1971). https://inis.iaea.org/ search/searchsinglerecord.aspx?recordsFor=SingleRecord& RN=4049614 46. H. Abramowicz, A. Levy, arXiv:hep-ph/9712415 47. L.D. Landau, I.J. Pomeranchuk, Dokl. Akad. Nauk SSSR 92, 535 (1953). http://cds.cern.ch/record/436540 48. A. Fedynitch, W. Woodley, M.-C. Piro, Astrophys. J. 928,27 (2022). https://doi.org/10.3847/1538-4357/ac5027 49. A. Okada, Astropart. Phys. 2, 4 (1994). https://doi.org/10.1016/ 0927-6505(94)90028-0 50. S. Adrián-Martínez et al., (the ANTARES Collaboration), Astropart. Phys. 35, 9 (2012). https://doi.org/10.1016/j. astropartphys.2011.12.003 51. S. Sparnocchia, G.P. Gasparini, K. Schroeder, M. Borghini, Nucl. Instrum. Methods Phys. Res. Sect. A 626–627 (2011). https://doi. org/10.1016/j.nima.2010.06.231 52. F.J. Millero, R. Feistel, D.G. Wright, T.J. McDougall, Deep Sea Res. Part I(55), 50–72 (2008). https://doi.org/10.1016/j.dsr.2007. 10.001 53. S.H.D. Haddock, M.A. Moline, J. Case, Annu. Rev. Mar. Sci. 2, 443–493 (2010). https://doi.org/10.1146/ annurev-marine-120308-081028 54. A. Albert et al., (the ANTARES Collaboration). JCAP 01, 064 (2021). https://doi.org/10.1088/1475-7516/2021/01/064 55. Y. Becherini, A. Margiotta, M. Sioli, M. Spurio, Astropart. Phys. 25, 1 (2006). https://doi.org/10.1016/j.astropartphys.2005.10.005 56. G. Carminati, A. Margiotta, M. Spurio, Comput. Phys. Commun. 179, 915 (2008). https://doi.org/10.1016/j.cpc.2008.07.014 57. C.Fortietal.,Phys.Rev.D42, 3668 (1990). https://doi.org/10. 1103/PhysRevD.42.3668 58. M. Ambrosio et al., (the MACRO Collaboration), Phys. Rev. D 60, 032001 (1999). https://doi.org/10.1103/PhysRevD.60.032001 59. O. Adriani et al., (the CALET Collaboration), Phys. Rev. Lett. 129, 101102 (2022). https://doi.org/10.1103/PhysRevLett.129.101102 60. F. Alemanno et al., (the DAMPE Collaboration), Phys. Rev. Lett. 126, 201102 (2021). https://doi.org/10.1103/PhysRevLett. 126.201102 61. G.H. Choi et al., (the ISS-CREAM Collaboration), Astrophys. J. 940, 107 (2022). https://doi.org/10.3847/1538-4357/ac9d2c 62. A. Albert et al., (the HAWC Collaboration), Phys. Rev. D 105, 063021 (2022). https://doi.org/10.1103/PhysRevD.105.063021 63. G. Riccobene et al., (the NEMO Collaboration), Astropart. Phys. 27, 1-9 (2007) https://doi.org/10.1016/j.astropartphys.2006.08. 006 64. X. Durrieu de Madron et al., J. Geophys. Res. Oceans 122, 2291– 2318 (2017). https://doi.org/10.1002/2016JC012062 65. M. Ageron et al., (the KM3NeT Collaboration), Eur. Phys. J. C 80, 99 (2020). https://doi.org/10.1140/epjc/s10052-020-7629-z 66. A. De Benedittis for the KM3NeT Collaboration, Phys. Sci. Forum 8(1), 44 (2023). https://doi.org/10.3390/psf2023008044 67. A. Fedynitch, R. Engel, T. K. Gaisser, F. Riehn, and T. Stanev, PoS(ICRC2015)1129. https://doi.org/10.22323/1.236.1129 68. R. Ulrich, A, Fedynitch, T, Pierog, M. Reininghaus, F. Riehn for the CORSIKA8 Collaboration, PoS(ICRC2021)474. https://doi. org/10.22323/1.395.0474 69. T. Pierog, I. Karpenko, J.M. Katzy, E. Yatsenko, K. Werner, Phys. Rev. C 92, 034906 (2015). https://doi.org/10.1103/PhysRevC.92. 034906 70. S. Ostapchenko, Phys. Rev. D 83, 014018 (2011). https://doi.org/ 10.1103/PhysRevD.83.014018 71. A. Fedynitch, F. Riehn, R. Engel, T.K. Gaisser, T. Stanev, Phys. Rev. D 100, 103018 (2019). https://doi.org/10.1103/PhysRevD. 100.103018 72. J. Mulder, R. Bruijn for the KM3NeT Collaboration, PoS(ICRC2023)355. https://doi.org/10.22323/1.444.0355 73. E.-J. Ahn, R. Engel, T.K. Gaisser, P. Lipari, T. Stanev, Phys. Rev. D 80, 094003 (2009). https://doi.org/10.1103/PhysRevD.80.094003 74. P. Achard et al., (the L3 + C Collaboration). Phys. Lett. B 598,15 (2004). https://doi.org/10.1016/j.physletb.2004.08.003 75. S.R. Baber, W.F. Nash, B.C. Rastin, Nucl. Phys. B 4, 539 (1968). https://doi.org/10.1016/0550-3213(68)90112-0 76. B.C. Rastin, J. Phys. G Nucl. Phys. 10, 1609 (1984). https://doi. org/10.1088/0305-4616/10/11/017 77. O.C. Allkofer, K. Carstensen, D.W. Dau, Phys. Lett. B 36, 425 (1971). https://doi.org/10.1016/0370-2693(71)90741-6 78. E.V. Bugaev, A. Misaki, V.A. Naumov, T.S. Sinegovskaya, S.I. Sinegovsky, N. Takahashi, Phys. Rev. D 58, 054001 (1998). https:// doi.org/10.1103/PhysRevD.58.054001 79. A.Aabetal.,(thePierreAugerCollaboration),Phys.Rev.Lett.126, 152002 (2021). https://doi.org/10.1103/PhysRevLett.126.152002 80. L. Cazon, R. Conceição, F. Riehn, Phys. Lett. B 784, 68 (2018). https://doi.org/10.1016/j.physletb.2018.07.026 123