Prospects for combined analyses of hadronic emission from gamma-ray sources in the Milky Way with CTA and KM3NeT

Eur. Phys. J. C (2024) 84:112 https://doi.org/10.1140/epjc/s10052-023-12279-z Regular Article - Experimental Physics Prospects for combined analyses of hadronic emission from γ-ray sources in the Milky Way with CTA and KM3NeT T. Unbehaun, L. Mohrmann, S. Funk (members of the CTA Consortium) and the KM3NeT Collaboration,a Received: 8 August 2023 / Accepted: 21 November 2023 / Published online: 2 February 2024 © The Author(s) 2024 Abstract TheCherenkovTelescopeArrayandtheKM3NeT neutrino telescopes are major upcoming facilities in the fields of γ-ray and neutrino astronomy, respectively. Possible simultaneous production of γrays and neutrinos in astrophysical accelerators of cosmic-ray nuclei motivates a combination of their data. We assess the potential of a combined analysis of CTA and KM3NeT data to determine the contribution of hadronic emission processes in known Galactic γ-ray emitters, comparing this result to the cases of two separate analyses. In doing so, we demonstrate the capability of Gammapy, an open-source software package for the analysisof γ-raydata, toalso processdata fromneutrino telescopes. For a selection of prototypical γ-ray sources within our Galaxy, we obtain models for primary proton and electron spectra in the hadronic and leptonic emission scenario, respectively, by fitting published γ-ray spectra. Using these models and instrument response functions for both detectors, we employ the Gammapy package to generate pseudo data sets, where we assume 200h of CTA observations and 10 years of KM3NeT detector operation. We then apply a three-dimensional binned likelihood analysis to these data sets, separately for each instrument and jointly for both. We find that the largest benefit of the combined analysis lies in the possibility of a consistent modelling of the γ-ray and neutrino emission. Assuming a purely leptonic scenario as input, we obtain, for the most favourable source, an average expected 68% credible interval that constrains the contribution of hadronic processes to the observed γ-ray emission to below 15%. 1 Introduction We live in the era of multi-messenger astrophysics [1]. Long anticipated, this paradigm advocates that unique insights about astrophysical objects and processes may be gained ae-mail: [email protected];[email protected]; [email protected] through the joint consideration of information carried by different messengers: photons, neutrinos, cosmic rays (CRs), and gravitational waves. In the past years, it has truly come to fruition, yielding the first promising results [2,3]. AstrophysicalobjectsintheMilkyWayarenotexpectedto producegravitationalwavesdetectablebycurrent-generation instruments (but by next-generation detectors, see [4]). GalacticCRsprovideimportantenergeticconstraintsontheir sourcepopulation, butcannotbeusedtodirectlystudyGalactic objects because they are deflected by magnetic fields. Hence,oneneedstoresorttophotonsandneutrinostoemploy multi-messenger astrophysics for the study of individual Galactic objects. Indeed, besides its conceptual attractiveness, the combined study of very-high-energy (VHE; E>100 GeV) γ rays and TeV–PeV neutrinos from Galactic sources is well motivated: they are expected to be produced simultaneously in‘hadronicaccelerators’,whereacceleratedCRnucleiinteract with ambient gas producing pions (and other mesons) that subsequently decay into γrays and neutrinos. In the following, this process is labelled ‘PD’, for pion decay. The situation is different in ‘leptonic accelerators’, where Inverse Compton(IC)up-scattering ofphotonsby CR electronsleads to VHE γ-ray emission without neutrino production. This implies that the detection of high-energy neutrinos coming from astrophysical objects is decisive in identifying them as hadronicaccelerators[5].Nevertheless,observationsofVHE γ-ray emission from the same objects are indispensable, as they provide a much higher detection sensitivity and allow a measurement of the spectrum and morphology of the emission in greater detail. This, in turn, enables a realistic estimation of the expected flux of neutrinos, and the possibility of detecting it with current (or planned) neutrino telescopes. The latter exercise has been carried out by various authors in the past (see, e.g. [6–12]). In this work, we focus on the Cherenkov Telescope Array (CTA) [13,14] and the KM3NeT neutrino telescopes [15], as major upcoming facilities for VHE γ-ray and neutrino astronomy, respectively. CTA will be built at La Palma, 123 112 Page 2 of 19 Eur. Phys. J. C (2024) 84 :112 Spain, and Paranal, Chile, covering the Northern and Southern sky, respectively. Our prime targets of interest being Galactic γ-ray sources, which are more easily observed from the Southern hemisphere, we consider only the site in Chile (CTA-South) for our study. At this site, the installation of ∼50 imaging atmospheric Cherenkov telescopes (IACTs) of two different sizes, covering the energy range between 100 GeV and 300 TeV, is foreseen. IACTs detect γrays by measuringthefaintflashof Cherenkovlightthat isemittedby secondary particles in the air shower that is launched when theprimaryγray hits the atmosphere of the Earth. Compared to current-generation arrays of IACTs, CTA is projected to provide a ten-fold increase in sensitivity. KM3NeT is a research infrastructure in the Mediterranean, consisting of neutrino telescopes installed in the deep sea at different locations. The ‘Oscillation Research with Cosmics in the Abyss’ (ORCA) detector, with its dense instrumentation, will focus on the study of neutrino properties measuring atmospheric neutrinos [16]. Here, we only consider the ‘Astroparticle Research with Cosmics in the Abyss’(ARCA)detector, whichtargetsthedetection ofhighenergyastrophysicalneutrinoswith TeV-PeVenergies.Hereafter, we will use the term ‘KM3NeT’ to refer to the ARCA telescope. It is currently under construction off-shore Sicily, Italy, and will ultimately consist of two building blocks comprising 115 vertical detection units each. Each detection unit carries 18 digital optical modules (DOMs) [17–19] with a vertical spacing of 36 m and is about 700 m tall. The units are arranged on a grid with about 90 m spacing between them. The DOMs contain light sensors that detect Cherenkov light radiated by secondary particles created in interactions of high-energy neutrinos in or near the detector. Of particular interest for this work are muons created in charged-current interactions of muon neutrinos, as their long propagation distances of up to several km in the water allow a precise reconstructionofthedirectionoftheincomingneutrino.Compared to the largest existing neutrino telescope, the IceCube Neutrino Observatory [20,21], KM3NeT utilises water instead of ice as detector medium. This reduces scattering of the Cherenkov light and is expected to lead to an improved angular resolution [15]. Its location in the Northern hemisphere implies that neutrinos potentially emitted by many Galactic γ-ray sources would reach KM3NeT through the Earth, which is advantageous for the suppression of atmospheric muon background events [22]. Therefore, the combination of CTA-South and KM3NeT data for the study of Galactic objects appears very natural. The discovery potential for extended Galactic sources by KM3NeT, in relation to the constraining power of CTA, has beeninvestigated in[23].Inthis paper,wedemonstrate howa combined analysis of CTA and KM3NeT data can be used to constrain physical properties of γ-ray sources in our Galaxy, with special attention to the contribution of hadronic emission processes. To this end, using Monte Carlo simulations as input, we have prepared instrument response functions (IRFs) for the KM3NeT detector and stored them in the same ‘GADF’ data format1used for the publicly available CTA IRFs. Then we have employed the Gammapy2package (version 0.17; [25,26]) to generate pseudo data sets based on the obtained IRFs and to perform a joint 3D likelihood fit on these data sets. Though Gammapy is still under development and application of the 3D likelihood method in IACT data analysis represents a recent approach, both have been validated using a public IACT data set [27]. Weapplytheanalysistoaselectionofprototypicalsources thatarepromising candidates fortheemissionof high-energy neutrinos (see Table 1). We motivate our choice briefly in the following: – Vela X is a pulsar wind nebula bright in γrays [28], associated with the well-known Vela pulsar. Pulsars being copious producers of electrons and positrons, the γ-ray emission from Vela X is expected to be largely due to IC emission, and no associated neutrino emission is expected. However, also mixed, lepto-hadronic models have been considered in the past (e.g. [29,30]), leaving room for some neutrino emission from this source. – RX J1713.7−3946 is a shell-type supernova remnant that emits γrays in excess of 10 TeV [31]. The emission has been modelled according to leptonic, hadronic, and mixed scenarios (see e.g. [32,33] for recent studies). The observation (but also non-observation) of neutrinos from RX J1713.7−3946 could therefore yield important clues about acceleration processes at play. – Westerlund 1 is the most massive young stellar cluster in the Milky Way [34], and is considered the most likely counterpart of the VHE γ-ray source HESS J1646−458 [35,36]. Massive stellar clusters have recently been hypothesized as PeVatrons [37], making Westerlund 1 a good candidate for high-energy neutrino emission. – eHWC J1907+063, also known as HESS J1908+063, is an unidentified γ-ray source [38] that was recently detected above energies of 100TeV by the High Altitude Water Cherenkov Observatory (HAWC) [39]aswell as by the Large High Altitude Air Shower Observatory (LHAASO) [40]. Like in the case of RX J1713.7−3946, observations with neutrino telescopes could help to constrain the nature of the source. Our selection of sources is neither a complete list of promising targets for the emission of neutrinos in our Galaxy, nor does it comprise only the most promising ones. Rather, we have aimed for a selection of different types of γ-ray sources 1See https://gamma-astro-data-formats.readthedocs.io and [24]. 2https://gammapy.org. 123 Eur. Phys. J. C (2024) 84 :112 Page 3 of 19 112 Fig. 1 Source visibility with KM3NeT. Shown is the fraction of time that each source is visible under a zenith angle θ, over the course of one year. The green-shaded area indicates the zenith angle range used in the analysis and the percentage value in parentheses specifies the fraction that each source is visible within this range that furthermore exhibit favourable locations for the observation with CTA-South and KM3NeT. Figure1shows the visibility of all sources for KM3NeT, as a function of the local zenith angle.3For CTA, within one year, the sources are observable above an altitude angle of 50◦for a maximum time of ∼400 h (Vela X), ∼510 h (RX J1713.7−3946), ∼500 h (Westerlund 1), and ∼380 h (eHWC J1907+063). Finally, we note that the IRFs we derived and used in this work are not representative of the final sensitivity of CTA and KM3NeT. They are based on preliminary simulations and event selections that will likely be improved in the future. In particular, the IRFs for KM3NeT are based on a pointsource analysis as presented in [41]. This does not impact the conclusions drawn in this work, which are focused more on the conceptual benefits of a combined analysis rather than on numerical results. The paper is structured as follows. In Sect. 2,weintroduce our methodology: the computation of IRFs (Sect. 2.1), the preparation of input models (Sect. 2.2), the generation of pseudo data sets (Sect. 2.3), the combined likelihood analy- 3The zenith angle θrefers to the location of the source. For θ=0◦ the source is above the detector and produces vertically down-going neutrinos, while for θ=180◦the source is located on the opposite side of the Earth and produces vertically up-going neutrinos. sis (Sect. 2.4), and the derivation of constraints on hadronic contributions (Sect. 2.5). The results of the analysis are then presented and discussed in Sect. 3, before we conclude the paper in Sect. 4. 2 Methodology 2.1 Instrument response functions Given a physical source model, the IRFs for each experiment allow us to compute how this source would appear in thedetector.Inour case,the relevant IRFscomprisetheeffective area, the energy dispersion, and the point spread function (PSF), reflecting the sensitivity, energy resolution, and angular resolution of the instrument, respectively. Additionally, background templates that yield the expected number of mis-classified background events, arising from CR-induced atmospheric air showers, are necessary. For IACTs, the IRFs are typically stored as a function of the true γ-ray energy and of the true angular offset of the events with respect to the pointing direction of the telescopes (‘offset angle’). The IRFs also depend on the angle with respect to zenith of the pointing position of the telescopes. However, because the variation within a specific observation run (of typically 30 min duration and a field of view of ∼5◦×5◦) is small, IRFs for the average zenith angle of the observation run are commonly employed. For CTA we use the publicly available ‘Prod 5’ IRFs4for the southern array at 20◦zenith angle and averaged over azimuth angle [42]. The IRFs for KM3NeT have been custom-generated for this study, as detailed in the following section. 2.1.1 Generation of KM3NeT IRFs The KM3NeT IRFs are based on extensive simulations of neutrinos and anti-neutrinos5that interact in or near the detector. We focus on charged-current interactions of muon neutrinos only, as they give rise to long-range muons that appearascharacteristic, track-likeevents inthedetector.This leads to a good angular resolution (<0.3◦for energies >10 TeV), which helps in suppressing background events (but see also [43]). The neutrino events have been simulated with the gSeaGen software [44] based on the GENIE neutrino generator [45], which allows the simulation of interactions of all neutrino flavours in the media around the detector. On the level of a single optical module, the decay of 40Kas well as bioluminescence are relevant sources of noise. Due to the design of the optical modules, which contain multi- 4See https://www.cta-observatory.org/science/ctao-performance. 5Hereafter, we will use the term ‘neutrinos’ to refer to both neutrinos and anti-neutrinos, unless explicitly stated otherwise. 123 112 Page 4 of 19 Eur. Phys. J. C (2024) 84 :112 Table 1 Galactic γ-ray sources investigated in this work Designation Type Spatial model r(deg) Declination (deg) Distance (kpc) References Vela X PWN disk 0.8 −45.19 0.29 [28] RX J1713.7−3946 SNR disk 0.6 −39.69 1 [31] Westerlund 1 SC disk 1.1 −45.85 3.9 [35] eHWC J1907+063 UNID Gaussian 0.67 +06.18 2.37 [39] ‘Spatial model’ specifies which type of spatial model is used in the analysis (cf. Sect. 2.3). rdenotes the radius of the disk in case of a disk model and the width of the Gaussian in case of a Gaussian model ‘Type’ refers to the source type, PWN pulsar wind nebula, SNR supernova remnant, SC stellar cluster, UNID unidentified ple photo-sensors each, these backgrounds can however be suppressed very efficiently by requiring a local coincidence between the photo-sensors [18]. On the analysis level, two types of background events are relevant for KM3NeT: neutrinos and muons, both created in CR-induced atmospheric air showers. Both can be further classified as ‘conventional’ – resulting mostly from the decays of pions and kaons – and ‘prompt’ – resulting from the decays of heavy hadrons and light vector mesons. The former exhibit a steeper energy spectrum,becausetheirparentparticleshaveanon-negligible chance to re-interact with air molecules, rather than to decay. To predict the rate of atmospheric neutrino events, we use the ‘HKKMS’ model [46] for conventional neutrinos and the ‘ERS’ model [47] for prompt neutrinos. Both models are based on outdated parametrisations of the primary CR flux and have been corrected as described in [48] to conform with the ‘H3a’ parametrisation from [49]. We note that there is also the possibility of a CR composition around the ‘knee’ feature in the CR spectrum that is heavier than predicted by the H3a model. This scenario is discussed in more detail in [50], but not investigated further here. The resulting event rates of conventional and prompt atmosphericneutrinos,integratedoverrelevantzenithangles, are shown in Fig. 2. The background of atmospheric muons has been estimated using dedicated simulations of muons using the MUPAGE package [51–53]. While there is in principle also a potential background due to diffuse astrophysical neutrinos not connected to the studied source itself [54], this background can be safely neglected here. Allsimulatedevents arereconstructedusing a trackreconstruction algorithm and subsequently undergo a selection procedure based on the reconstruction quality and a classification algorithm using boosted decision trees (BDTs), as detailed in [41]. In order to suppress the background of atmospheric muons, which always arrive from above the detector, we restrict the analysis region to reconstructed zenith angles θreco >80◦. Thus, only very few atmospheric muon events remain in the final sample, almost all concentrated close to the horizon region (80◦<θ reco <90◦), see the blue histogram in Fig. 2. In order to avoid interpolation problems due to empty bins in the histogram, we fit a spline curve to the histogram and use this curve to predict the expected rate of atmospheric muon events (black line). We note that due to insufficient simulation statistics, the exact shape of the distribution at energies below 1 TeV should not be trusted. Because the muon background is sub-dominant compared to the atmospheric neutrino background by several orders of magnitude at these energies, however, this does not affect our results. The KM3NeT IRFs mainly depend on neutrino energy and zenith angle. To be able to store the IRFs in the (IACT- centred) GADF data format, we utilise the offset-angle axis defined there to describe the dependence of the KM3NeT IRFs on the zenith angle. The IRFs are then generated by creating histograms of the appropriate event properties (e.g. the angle between the reconstructed and true neutrino direction in case of the PSF) and applying corresponding normalisation factors. For the effective area IRF, we use 48 logarithmic bins in true energy between 100 GeV and 100 PeV and 12 zenith angle bins linear in true cos(θ). The energy dispersion and PSF IRFs – featuring one more dimension than the effective area – are created with twice the bin size, in order to ensure sufficient statistics in each bin. In the analysis, a linear interpolation between the individual bins is performed. IRFs for neutrinos and anti-neutrinos are derived separately and subsequently averaged, assuming equipartition of the source flux between the two. 2.1.2 Comparison of IRFs In this section, we provide a comparison of the IRFs of CTA and KM3NeT. While the KM3NeT IRFs generated by ourselves are defined up to an energy of 100 PeV, the public CTA IRFs are valid up to an energy of only ∼300 TeV. This is because, given its limited duty cycle and need for pointing, CTA is not expected to be able to effectively measure fluxes beyond that energy. Figure 3shows a comparison of the effective areas. The effective area of CTA rises sharply at the threshold energy of the instrument (around 0.1 TeV), before the curve gradually flattens as γrays are detected more and more efficiently. The effective area of KM3NeT for neutrinos is much lower 123 Eur. Phys. J. C (2024) 84 :112 Page 5 of 19 112 Fig. 2 Background event rates in KM3NeT as a function of reconstructed energy Ereco. The atmospheric neutrino rates are integrated over all zenith angles in the analysis region (i.e. θreco >80◦). The atmospheric muon rate is shown for the zenith angle bin 80◦–90◦only, since it is completely negligible for larger angles. The black line displays the smoothed curve used in the analysis Fig. 3 Comparison of effective areas. The CTA effective area is shown for a zenith angle of θ=20◦and an offset angle from the pointing direction of ϑ=1◦. For KM3NeT, average effective areas for the full analysis region as well as for different sub-ranges in zenith angle are shown than that of CTA for γrays because of the low interaction probability of neutrinos. The increase in neutrino effective area with increasing energy reflects a corresponding increase of the interaction cross section and detection efficiency. At the highest energies, the interaction cross section becomes large enough for the Earth to become opaque to neutrinos, leading to a decrease in the effective area for neutrinos that traverse large amounts of matter (green dashed-dotted and red dotted line in Fig. 3). Fig. 4 Comparison of the directional reconstruction accuracy of CTA and KM3NeT. Shown are the 68% and 95% containment radii of the PSF as a function of the true γ-ray/neutrino energy. Possible differences to previous publications may arise from the finite binning of the PSF applied here The angular resolutions of the two instruments – here expressed in terms of the 50%, 68%, and 95% quantiles of the respective PSFs – are compared in Fig. 4. While the angular resolution of CTA is clearly superior to that of KM3NeT, the selection of track-like events for the KM3NeT analysis still leads to a median resolution of better than 0.3◦above ∼10 TeV. The PSF strongly affects the sensitivity to pointlike or marginally extended sources, as the contribution of background events increases quadratically with the radius of the source after PSF convolution. For a comparison of the KM3NeT PSF with that of the IceCube neutrino telescope, see for example [41]. Figure 5provides a comparison of the energy resolution of the two instruments, here indicated by the 10%, 50%, and 90% quantiles of the ratio between reconstructed and true energy. In the case of KM3NeT, muons created in muon neutrino interactions may lose energy before entering the detector or carry away energy when leaving it. Because only the energy deposited inside the detector can reliably be estimated, the resulting reconstructed energy is on average lower than the true neutrino energy. This effect becomes more and more prominent as the neutrino energy – and hence the track length of the resulting muon – increases. 2.2 Input models Input models are needed for the likelihood analysis of each analysed source (cf. Table 1). During the analysis procedure a spatial model is chosen for each source according to the description found in the literature: either a uniform disk or a two-dimensional Gaussian model. These spatial models are not varied or fitted (i.e. remain fixed) during the entire analysis procedure. In addition, two different spectral models have been considered for each source in order to study the 123 112 Page 6 of 19 Eur. Phys. J. C (2024) 84 :112 Fig. 5 Comparison of the energy reconstruction accuracy of CTA and KM3NeT. Shown are the 10%, 50%, and 90% quantiles of the ratio between the reconstructed and true γ-ray/neutrino energy, as a function of the true energy fraction of hadronic γ-ray emission: an IC model, assuming a purely leptonic emission and a PD model, to describe a purely hadronic emission. Both models have independently beenfittedtopublishedγ-rayspectraof thesources (cf.references in Table 1). For the IC model, we have used the implementation of the InverseCompton model in the naima package [55], which provides one-zone, time-independent radiative models. For the PD model, we have implemented a corresponding model based on the parametrisation in [56]6. For both the IC and PD models, we assume a power-law model with an exponential cut-off for the primary electron and proton distributions, Φ(E)=A·E E0−Γ exp −E Ecut β,(1) where Adenotes the amplitude, E0the reference energy, Γ the spectral index, Ecut the cut-off energy, and βthe cut-off strength. For each source and for both the IC and PD models, we have adjusted the amplitude, spectral index, and cut-off energy using a simple χ2fit, keeping the reference energy and cut-off strength fixed (at E0=10 TeV and β=1, respectively). The fit for Vela X is shown in Fig. 6, whereas those for the other sources can be found in Appendix A. Both models describe the γ-ray flux equally well, illustrating the difficulty to distinguish between the leptonic and hadronic scenariosbasedonγ-raydata alone.Wenotethattheaddition 6We are aware of more recent parametrisations such as that provided by [57],whichisalsousedinthePionDecay model implementation in naima. That parametrisation, however, does not provide a prediction for the expected neutrino flux, required for our analysis. The focus of this study lying on the technical feasibility of a joint γ-ray/neutrino analysis, the choice of parametrisation for the hadronic model is not relevant here. Fig. 6 Fit of hadronic (PD) and leptonic (IC) input models for Vela X. The muon-neutrino prediction based on the best-fit PD model is shown as dashed line. The data points are taken from [28], and based on 53 h of H.E.S.S. observations of lower-energy radio or X-ray data would further constrain the fit, but regard this as beyond the scope of this work. 2.3 Generation of pseudo data sets With the IRFs and input models in hand, we generated 100 pseudodata setsfor eachsourceandeachinstrument,bothfor the PD and IC models. We note that for both instruments, we do not take into account diffuse astrophysical γ-ray or neutrino emission that is unrelated to the source itself. For the CTA data sets we used an analysis setup with 16 energy bins per decade between 0.1 TeV and 154 TeV and spatial bins of 0.02◦×0.02◦size. For each pseudo data set, we assumed a total observation time of 200h, split equally between four pointing positions with 1◦offset with respect to the source position. The predicted number of source and background events are summed for each pixel and Poisson-distributed random counts are drawn based on those values. As an example, in Figs. 7and 8, we show projections of one generated pseudo CTA data set based on the PD model for the source Vela X onto the spatial and energy axes, respectively. A similar procedure has been used to generate the KM3NeT pseudo data sets, for which we assume a total detector operation time of 10 years. For each zenith angle bin (cf. Fig. 1), we generate an observation set evaluating the associated IRFs, according to the corresponding fraction of the total observation time. This results in 6 or 7 observation setsforeachsource,dependingonthevisibility.Theexposure and expected background for every data set are computed in equatorial coordinates by integrating over time the respective IRFs (defined in terms of the local zenith angle). Data are binned using spatial pixels of 0.1◦×0.1◦size and 4 bins per decade in energy, between 100 GeV and 1 PeV. The IRFs are evaluated using a finer binning in energy, with 16 bins per decade between 100 GeV and 10 PeV. Several tests have 123 Eur. Phys. J. C (2024) 84 :112 Page 7 of 19 112 Fig. 7 Counts map of a pseudo CTA data set based on the PD model for Vela X with 200h of observation time. The counts are Poissonrandomised based on the model prediction for Vela X (modelled as a disk with radius 0.8◦) and the residual hadronic background, summed over all energies and smoothed with a 0.05◦Gaussian. The blue circle and ‘×’ markers denote the source and pointing positions, respectively. The white dashed circle shows the source region for which the counts spectra shown in Fig. 8have been extracted Fig. 8 Counts spectra for the Vela X CTA PD data set, extracted for a region encompassing the source (cf. Fig. 7). The coloured lines denote the number of predicted counts within the source region for an observation time of 200h. The black data points visualise one random Poisson realisation, drawn from the model predictions been performed with finer zenith, energy, and spatial binning, yielding consistent results and no significant change in sensitivity. An example of a KM3NeT pseudo data set is visualised in Figs. 9and 10, respectively. 2.4 Likelihood analysis The analysis is performed using the binned likelihood formalism implemented in the Gammapy package.7Leptonic 7We note that analyses carried out with the analysis tools normally employed by the KM3NeT Collaboration are typically performed using Fig. 9 Counts map of a pseudo KM3NeT data set based on the PD model for Vela X with 10 years of observation time. The counts are Poisson-randomised based on the model prediction for Vela X, summed over all energies and smoothed with a 0.25◦Gaussian. The white dashed circle shows the source region for which the counts spectra shown in Fig. 10 have been extracted (same region as in Fig. 7) Fig. 10 Counts spectra for the Vela X KM3NeT PD data set, extracted for a region encompassing the source (cf. Fig. 9). The coloured lines denote the number of predicted counts within the source region for an observationtimeof10 years.Theblackdatapoints visualiseone random Poisson realisation, drawn from the model predictions and hadronic models are fitted to the generated pseudo data sets8by minimising the ‘Cash statistic’ [58] C(ξ) =−2lnL(ξ), (2) where L(ξ) = N  i=1 P(ni|νi(ξ)) (3) an unbinned likelihood formalism, which can enhance the sensitivity. An unbinned analysis is however not yet implemented in Gammapy. 8The IC model, predicting γrays but no neutrinos, is of course fitted to the CTA data sets only. 123 112 Page 8 of 19 Eur. Phys. J. C (2024) 84 :112 denotes the total likelihood to observe the generated data, and P(ni|νi(ξ)) =νni i(ξ) ni!·exp(−νi(ξ)) (4) is the Poisson probability to measure nievents in pixel i, given a model prediction νi(ξ) that depends on the model parameters ξand is computed taking into account the IRFs of theinstruments.Severaldatasetscanbesimultaneouslyfitted by multiplying their respective likelihood values. Because the data sets are analysed with the same IRFs that were also used to create them, systematic uncertainties related to the generation of the IRFs are not taken into account here. Confidence intervals for specific model parameters can be obtained by means of a profile likelihood scan. In the scan, the parameter of interest xis consecutively fixed to values xscan around the optimum value ˆx, while the other model parameters are optimised in each step. A confidence interval can then be derived from the difference in ‘test statistic’, ΔTS(xscan)=−2lnL(xscan,ξ) L(ˆ ξ) ,(5) where ˆ ξare the parameter values for which Lis maximal. As a result of the re-optimisation of the other parameters ξ the test statistic is effectively only a function of the parameter of interest. 2.5 Constraining the hadronic contribution In this work, we derive credible intervals for the contribution of hadronic emission processes (i.e. the PD model) to the total γ-ray emission of the investigated sources. To this end, we simultaneously fit a PD model and an IC model to each of the pseudo data sets (i.e. the total γ-ray emission is given by the sum of the two models), which were generated with either the PD model or the IC model as input. The free parameters ξof this composite model are the parameters ξpdescribing the proton population and the parameters ξefor the electron distribution (cf. Eq. 1). Our parameter of interest is then the ‘hadronic fraction’ f(ξ) =Ihad(ξp) Ihad(ξp)+Ilep(ξe),(6) where Ihad and Ilep denote the integrated γ-ray flux between 100 GeV and 100 TeV for the best-fit hadronic (PD) and leptonic (IC) model, respectively. Since the hadronic fraction is not a direct parameter of the model and we cannot simply fix it to certain values, we add a penalty term Pto the Cash statistic, Ctot =C(ξ) +P(ξ, fscan), (7) where P(ξ, fscan)=Ap·(f(ξ) −fscan)2/Δ f2(8) and fscan is the value of fthat we want to probe. This penalty term allows us to fully re-optimise the model (i.e. all its direct parameters ξ) while maintaining a hadronic fraction close to fscan ±Δf. It is thus mostly a technical tool that enables us to carry out profile likelihood scans for f. We scan 21 values equally spaced between 0 and 1.9Values of Ap=0.1 and Δf=0.01wereempiricallyfoundtoleadtoastrong-enough constraint – that is, to ensure that the allowed variation in f is small compared to the spacing of the scan values – while yielding stable results. In the following, we denote with ˆ ξscan the best-fit parameter values for a given fscan, and with ˆ ξthe parameter values corresponding to the overall best fit (i.e. without a penalty term that constrains f). In the limit of sufficient statistics and for parameter values far enough from parameter boundaries, Wilk’s theorem [59] states that ΔTS follows a χ2distribution and can therefore directly be used to deduce confidence intervals [60]. However, as the parameter fis bounded between 0 and 1, we cannot invoke Wilk’s theorem here. An alternative method would be to derive the expected distribution of ΔTS by generating and fitting a large number (100) of pseudo data sets. Unfortunately, the profile likelihood scan with full reoptimisation of all model parameters is rather computingintensive, implying that this approach is also not feasible here. We therefore adopt a Bayesian approach, in which we infer a posterior probability density function (PDF) Φ( f) from the likelihood ratio L(ˆ ξscan)/L(ˆ ξ). By definition (cf. Eq. 5), L(ˆ ξscan) L(ˆ ξ) =L(fscan,ξ) L(ˆ ξ) =exp −1 2ΔTS(fscan),(9) which we use to derive the PDF as Φ( f)=c·exp −1 2ΔTS(f),(10) where f≡f(ˆ ξscan)and cis a normalisation constant that can be determined by requiring 1 0Φ( f)df=1. To obtain a smooth curve, we fit the ΔTS values obtained from the profile likelihood scan with a cubic spline function. A central Bayesian credible interval for fcan then be derived by integrating the PDF around the best-fit value up to a certain probability (e.g. 68% or 90%) – this is also known as the 9We note that the contribution of the penalty term is not included in the further evaluation of the likelihood. Furthermore, in order to take into account the small possible variations of faround fscan, the exact hadronic fractions resulting from the optimised model are used in the subsequent analysis. 123 Eur. Phys. J. C (2024) 84 :112 Page 9 of 19 112 construction of a highest posterior density interval [61]. We assume a flat prior distribution for fin this procedure. An exemplary PDF together with its 90% credible interval is shown in Fig.16 in Appendix 1. 3 Results In this section, we will introduce the different analysis scenarios we have considered (Sect. 3.1), present the results (Sect. 3.2) and discuss their implications (Sect. 3.3). 3.1 Analysis scenarios For every source, input model (leptonic or hadronic), and pseudo data set, we begin by fitting the PD and IC model to only the CTA data sets. Here we obtain the optimal parameters ξpand ξefor each model, which serve as starting parameters for the following, more complicated, composite model. In total we perform the analysis in three different scenarios, each with the goal of recovering the hadronic fraction and its uncertainty (cf. Sect. 2.5). The first scenario is a “CTA only” analysis, in which we analyse only the γ-ray data provided by CTA and ignore the KM3NeT neutrino data. This scenario is intended to illustrate to which degree CTA alone can differentiate between the two models, based on the γ-ray energy spectrum. As mentioned before, we perform a profile likelihood scan of the hadronic fraction by re-optimising the composite model (the sum of PD and IC model) at different ratios of hadronic to leptonic γ-ray flux prediction. Second, we perform a “KM3NeT only” analysis, where we include onlytheneutrinodataprovidedby KM3NeT.However,asthe primary CR energy spectrum can presently not be measured well with neutrinos alone, we constrain the parameters of the PD model to be consistent with those derived in the fit of the pure PD model to the CTA data sets, using the same concept of penalty terms as for the hadronic fraction f(cf. the previous section). Specifically, for each free parameter pof the PD model, we add a term (p−ˆp)2/Δ ˆp2, where ˆpand Δˆp are the best-fit parameter value and its uncertainty, respectively. Thus, the second scenario shows how well the leptonic and hadronic scenario can be distinguished with KM3NeT data if the primary CR energy spectrum is known to a certain degree. As the IC model is irrelevant for the neutrino flux, the hadronic fraction now corresponds to the ratio of γ-ray flux expected from the fitted proton distribution to the total γ-ray flux measured with CTA. Finally, in a third scenario we combine the γ-ray and neutrino data and perform a joint analysis of the CTA and KM3NeT data sets. As the primary CR energy spectra are now directly constrained by the CTA data, the prior terms added for the second scenario are removed again. This scenario demonstrates the benefits of the combined analysis. We note that, if the best-fit models Fig. 11 Profile likelihood scan results for Vela X and a hadronic input model ( fin =1). The displayed curves show the average of the curves obtained for the 100 generated pseudo data sets. The red and blue ‘×’ markers correspond to the analysis scenarios “CTA only” and “KM3NeT only”, respectively. The green line shows the result for the combined analysis, together with 68% and 95% quantile intervals to indicate the statistical spread. The dashed red and blue lines show the respective contributions of the CTA and KM3NeT data sets to the combined result obtained in the three scenarios are similar, a combination of the profile likelihood scans performed in the first two scenarios will lead to the same constraints on the hadronic fraction as the combined analysis in the third scenario. While this is often the case for the relatively simple models employed here, the situation can be different for more complex models, for which the combined analysis may be able to break degeneracies between model parameters. 3.2 Analysis results As an example, we show in Fig. 11 average profile likelihood scans of the hadronic fraction ffor Vela X, where the hadronic PD model has been used as input when generating the pseudo data sets. Corresponding plots for the other sources and for a purely leptonic (IC) input model can be found in Appendix B. By definition, ΔTS =0 at the minimum of the curves, which (as expected) occurs at the value of fthat corresponds to the input model (i.e. fin =0fora purelyleptonicinput modeland fin =1for apurely hadronic input model). Moving away from the minimum, larger values of ΔTS imply a stronger rejection of the corresponding hadronic fraction f. In Fig. 12, we display the 68% and 90% quantile intervals of the distribution of best-fit values of the hadronic fraction ffor all 100 pseudo data sets. Furthermore, as described in Sect. 2.5, we obtain credible intervals for ffrom the profile likelihood scans. The average expected 68% credible intervals are indicated by the black bars. We note that, for indi- 123 112 Page 16 of 19 Eur. Phys. J. C (2024) 84 :112 References 1. P. Mészáros, D.B. Fox, C. Hanna, K. Murase, Nat. Rev. Phys. 1, 585 (2019). https://doi.org/10.1038/s42254-019-0101-z 2. B.P. Abbott, R. Abbott, T.D. Abbott et al., ApJ 848, L12 (2017). https://doi.org/10.3847/2041-8213/aa91c9 3. M.G. Aartsen, M. Ackermann, J. Adams et al., Science 361, eaat1378 (2018). https://doi.org/10.1126/science.aat1378 4. S.E. Gossan, E.D. Hall, S.M. Nissanke, ApJ 926, 231 (2022). https://doi.org/10.3847/1538-4357/ac4164 5. F.Halzen,Astropart.Phys.43,155(2013).https://doi.org/10.1016/ j.astropartphys.2011.10.003 6. F. Vissani, Astropart. Phys. 26, 310 (2006). https://doi.org/10. 1016/j.astropartphys.2006.07.005 7. M.D. Kistler, J.F. Beacom, PRD 74, 063007 (2006). https://doi. org/10.1103/PhysRevD.74.063007 8. F.L. Villante, F. Vissani, PRD 78, 103007 (2008). https://doi.org/ 10.1103/PhysRevD.78.103007 9. A. Kappes, J. Hinton, C. Stegmann, F. Aharonian, ApJ 656, 870 (2007). https://doi.org/10.1086/508936 10. M.C. Gonzalez-Garcia, F. Halzen, V. Niro, Astropart. Phys. 57–58, 39 (2014). https://doi.org/10.1016/j.astropartphys.2014.04.001 11. F. Halzen, A. Kheirandish, V. Niro, Astropart. Phys. 86, 46 (2017). https://doi.org/10.1016/j.astropartphys.2016.11.004 12. S. Celli, A. Palladino, F. Vissani, Eur. Phys. J. C 77, 66 (2017). https://doi.org/10.1140/epjc/s10052-017-4635-x 13. B.S. Acharya, I. Agudo, I. Al Samarai, et al., Science with the Cherenkov Telescope Array (World Scientific Publishing, 2018). https://doi.org/10.1142/10986 14. W. Hofmann, R. Zanin, arXiv e-prints (2023). arXiv:2305.12888 15. S. Adrián-Martínez, M. Ageron, F. Aharonian et al., J. Phys. G Nucl. Part. 43, 084001 (2016) 16. S. Aiello, A. Albert, S. Alves Garre et al., Eur. Phys. J. C 82,26 (2022). https://doi.org/10.1140/epjc/s10052-021-09893-0 17. S. Adrián-Martínez, M. Ageron, F. Aharonian et al., Eur. Phys. J. C 74, 3056 (2014). https://doi.org/10.1140/epjc/s10052-014-3056-3 18. S. Adrián-Martínez, M. Ageron, F. Aharonian et al., Eur. Phys. J. C76, 54 (2016). https://doi.org/10.1140/epjc/s10052-015-3868-9 19. S. Aiello, A. Albert, M. Alshamsi et al., J. Instrum. 17, P07038 (2022). https://doi.org/10.1088/1748-0221/17/07/P07038 20. M. Ahlers, F. Halzen, Prog. Theor. Exp. Phys. 2017, 12A105 (2017). https://doi.org/10.1093/ptep/ptx021 21. M.Ahlers,F.Halzen,Prog.Part.Nucl.Phys.102,73(2018).https:// doi.org/10.1016/j.ppnp.2018.05.001 22. S. Aiello, S.E. Akrame, F. Ameli et al., Astropart. Phys. 111, 100 (2019). https://doi.org/10.1016/j.astropartphys.2019.04.002 23. L.Ambrogi,S.Celli, F.Aharonian,Astropart.Phys.100, 69(2018). https://doi.org/10.1016/j.astropartphys.2018.03.001 24. C. Deil, M. Wood, T. Hassan et al., Data formats for gamma-ray astronomy—version 0.3 (2022). https://doi.org/10.5281/zenodo. 7304668 25. C. Deil, R. Zanin, J. Lefaucheur, et al., in Proc. 35th Int. Cosmic Ray Conf., ICRC2017 (2017), p. 766. arXiv:1709.01751 26. C. Deil, A. Donath, R. Terrier et al., Zenodo (2020). https://doi. org/10.5281/zenodo.4701492 27. L. Mohrmann, A. Specovius, D. Tiziani, S. Funk, D. Malyshev, K. Nakashima, C. van Eldik, A&A 632, A72 (2019). https://doi.org/ 10.1051/0004-6361/201936452 28. A. Abramowski, F. Acero, F. Aharonian et al., A&A 548,A38 (2012). https://doi.org/10.1051/0004-6361/201219919 29. D. Horns, F. Aharonian, A. Santangelo, A.I.D. Hoffmann, C. Masterson, A&A 451, L51 (2006). https://doi.org/10.1051/0004-6361: 20065116 30. L. Zhang, X.C. Yang, ApJ 699, L153 (2009). https://doi.org/10. 1088/0004-637X/699/2/L153 31. H. Abdalla, A. Abramowski, F. Aharonian et al., A&A 612,A6 (2018). https://doi.org/10.1051/0004-6361/201629790 32. Y. Fukui, H. Sano, Y. Yamane, T. Hayakawa, T. Inoue, K. Tachihara, G. Rowell, S. Einecke, ApJ 915, 84 (2021). https://doi.org/ 10.3847/1538-4357/abff4a 33. P. Cristofari, V. Niro, S. Gabici, MNRAS 508, 2204 (2021). https:// doi.org/10.1093/mnras/stab2380 34. J.S. Clark, I. Negueruela, P. Crowther, S.P. Goodwin, A&A 434, 949 (2005). https://doi.org/10.1051/0004-6361:20042413 35. A. Abramowski, F. Acero, F. Aharonian et al., A&A 537, A114 (2012). https://doi.org/10.1051/0004-6361/201117928 36. F. Aharonian, H. Ashkar, M. Backes et al., A&A 666, A124 (2022). https://doi.org/10.1051/0004-6361/202244323 37. F. Aharonian, R. Yang, E. de Oña Wilhelmi, Nat. Astron. 3, 561 (2019). https://doi.org/10.1038/s41550-019-0724-0 38. F. Aharonian, A.G. Akhperjanian, G. Anton et al., A&A 499, 723 (2009). https://doi.org/10.1051/0004-6361/200811357 39. A.U. Abeysekara, A. Albert, R. Alfaro et al., PRL 124, 021102 (2020). https://doi.org/10.1103/PhysRevLett.124.021102 40. Z. Cao, F.A. Aharonian, Q. An et al., Nature 594, 33 (2021). https:// doi.org/10.1038/s41586-021-03498-z 41. R. Muller, A. Heijboer, A. Garcia Soto, et al., in Proc. 37th Int. Cosmic Ray Conf., ICRC2021 (2021), p. 1077. https://doi.org/10. 22323/1.395.1077 42. CTAO & CTAC, CTAO Instrument Response Functions—prod5 version v0.1 (2021). https://doi.org/10.5281/zenodo.5499840 43. T. Sudoh, J.F. Beacom, PRD 108, 043016 (2023). https://doi.org/ 10.1103/PhysRevD.108.043016 44. S. Aiello, A. Albert, S. AlvesGarre et al., Comput. Phys. Commun. 256, 107477 (2020). https://doi.org/10.1016/j.cpc.2020.107477 45. C. Andreopoulos, A. Bell, D. Bhattacharya, F. Cavanna, J. Dobson, S. Dytman, H. Gallagher, P. Guzowski, R. Hatcher, P. Kehayias, A. Meregaglia, D. Naples, G. Pearce, A. Rubbia, M. Whalley, T. Yang, Nucl. Instrum. Methods Phys. Res. A 614, 87 (2010). https://doi. org/10.1016/j.nima.2009.12.009 46. M. Honda, T. Kajita, K. Kasahara, S. Midorikawa, T. Sanuki, PRD 75, 043006 (2007). https://doi.org/10.1103/PhysRevD.75.043006 47. R.Enberg,M.H. Reno, I. Sarcevic, PRD 78,043005(2008). https:// doi.org/10.1103/PhysRevD.78.043005 48. M.G. Aartsen, R. Abbasi, M. Ackermann et al., PRD 89, 062007 (2014). https://doi.org/10.1103/PhysRevD.89.062007 49. T.K. Gaisser, Astropart. Phys. 35, 801 (2012). https://doi.org/10. 1016/j.astropartphys.2012.02.010 50. C. Mascaretti, P. Blasi, C. Evoli, Astropart. Phys. 114, 22 (2020). https://doi.org/10.1016/j.astropartphys.2019.06.002 51. Y. Becherini, A. Margiotta, M. Sioli, M. Spurio, Astropart. Phys. 25, 1 (2006). https://doi.org/10.1016/j.astropartphys.2005.10.005 52. G.Carminati,M.Bazzotti,A.Margiotta,M. Spurio, Comput. Phys. Commun. 179, 915 (2008). https://doi.org/10.1016/j.cpc.2008.07. 014 53. A. Albert, M. André, M. Anghinolfi et al., JCAP 01, 064 (2021). https://doi.org/10.1088/1475-7516/2021/01/064 54. R. Abbasi, M. Ackermann, J. Adams et al., ApJ 928, 50 (2022). https://doi.org/10.3847/1538-4357/ac4d29 55. V. Zabalza, in Proc. 34th Int. Cosmic Ray Conf., ICRC2015 (2015), p. 922. arXiv:1509.03319 56. S.R. Kelner, F.A. Aharonian, V.V. Bugayov, PRD 74, 034018 (2006) 57. E.Kafexhiu,F.Aharonian,A.M.Taylor,G.S.Vila,PRD 90,123014 (2014). https://doi.org/10.1103/PhysRevD.90.123014 58. W. Cash, ApJ 228, 939 (1979). https://doi.org/10.1086/156922 59. S.S. Wilks, Ann. Math. Stat. 9, 60 (1938). https://doi.org/10.1214/ aoms/1177732360 60. R.L.Workmanetal., Prog.Theor. Exp.Phys.2022,083C01 (2022). https://doi.org/10.1093/ptep/ptac097 123 Eur. Phys. J. C (2024) 84 :112 Page 17 of 19 112 61. J.O. Berger, Statistical Decision Theory: Foundations, Concepts, and Methods (Springer, New York, 1980). https://doi.org/10.1007/ 978-1-4757-1727-3 62. L. Tibaldo, R. Zanin, G. Faggioli, J. Ballet, M.H. Grondin, J.A. Hinton, M. Lemoine-Goumard, A&A 617, A78 (2018). https:// doi.org/10.1051/0004-6361/201833356 63. A.U. Abeysekara, A. Albert, R. Alfaro et al., ApJ 843, 39 (2017). https://doi.org/10.3847/1538-4357/aa7555 64. A. Albert, R. Alfaro, J.C. Arteaga-Velázquez et al., A&A 667,A36 (2022). https://doi.org/10.1051/0004-6361/202243527 65. R.Abbasi,M.Ackermann,J.Adamsetal.,Science 378,538(2022). https://doi.org/10.1126/science.abg3395 66. T. Unbehaun, L. Mohrmann, M. Smirnov, J. Schnabel, T. Gál, Prospects for combined galactic source searches with CTA and KM3NeT (2023). https://doi.org/10.5281/zenodo.8298464 67. T.P. Robitaille, E.J. Tollerud, P. Greenfield et al., A&A 558,A33 (2013). https://doi.org/10.1051/0004-6361/201322068 68. A.M. Price-Whelan, B.M. Sip˝ocz, H.M. Günther, et al., AJ 156, 123 (2018). https://doi.org/10.3847/1538-3881/aabc4f 69. J.D. Hunter, Comput. Sci. Eng. 9, 90 (2007). https://doi.org/10. 1109/MCSE.2007.55 70. N.P. Chue Hong, D.S. Katz, M. Barker, et al., FAIR Principles for Research Software (FAIR4RS Principles) (2022). https://doi.org/ 10.15497/RDA00068 Members of the CTA Consortium T. Unbehaun1,*, L. Mohrmann2,*, S. Funk1 KM3NeT Collaboration* S. Aiello3, A. Albert4,58, S. Alves Garre5,Z.Aly 6, A. Ambrosone8,7, F. Ameli9, M. Andre10, E. Androutsou11, M. Anghinolfi12, M. Anguita13, L. Aphecetche14,M.Ardid 15,S.Ardid 15, H. Atmani16, J. Aublin17, C. Bagatelas11, L. Bailly-Salins18, B. Baret17, Z. Barda˘cová40,S.BasegmezduPree 19, Y. Becherini17, M. Bendahman16,17, F. Benfenati20,21, M. Benhassi7,22, D. M. Benoit23, E. Berbee19, V. Bertin6, S. Biagi24, M. Boettcher25, M. Bou Cabo26, J. Boumaaza16, M. Bouta27, M. Bouwhuis19, C. Bozza7,28, R. M. Bozza7,8, H. Brânza¸s29, F. Bretaudeau14, R. Bruijn19,30, J. Brunner6, R. Bruno3,E.Buis 19,31, R. Buompane7,22,J.Busto 6,B.Caiffi 12,D.Calvo 5, S. Campion9,32, A. Capone9,32, F. Carenini20,21, V. Carretero5, T. Cartraud17, P. Castaldi20,33, V. Cecchini5, S. Celli9,32,L.Cerisy 6, M. Chabab34, M. Chadolias1, A. Chen35, S. Cherubini24,36, T. Chiarusi20, M. Circella37, R. Cocimano24,J.A.B.Coelho 17, A. Coleiro17, R. Coniglione24,P.Coyle 6, A. Creusot17,A.Cruz 38, G. Cuttone24, R. Dallier14,Y.Darras 1, A. De Benedittis7, B. De Martino6, V. Decoene14,R.DelBurgo 7, L.S.DiMauro 24,I.DiPalma 9,32, A.F.Díaz 13, D. Diego-Tortosa24, C. Distefano24,A.Domi 19,30, C. Donzaud17, D. Dornic6,M.Dörr 39, E. Drakopoulou11, D. Drouhin4,58,R.Dvornický 40, T. Eberl1, E. Eckerová55, A. Eddymaoui16, T. van Eeden19,M.Eff 1, D. van Eijk19, I. El Bojaddaini27, S. El Hedri17, A. Enzenhöfer6,G. Ferrara24,36,M. D. Filipovi´c41,F. Filippini20,21,L. A. Fusco28,J. Gabriel42,T.Gal1,J. García Méndez15, A. Garcia Soto5, C. Gatius Oliver19, N. Geißelbrecht1, H. Ghaddari27, L. Gialanella7,22, B.K.Gibson 23, E. Giorgio24, A. Girardi9, I. Goos17, S. R. Gozzini5, R. Gracia1,K.Graf 1, D. Guderian59, C. Guidi12,43, B. Guillon18, M. Gutiérrez44, H. van Haren45, A. Heijboer19, A. Hekalo39, L. Hennig1, J. J. Hernández-Rey5, F. Huang6, W. Idrissi Ibnsalih7, G. Illuminati20,21,C.W.James 38, M. de Jong19,46, P. de Jong19,30,B.J.Jung 19, P. Kalaczy´nski47, O. Kalekin1,U.F.Katz 1, N. R. Khan Chowdhury5, A. Khatun40, G. Kistauri48,49, F. van der Knaap31, A. Kouchner17,50, V. Kulikovskiy12, R. Kvatadze49, M. Labalme18, R. Lahmann1, G. Larosa24,C.Lastoria 6, A. Lazo5,S.LeStum 6, G. Lehaut18, E. Leonora3, N. Lessing5,G.Levi 20,21, M. Lindsey Clark17, F. Longhitano3, J. Majumdar19, L. Malerba12,J.Ma´nczak5, A. Manfreda7, M. Marconi12,43, A. Margiotta20,21, A. Marinelli7,8,C.Markou 11, L. Martin14, F. Marzaioli7,22, M. Mastrodicasa9,32, S. Mastroianni7, S. Miccichè24, G. Miele7,8, P. Migliozzi7, E. Migneco24, P. Mijakowski47, M. L. Mitsou7, C. M. Mollo7, L. Morales-Gallegos7, C. Morley-Wong38, A. Mosbrugger1, A. Moussa27, I. Mozun Mateo51,52, R. Muller19, M. R. Musone7,22, M. Musumeci24, L. Nauta19,S.Navas 44, A. Nayerhoda37, C. A. Nicolau9, B. Nkosi35, B. Ó Fearraigh19,30, V. Oliviero7,8, A. Orlando24, E. Oukacha17, J. Palacios González5, G. Papalashvili48, E. J. Pastor Gomez5, A.M.P˘aun29, G.E.P˘av˘ala¸s29, S. Peña Martínez17, M. Perrin-Terrin6, J. Perronnel18, V. Pestel52, R. Pestes17, P. Piattelli24, C. Poirè28, V. Popa29, T. Pradier4, S. Pulvirenti24, G. Quéméner18,C.Quiroz 15, U. Rahaman5, N. Randazzo3, S. Razzaque53,I.C.Rea 7, D. Real5, S. Reck1, G. Riccobene24, J. Robinson25, A. Romanov12,43, L. Roscilli7, A. Saina5, F. Salesa Greus5, D. F. E. Samtleben19,46, A. Sánchez Losa5,37, M. Sanguineti12,43, C. Santonastaso7,54, D. Santonocito24, P. Sapienza24, J. Schnabel1, M. F. Schneider1, J. Schumann1, H. M. Schutte25, J. Seneca19, N. Sennan27, B. Setter1,I. Sgura37,R. Shanidze48,Y. Shitov55,F.Šimkovic40,A. Simonelli7,A. Sinopoulou3,M.V.Smirnov1,B. Spisso7, M. Spurio20,21, D. Stavropoulos11, I. Štekl55, M. Taiuti12,43, Y. Tayalati16, H. Tedjditi12, H. Thiersen25, I. Tosta e Melo3,36, B. Trocme17, S. Tsagkli11, V. Tsourapis11, E. Tzamariudaki11, A. Vacheret18, V. Valsecchi24, V. Van Elewyck17,50, G. Vannoye6, G. Vasileiadis56, F. Vazquez de Sola19, C. Verilhac17, A. Veutro9,32, S. Viola24,D.Vivolo 7,22, H. Warnhofer1, J. Wilms57,E.deWolf 19,30, T. Yousfi27, G. Zarpapis11, S. Zavatarelli12, A. Zegarelli9,32, D. Zito24, J. D. Zornoza5, J. Zúñiga5, N. Zywucka25 123 112 Page 18 of 19 Eur. Phys. J. C (2024) 84 :112 1Erlangen Centre for Astroparticle Physics, University of Erlangen-Nuremberg, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany 2Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, Heidelberg 69117, Germany 3INFN, Sezione di Catania, Via Santa Sofia 64, 95123 Catania, Italy 4CNRS, IPHC UMR 7178, Université de Strasbourg, 67000 Strasbourg, France 5IFIC-Instituto de Física Corpuscular (CSIC-Universitat de València), c/Catedrático José Beltrán 2, 46980 Paterna, Valencia, Spain 6CNRS/IN2P3, CPPM, Aix Marseille Univ, Marseille, France 7INFN, Sezione di Napoli, Complesso Universitario di Monte S. Angelo, Via Cintia ed. G, 80126 Naples, Italy 8Università di Napoli “Federico II”, Dip. Scienze Fisiche “E. Pancini”, Complesso Universitario di Monte S. Angelo, Via Cintia ed. G, 80126 Naples, Italy 9INFN, Sezione di Roma, Piazzale Aldo Moro 2, 00185 Rome, Italy 10 Laboratori d’Aplicacions Bioacústiques, Centre Tecnològic de Vilanova i la Geltrú, Universitat Politècnica de Catalunya, Avda. Rambla Exposició s/n, Vilanova i la Geltrú 08800, Spain 11 Institute of Nuclear and Particle Physics, NCSR Demokritos, Ag. Paraskevi Attikis, 15310 Athens, Greece 12 INFN, Sezione di Genova, Via Dodecaneso 33, 16146 Genoa, Italy 13 Department of Computer Architecture and Technology/CITIC, University of Granada, 18071 Granada, Spain 14 Subatech, IMT Atlantique, IN2P3-CNRS, Université de Nantes, 4 rue Alfred Kastler-La Chantrerie, BP 20722, 44307 Nantes, France 15 Instituto de Investigación para la Gestión Integrada de las Zonas Costeras, Universitat Politècnica de València, C/ Paranimf 1, Gandia 46730, Spain 16 Faculty of Sciences, University Mohammed V in Rabat, 4 av. Ibn Battouta, B.P. 1014, R.P. 10000 Rabat, Morocco 17 CNRS, Astroparticule et Cosmologie, Université Paris Cité, 75013 Paris, France 18 LPC CAEN, Normandie Univ, ENSICAEN, UNICAEN, CNRS/IN2P3, 6 Boulevard Maréchal Juin, 14050 Caen, France 19 Nikhef, National Institute for Subatomic Physics, PO Box 41882, 1009 DB Amsterdam, The Netherlands 20 INFN, Sezione di Bologna, v.le C. Berti-Pichat 6/2, 40127 Bologna, Italy 21 Dipartimento di Fisica e Astronomia, Università di Bologna, v.le C. Berti-Pichat 6/2, Bologna 40127, Italy 22 Dipartimento di Matematica e Fisica, Università degli Studi della Campania “Luigi Vanvitelli”, Viale Lincoln 5, 81100 Caserta, Italy 23 E.A. Milne Centre for Astrophysics, University of Hull, Hull HU6 7RX, UK 24 INFN, Laboratori Nazionali del Sud, Via S. Sofia 62, Catania 95123, Italy 25 Centre for Space Research, North-West University, Private Bag X6001, Potchefstroom 2520, South Africa 26 Instituto Español de Oceanografía, Unidad Mixta IEO-UPV, C/ Paranimf 1, 46730 Gandia, Spain 27 Faculty of Sciences, University Mohammed I, BV Mohammed VI, B.P. 717, R.P. 60000 Oujda, Morocco 28 Dipartimento di Fisica, Università di Salerno e INFN Gruppo Collegato di Salerno, Via Giovanni Paolo II 132, 84084 Fisciano, Italy 29 ISS, Atomistilor 409, 077125 M˘agurele, Romania 30 Institute of Physics/IHEF, University of Amsterdam, PO Box 94216, 1090 GE Amsterdam, The Netherlands 31 TNO, Technical Sciences, PO Box 155, 2600 AD Delft, The Netherlands 32 Dipartimento di Fisica, Università La Sapienza, Piazzale Aldo Moro 2, 00185 Rome, Italy 33 Dipartimento di Ingegneria dell’Energia Elettrica e dell’Informazione ”Guglielmo Marconi”, Università di Bologna, Via dell’Università 50, 47521 Cesena, Italy 34 Physics Department, Faculty of Science Semlalia, Cadi Ayyad University, Av. My Abdellah, P.O.B. 2390, 40000 Marrakech, Morocco 35 School of Physics, University of the Witwatersrand, Private Bag 3, 2050 Johannesburg, Wits, South Africa 36 Dipartimento di Fisica e Astronomia “Ettore Majorana”, Università di Catania, Via Santa Sofia 64, Catania 95123, Italy 37 INFN, Sezione di Bari, Via Orabona 4, 70125 Bari, Italy 38 International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia 39 University Würzburg, Emil-Fischer-Straße 31, 97074 Würzburg, Germany 40 Comenius University in Bratislava, Department of Nuclear Physics and Biophysics, Mlynska dolina F1, Bratislava 842 48, Slovak Republic 123 Eur. Phys. J. C (2024) 84 :112 Page 19 of 19 112 41 School of Computing, Engineering and Mathematics, Western Sydney University, Locked Bag 1797, Penrith, NSW 2751, Australia 42 IN2P3, LPC, Campus des Cézeaux 24, Avenue des Landais, BP 80026, 63171 Aubière Cedex, France 43 Università di Genova, Via Dodecaneso 33, 16146 Genoa, Italy 44 Dpto. de Física Teórica y del Cosmos & C.A.F.P.E., University of Granada, 18071 Granada, Spain 45 NIOZ (Royal Netherlands Institute for Sea Research), PO Box 59, Den Burg, 1790 AB Texel, The Netherlands 46 Leiden Institute of Physics, Leiden University, PO Box 9504, 2300 RA Leiden, The Netherlands 47 National Centre for Nuclear Research, 02-093, Warsaw, Poland 48 Department of Physics, Tbilisi State University, 3, Chavchavadze Ave., 0179 Tbilisi, Georgia 49 Institute of Physics, The University of Georgia, Kostava str. 77, 0171 Tbilisi, Georgia 50 Institut Universitaire de France, 1 Rue Descartes, 75005 Paris, France 51 IN2P3, 3, Rue Michel-Ange, 75794 Paris 16, France 52 LPC, Campus des Cézeaux 24, Avenue des Landais, BP 80026, 63171 Aubière Cedex, France 53 Department Physics, University of Johannesburg, PO Box 524, 2006 Auckland Park, South Africa 54 Laboratorio CIRCE, Dip. Di Matematica e Fisica, Università degli Studi della Campania “Luigi Vanvitelli”, CAPACITY, Viale Carlo III di Borbone 153, 81020 San Nicola La Strada, Italy 55 Czech Technical University in Prague, Institute of Experimental and Applied Physics, Husova 240/5, Prague 110 00, Czech Republic 56 Laboratoire Univers et Particules de Montpellier, Place Eugène Bataillon-CC 72, 34095 Montpellier Cédex 05, France 57 Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Remeis Sternwarte, Sternwartstraße 7, 96049 Bamberg, Germany 58 Université de Haute Alsace, Rue des Frères Lumière, 68093 Mulhouse Cedex, France 59 Institut für Kernphysik, University of Münster, Wilhelm-Klemm-Str. 9, 48149 Münster, Germany 123