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

T Unbehaun, L Mohrmann, S Funk, S Aiello, A Albert, S. Alves Garre, Z Aly, A Ambrosone, F Ameli, M Andre, et al.

► To cite this version:

T Unbehaun, L Mohrmann, S Funk, S Aiello, A Albert, et al. Prospects for combined analyses of hadronic emission from γ-ray sources in the Milky Way with CTA and KM3NeT. European Physical Journal C: Particles and Fields, 2024, 84 (2), pp.112. 10.1140/epjc/s10052-023-12279-z. hal-04212497

HAL Id: hal-04212497

https://hal.science/hal-04212497v1

Submitted on 16 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 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** The Cherenkov Telescope Array and the KM3NeT 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 analysis of γ-ray data, to also process data from neutrino 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 200 h 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

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].

Astrophysical objects in the Milky Way are not expected to produce gravitational waves detectable by current-generation instruments (but by next-generation detectors, see [4]). Galactic CRs provide important energetic constraints on their source population, but cannot be used to directly study Galactic objects because they are deflected by magnetic fields. Hence, one needs to resort to photons and neutrinos to employ 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 ‘hadronic accelerators’, where accelerated CR nuclei interact 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 of photons by CR electrons leads 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 hadronic accelerators [5]. Nevertheless, observations of VHE γ-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,

<footnote>a e-mail: tim.unbehaun@fau.de; lars.mohrmann@mpi-hd.mpg.de; km3net-pc@km3net.de</footnote>
Spain, and Paranal, Chile, covering the Northern and South-
ern 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
measuring the faint flash of Cherenkov light that is emitted by
secondary particles in the air shower that is launched when
the primary γ 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 Mediter-
ranean, 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 proper-
ties measuring atmospheric neutrinos [16]. Here, we only
consider the ‘Astroparticle Research with Cosmics in the
Abyss’ (ARCA) detector, which targets the detection of high-
energy astrophysical neutrinos with TeV-PeV energies. Here-
after, 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 com-
prising 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 dis-
tances of up to several km in the water allow a precise recon-
struction of the direction of the incoming neutrino. Compared
to the largest existing neutrino telescope, the IceCube Neu-
trino 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 angu-
lar resolution [15]. Its location in the Northern hemisphere
implies that neutrinos potentially emitted by many Galac-
tic γ-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 been investigated in [23]. In this paper, we demonstrate how a 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 emis-

sion 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 format<fnref n="1" /> used for the publicly available
CTA IRFs. Then we have employed the GAMMAPY<fnref n="2" /> package
(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 develop-
ment 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].

We apply the analysis to a selection of prototypical sources
that are promising candidates for the emission of 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], asso-
  ciated 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] as well 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 promis-
ing 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

<footnote n="1">¹ See https://gamma-astro-data-formats.readthedocs.io and [24].</footnote>

<footnote n="2">² https://gammapy.org.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Source visibility with KM3NeT. Shown is the fraction of time that each source is visible under a zenith angle $\theta$, 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</figcaption></figure>

that furthermore exhibit favourable locations for the observation with CTA-South and KM3NeT. Figure 1 shows the visibility of all sources for KM3NeT, as a function of the local zenith angle.<fnref n="3" /> For 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 point-source 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, we introduce 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 analysis (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.

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 the detector. In our case, the relevant IRFs comprise the effective 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 $\gamma$-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' IRFs<fnref n="4" /> for 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-neutrinos<fnref n="5" /> that 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 appear as characteristic, track-like events in the detector. 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 $^{40}$K as well as bioluminescence are relevant sources of noise. Due to the design of the optical modules, which contain multi-

<footnote n="3">³ The zenith angle $\theta$ refers to the location of the source. For $\theta = 0°$ the source is above the detector and produces vertically down-going neutrinos, while for $\theta = 180°$ the source is located on the opposite side of the Earth and produces vertically up-going neutrinos.</footnote>

<footnote n="4">⁴ See https://www.cta-observatory.org/science/ctao-performance.</footnote>

<footnote n="5">⁵ Hereafter, we will use the term 'neutrinos' to refer to both neutrinos and anti-neutrinos, unless explicitly stated otherwise.</footnote>
<table><caption>**Table 1** Galactic γ-ray sources investigated in this work</caption><thead><tr><td>Designation</td><td>Type</td><td>Spatial model</td><td>r (deg)</td><td>Declination (deg)</td><td>Distance (kpc)</td><td>References</td></tr></thead><tbody><tr><td>Vela X</td><td>PWN</td><td>disk</td><td>0.8</td><td>-45.19</td><td>0.29</td><td>[28]</td></tr><tr><td>RX J1713.7-3946</td><td>SNR</td><td>disk</td><td>0.6</td><td>-39.69</td><td>1</td><td>[31]</td></tr><tr><td>Westerlund 1</td><td>SC</td><td>disk</td><td>1.1</td><td>-45.85</td><td>3.9</td><td>[35]</td></tr><tr><td>eHWC J1907+063</td><td>UNID</td><td>Gaussian</td><td>0.67</td><td>+06.18</td><td>2.37</td><td>[39]</td></tr></tbody></table>

‘Spatial model’ specifies which type of spatial model is used in the analysis (cf. Sect. 2.3). *r* denotes 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, because their parent particles have a non-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 atmospheric neutrinos, integrated over relevant zenith angles, 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.

All simulated events are reconstructed using a track reconstruction 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 $\theta_{reco} > 80^{\circ}$. Thus, only very few atmospheric muon events remain in the final sample, almost all concentrated close to the horizon region ($80^{\circ} < \theta_{reco} < 90^{\circ}$), 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 normalization 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(\theta)$. 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 3 shows 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
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Background event rates in KM3NeT as a function of reconstructed energy $E_{reco}$. The atmospheric neutrino rates are integrated over all zenith angles in the analysis region (i.e. $\theta_{reco} > 80^{\circ}$). The atmospheric muon rate is shown for the zenith angle bin $80^{\circ}-90^{\circ}$ only, since it is completely negligible for larger angles. The black line displays the smoothed curve used in the analysis</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Comparison of effective areas. The CTA effective area is shown for a zenith angle of $\theta = 20^{\circ}$ and an offset angle from the pointing direction of $\vartheta = 1^{\circ}$. For KM3NeT, average effective areas for the full analysis region as well as for different sub-ranges in zenith angle are shown</figcaption></figure>

than that of CTA for $\gamma$ 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).

<figure><img src="image_5.png" /><figcaption>**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 $\gamma$-ray/neutrino energy. Possible differences to previous publications may arise from the finite binning of the PSF applied here</figcaption></figure>

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^{\circ}$ above ~10 TeV. The PSF strongly affects the sensitivity to point-like 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 5 provides 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
<figure><img src="image_6.png" /><figcaption>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</figcaption></figure>

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 been fitted to published γ-ray spectra of the sources (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]<fnref n="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,

$$ \Phi(E) = A \cdot \left(\frac{E}{E_0}\right)^{-\Gamma} \exp\left[-\left(\frac{E}{E_{\text{cut}}}\right)^{\beta}\right], \quad (1) $$

where A denotes the amplitude, $E_0$ the reference energy, $\Gamma$ the spectral index, $E_{\text{cut}}$ the cut-off energy, and $\beta$ 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 $\chi^2$ fit, keeping the reference energy and cut-off strength fixed (at $E_0 = 10$ TeV and $\beta = 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 scenarios based on γ-ray data alone. We note that the addition

<figure><img src="image_7.png" /><figcaption>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</figcaption></figure>

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 pseudo data sets for each source and each instrument, both for 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 200 h, 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. 7 and 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 sets for each source, depending on the visibility. The exposure 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

<footnote n="6">⁶ We are aware of more recent parametrisations such as that provided by [57], which is also used in the PionDecay 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.</footnote>
<figure><img src="image_8.png" /><figcaption>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 Poisson-randomised 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. 8 have been extracted</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>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</figcaption></figure>

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. 9 and 10, respectively.

## 2.4 Likelihood analysis

The analysis is performed using the binned likelihood formalism implemented in the GAMMAPY package.<fnref n="7" /> Leptonic

<figure><img src="image_10.png" /><figcaption>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)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>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 observation time of 10 years. The black data points visualise one random Poisson realisation, drawn from the model predictions</figcaption></figure>

and hadronic models are fitted to the generated pseudo data sets<fnref n="8" /> by minimising the ‘Cash statistic’ [58]

$$ \mathcal{C}(\xi) = -2 \ln \mathcal{L}(\xi), \qquad (2) $$

where

$$ \mathcal{L}(\xi) = \prod_{i=1}^{N} \mathcal{P}(n_i | v_i(\xi)) \qquad (3) $$

<footnote n="7">⁷ We note that analyses carried out with the analysis tools normally employed by the KM3NeT Collaboration are typically performed using</footnote>

<footnote>an unbinned likelihood formalism, which can enhance the sensitivity. An unbinned analysis is however not yet implemented in GAMMAPY.</footnote>

<footnote n="8">⁸ The IC model, predicting γ rays but no neutrinos, is of course fitted to the CTA data sets only.</footnote>
denotes the total likelihood to observe the generated data, and

$$ \mathcal{P}(n_i|v_i(\xi)) = \frac{v_i^{n_i}(\xi)}{n_i!} \cdot \exp(-v_i(\xi)) \quad (4) $$

is the Poisson probability to measure $n_i$ events in pixel $i$, given a model prediction $v_i(\xi)$ that depends on the model parameters $\xi$ and is computed taking into account the IRFs of the instruments. Several data sets can be simultaneously fitted 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 $x$ is consecutively fixed to values $x_{scan}$ around the optimum value $\hat{x}$, while the other model parameters are optimised in each step. A confidence interval can then be derived from the difference in ‘test statistic’,

$$ \Delta \text{TS}(x_{\text{scan}}) = -2 \ln \left( \frac{\mathcal{L}(x_{\text{scan}}, \xi)}{\mathcal{L}(\hat{\xi})} \right), \quad (5) $$

where $\hat{\xi}$ are the parameter values for which $\mathcal{L}$ is maximal. As a result of the re-optimisation of the other parameters $\xi$ 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 $\gamma$-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 $\gamma$-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 $\xi$ of this composite model are the parameters $\xi_p$ describing the proton population and the parameters $\xi_e$ for the electron distribution (cf. Eq. 1). Our parameter of interest is then the ‘hadronic fraction’

$$ f(\xi) = \frac{I_{\text{had}}(\xi_p)}{I_{\text{had}}(\xi_p) + I_{\text{lep}}(\xi_e)}, \quad (6) $$

where $I_{\text{had}}$ and $I_{\text{lep}}$ denote the integrated $\gamma$-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 $P$ to the Cash statistic,

$$ \mathcal{C}_{\text{tot}} = \mathcal{C}(\xi) + P(\xi, f_{\text{scan}}), \quad (7) $$

where

$$ P(\xi, f_{\text{scan}}) = A_p \cdot (f(\xi) - f_{\text{scan}})^2 / \Delta f^2 \quad (8) $$

and $f_{\text{scan}}$ is the value of $f$ that we want to probe. This penalty term allows us to fully re-optimise the model (i.e. all its direct parameters $\xi$) while maintaining a hadronic fraction close to $f_{\text{scan}} \pm \Delta 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.<fnref n="9" /> Values of $A_p = 0.1$ and $\Delta f = 0.01$ were empirically found to lead to a strong-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 $\hat{\xi}_{\text{scan}}$ the best-fit parameter values for a given $f_{\text{scan}}$, and with $\hat{\xi}$ 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 $\Delta \text{TS}$ follows a $\chi^2$ distribution and can therefore directly be used to deduce confidence intervals [60]. However, as the parameter $f$ is bounded between 0 and 1, we cannot invoke Wilk’s theorem here. An alternative method would be to derive the expected distribution of $\Delta \text{TS}$ by generating and fitting a large number ($\gg 100$) of pseudo data sets. Unfortunately, the profile likelihood scan with full re-optimisation of all model parameters is rather computing-intensive, 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) $\Phi(f)$ from the likelihood ratio $\mathcal{L}(\hat{\xi}_{\text{scan}})/\mathcal{L}(\hat{\xi})$. By definition (cf. Eq. 5),

$$ \frac{\mathcal{L}(\hat{\xi}_{\text{scan}})}{\mathcal{L}(\hat{\xi})} = \frac{\mathcal{L}(f_{\text{scan}}, \xi)}{\mathcal{L}(\hat{\xi})} = \exp \left( -\frac{1}{2} \Delta \text{TS}(f_{\text{scan}}) \right), \quad (9) $$

which we use to derive the PDF as

$$ \Phi(f) = c \cdot \exp \left( -\frac{1}{2} \Delta \text{TS}(f) \right), \quad (10) $$

where $f \equiv f(\hat{\xi}_{\text{scan}})$ and $c$ is a normalisation constant that can be determined by requiring $\int_0^1 \Phi(f) df = 1$. To obtain a smooth curve, we fit the $\Delta \text{TS}$ values obtained from the profile likelihood scan with a cubic spline function. A central Bayesian credible interval for $f$ can 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

<footnote n="9">⁹ We 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 $f$ around $f_{\text{scan}}$, the exact hadronic fractions resulting from the optimised model are used in the subsequent analysis.</footnote>
construction of a highest posterior density interval [61]. We assume a flat prior distribution for *f* in 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 $\xi_p$ and $\xi_e$ for 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 $\gamma$-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 $\gamma$-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 $\gamma$-ray flux prediction. Second, we perform a “KM3NeT only” analysis, where we include only the neutrino data provided by KM3NeT. However, as the 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 *p* of the PD model, we add a term $(p - \hat{p})^2/\Delta\hat{p}^2$, where $\hat{p}$ and $\Delta\hat{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 $\gamma$-ray flux expected from the fitted proton distribution to the total $\gamma$-ray flux measured with CTA. Finally, in a third scenario we combine the $\gamma$-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

<figure><img src="image_12.png" /><figcaption>**Fig. 11** Profile likelihood scan results for Vela X and a hadronic input model ($f_{in} = 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</figcaption></figure>

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 *f* for 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, $\Delta\text{TS} = 0$ at the minimum of the curves, which (as expected) occurs at the value of *f* that corresponds to the input model (i.e. $f_{in} = 0$ for a purely leptonic input model and $f_{in} = 1$ for a purely hadronic input model). Moving away from the minimum, larger values of $\Delta\text{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 *f* for all 100 pseudo data sets. Furthermore, as described in Sect. 2.5, we obtain credible intervals for *f* from the profile likelihood scans. The average expected 68% credible intervals are indicated by the black bars. We note that, for indi-
<figure><img src="image_13.png" /><figcaption>**Fig. 12** Summary of results for all sources and input models. The blue and red horizontal bars show the 68% and 90% quantile intervals of the distribution of best-fit values for the leptonic ($f_{in} = 0$) and hadronic ($f_{in} = 1$) input models, respectively. The black bars indicate the average position and sizes of the 68% credible intervals derived for each input model</figcaption></figure>

vidual pseudo data sets, it is possible that the 68% credible interval does not contain the true input value of $f$. Hence, this is also the case for the averaged interval: this is an artefact caused by the fact that the input values lie at the boundaries of the allowed values for $f$.

## 3.3 Discussion

It is evident that in essentially all cases, the “CTA only” and “KM3NeT only” scenarios yield results that correspond exactly to the contributions of the two instruments in the combined analysis. That is, in Figs. 11, 14, and 15, the red and blue crosses fall on top of the red and blue dashed lines, respectively. This means that consistent source models are fitted in all three scenarios and implies that, in principle, a sensitivity identical to that of the combined analysis can be achieved by combining the results of the profile likelihood scans of the single-instrument analyses. We note, however, that this only holds because we use the same source models in all analysis scenarios, and because we incorporate the CTA constraints on the primary CR spectrum in the “KM3NeT only” analysis. In practice, ensuring this consistency between two separate analyses may not always be possible – in particular when more sophisticated source models are considered, which could, for example, feature different spatial models for the leptonic and hadronic cases. The combined analysis, on the other hand, naturally ensures that the $\gamma$-ray and neutrino data are described with the same physical models. Furthermore, a combined analysis enables a consistent incorporation of systematic uncertainties: for example, it would be possible to add to the analysis a nuisance parameter that modifies the relative energy scale between the two instruments – something that would be much more difficult to take into account

when combining the results of two single-instrument analyses.

Investigating the profile likelihood scans for the different sources, we note that the curves obtained from the KM3NeT data sets are usually very similar for the leptonic and hadronic input model, except being flipped horizontally, and often yield stronger constraints on the hadronic fraction $f$ compared to the CTA curves. This is not unexpected, as the predicted neutrino fluxes differ fundamentally between the leptonic model (no neutrinos) and the hadronic model (neutrino flux approximately equal to $\gamma$-ray flux), whereas the predicted $\gamma$-ray fluxes can look very similar (cf. Figs. 6, 13). This is especially true here, as we use the same spatial models for the leptonic and hadronic scenario, implying that for CTA any separation power between the models originates from differences in the predicted $\gamma$-ray spectra. For more realistic models that also differ in their morphology, the CTA data will lead to better constraints than achieved here – the presented curves should therefore not be regarded as a generally valid estimate of the CTA sensitivity to differentiate between leptonic and hadronic models. The KM3NeT results, on the other hand, are more representative of the true expected sensitivity of the instrument. Here, the achieved constraints on $f$ mostly depend on the hardness of the fitted input spectra, that is, the predicted $\gamma$-ray flux at the highest energies (> 10 TeV). The large statistical spread in the achieved constraints (indicated by the green bands in Figs. 11, 14, and 15 and by the blue and red bars in Fig. 12) is caused by the small number of neutrinos that is expected to be detectable with KM3NeT, even for the most promising sources.

In the following, we briefly discuss the results obtained for the individual studied sources. A more general assessment of the sensitivity of KM3NeT to the sources studied in this work, except for Westerlund 1, can be found in [22].

**Vela X:** For Vela X we found the strongest discrimination potential for all sources investigated in this study. This is largely due to its full visibility for KM3NeT and its large $\gamma$-ray flux in the 10–100 TeV range. The large $\Delta$TS values obtained with the CTA data sets for the leptonic input model are a result of the strongly curved measured $\gamma$-ray spectrum, which is easier to reproduce with an IC model. While it appears clear that a purely hadronic origin of the $\gamma$-ray emission is already excluded by studies in the $\gamma$-ray domain (e.g. [28,62]), our results indicate that, in the case of a purely leptonic origin, a combined analysis of $\gamma$-ray and neutrino data may help in ruling out a potential small contribution to the $\gamma$-ray flux from hadronic processes. Specifically, the obtained average 68% credible interval constrains the hadronic fraction $f$ to below 15%.

**RX J1713.7–3946:** Compared to Vela X, the spectrum of RX J1713.7–3946 exhibits less curvature, which reflects in the flat $\Delta$TS curves obtained with CTA for this source. As the expected neutrino flux is slightly lower than for
Vela X, the resulting constraints on the hadronic fraction *f* are weaker and the statistical spread is larger. Nevertheless, as the physical processes responsible for the γ-ray emission from RX J1713.7-3946 are not yet totally clear, the addition of neutrino data from KM3NeT may yield important new constraints.

**Westerlund 1:** Compared to the other sources studied in this work, the stellar cluster Westerlund 1 has a relatively hard energy spectrum without a strong cut-off, which leads – in a hadronic scenario – to an expected flux of neutrinos that is detectable with KM3NeT. Therefore, despite Westerlund 1 being a significantly extended source (which implies a larger background contamination), a combined analysis can, on average, distinguish between the leptonic and hadronic models investigated here. However, as the statistical spread is very large, the limits obtained from the different pseudo data sets differ strongly.

**eHWC J1907+063:** In contrast to the other sources, eHWC J1907+063 is modelled using a Gaussian spatial model and has the largest 68% flux containment radius. Because it is furthermore visible for only 54% of the time for KM3NeT, the constraints obtained for this source are the weakest in this study. This is reflected in the overlapping 90% quantile intervals of the fitted hadronic fractions *f* for the leptonic and hadronic input models (see Fig. 12). For the case of a hadronic input model, the CTA data sets yield some separation power as the IC model, being intrinsically curved, cannot reproduce well the measured spectrum which shows only a small curvature. Nevertheless, we conclude that even with 200 h of CTA data and 10 years of KM3NeT data, a differentiation between a leptonic and a hadronic scenario is challenging.

## 4 Conclusion

In this work, we explore the prospects of a combined analysis of γ-ray and neutrino measurements with the upcoming CTA and KM3NeT facilities. To this end, we demonstrate the capability of the open-source γ-ray analysis package GAMMAPY to also analyse data from neutrino telescopes such as KM3NeT, even though these observe the sky in a conceptually different way compared to γ-ray telescope arrays such as CTA.<fnref n="10" /> For the combined analysis, we use publicly available IRFs for CTA and generate custom IRFs for KM3NeT. We fit measured γ-ray spectra of four prototypical Galactic γ-ray sources to obtain both a leptonic model and a hadronic

model. Using either of these models as well as the IRFs as input, we calculate the expected γ-ray and neutrino flux and generate pseudo data sets for both CTA and KM3NeT for all sources. Using a binned likelihood analysis approach, we derive from these pseudo data sets average expected constraints on the contribution of hadronic emission processes to the total emission.

We find that our combined analysis approach enables a consistent modelling of the γ-ray and neutrino observations. Due to simplifying assumptions (e.g. fixed spatial models), the sensitivity to constrain the physical mechanism responsible for the γ-ray emission is often dominated by the KM3NeT data; we note that this will not generally be the case when employing more sophisticated models. For sources that emit a sufficiently large flux of γ rays at high energies (> 10 TeV), our results indicate that a combined analysis of 200 h of CTA and 10 year of KM3NeT data will be able to differentiate between a leptonic and a hadronic emission scenario.

This work constitutes the first demonstration of a general framework for jointly analysing event-level data from γ-ray and neutrino telescopes in a likelihood formalism. Here, we have applied it to simulated data from the CTA and KM3NeT observatories, for a selection of Galactic γ-ray sources, with the aim of deriving expected constraints on the contribution of hadronic emission processes. This approach may be complemented with others, for example the analysis of multi-wavelength data (e.g. in the radio or X-ray domain). On the other hand, the combination of γ-ray and neutrino data can be applied to many more questions. For example, there is increasing evidence that active galactic nuclei – a prominent extragalactic γ-ray source class – are also neutrino sources, see for example [3,65]. In this regard, it is worth noting that the data recorded with both the CTA and KM3NeT observatories will be made publicly available – opening a bright future for multi-messenger analyses.

We release along with this paper a software repository that provides the necessary code, in the form of JUPYTER notebooks, to reproduce the presented results and figures [66]. It can be found at the following URL: https://zenodo.org/record/8298464.

**Acknowledgements** The KM3NeT Collaboration acknowledges the financial support of the funding agencies: Agence Nationale de la Recherche (contract ANR-15-CE31-0020), Centre National de la Recherche Scientifique (CNRS), Commission Européenne (FEDER fund and Marie Curie Program), LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001), Paris Île-de-France Region, France; Shota Rustaveli National Science Foundation of Georgia (SRNSFG, FR-22-13708), Georgia; The General Secretariat of Research and Innovation (GSRI), Greece; Istituto Nazionale di Fisica Nucleare (INFN), Ministero dell'Università e della Ricerca (MIUR), PRIN 2017 program (Grant NAT-NET 2017W4HA7S) Italy; 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); National

<footnote n="10">¹⁰ We note that, since the initiation of this work, several improvements to both the GAMMAPY package and the 'GADF' data format specifications have further eased the analysis of data from wide-field instruments, which besides neutrino telescopes also include γ-ray observatories that employ particle sampling detectors, like the HAWC experiment [63,64].</footnote>
Authority for Scientific Research (ANCS), Romania; Grants PID2021-124591NB-C41, -C42, -C43 funded by MCIN/AEI/ 10.13039/501100011033 and, as appropriate, by "ERDF A way of making Europe", by the "European Union" or by the "European Union NextGenerationEU/PRTR", Programa de Planes Complementarios I+D+I (refs. ASFAE/2022/023, ASFAE/2022/014), Programa Prometeo (PROMETEO/2020/019) and GenT (refs. CIDEGENT/2018/034, /2019/043, /2020/049, /2021/23) of the Generalitat Valenciana, Junta de Andalucía (ref. SOMM17/6104/UGR, P18-FR-5057), EU: MSC program (ref. 101025085), Programa María Zambrano (Spanish Ministry of Universities, funded by the European Union, NextGenerationEU), Spain; The European Union's Horizon 2020 Research and Innovation Programme (ChETEC-INFRA - Project no. 101008324). T. Unbehaun, L. Mohrmann, and S. Funk are members of the CTA Consortium. The paper has been reviewed by the CTA Consortium Speakers and Publications Office. This research has made use of the "Prod 5" CTA instrument response functions provided by the CTA Consortium and Observatory, see https://www.cta-observatory.org/science/cta-performance for more details. This research made use of the ASTROPY (http://www.astropy.org; [67,68]) and MATPLOTLIB (http://www.matplotlib.org; [69]) software packages. The authors gratefully acknowledge the compute resources and support provided by the Erlangen National High Performance Computing Center (NHR@FAU). The publication of this

analysis in the open science regime under FAIR principles [70] is pursued in the EOSC Future (https://eoscfuture.eu) and ESCAPE (https://projectescape.eu) projects that received funding from the European Union's Horizon Europe and 2020 programme under Grant Agreement nos. 101017536 and 824064, respectively.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, 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://creativecommons.org/licenses/by/4.0/. Funded by SCOAP³.

## Appendix A: Additional input model fits

Input model fits for all sources are shown in Fig. 13. The best-fit parameter values are summarised in Table 2. In all fits, we have adopted a distance to the source as listed in Table 1. For the PD model, we have assumed an ambient gas density of 1 cm⁻³. While this density is certainly not a valid assumption for all of the studied sources, this has no impact on the analysis results presented in Sect. 3.2, as the sum of the γ-ray emission predicted by the leptonic and hadronic models is always constrained to match the total observed flux. Different gas densities would thus only imply different normalisations of the primary particle spectra, but have no effect on the derived constraints on the hadronic fraction.
<figure><img src="image_14.png" /><figcaption>**Fig. 13** Summary of the γ-ray input models used for this study. **a** Vela X (already shown in Fig. 6), flux points taken from [28]; **b** RX J1713.7–3946, flux points taken from [31]; **c** Westerlund 1, flux points taken from [35]; **d** eHWC J1907+063, flux points taken from</figcaption></figure>

[39]. For every source, both a PD and an IC model are fitted to the flux points. The goodness-of-fit of the $\chi^2$ fits is indicated in the legend. The prediction of the muon-neutrino flux based on the best-fit PD model is shown as a dotted line

<table><caption>**Table 2** Best-fit parameter values of input models</caption><thead><tr><th>Model</th><th>Parameter</th><th>Unit</th><th>Vela X<sup>a</sup></th><th>RX J1713.7–3946</th><th>Westerlund 1</th><th>eHWC J1907+063</th></tr></thead><tbody><tr><td rowspan="3">PD</td><td>A</td><td>10<sup>35</sup>eV<sup>-1</sup></td><td>(4±45)×10<sup>-5</sup></td><td>1.50 ± 0.07</td><td>6.2 ± 1.5</td><td>6.4 ± 0.5</td></tr><tr><td>Γ</td><td>–</td><td>-4.9 ± 9.9</td><td>2.00 ± 0.08</td><td>1.92 ± 0.35</td><td>2.23 ± 0.09</td></tr><tr><td>E<sub>cut</sub></td><td>TeV</td><td>14 ± 22</td><td>110 ± 30</td><td>890 ± 2340</td><td>540 ± 220</td></tr><tr><td rowspan="3">IC</td><td>A</td><td>10<sup>33</sup>eV<sup>-1</sup></td><td>0.0059 ± 0.0008</td><td>0.283 ± 0.013</td><td>1.17 ± 0.19</td><td>1.20 ± 0.10</td></tr><tr><td>Γ</td><td>–</td><td>0.58 ± 0.40</td><td>2.74 ± 0.12</td><td>2.65 ± 0.49</td><td>3.15 ± 0.12</td></tr><tr><td>E<sub>cut</sub></td><td>TeV</td><td>24 ± 6</td><td>49 ± 12</td><td>960 ± 11200</td><td>380 ± 290</td></tr></tbody></table>

See Eq. 1 for a definition of the parameter values

## Appendix B: Likelihood scan results for all sources

Likelihood scan results for Vela X and RX J1713.7–3946 are shown in Fig. 14, while those for Westerlund 1 and eHWC J1907+063 are shown in Fig. 15.

<footnote><sup>a</sup>The value of Γ for the PD model for Vela X is extreme. Furthermore, the large uncertainties of all parameters of this model indicate a strong correlation between them. These results reflect that the γ-ray emission of Vela X is likely dominantly of leptonic origin, and the strongly curved spectrum can be fitted with the PD model with difficulty only</footnote>
<figure><img src="image_15.png" /></figure>

RX J1713.7-3946

<figure><img src="image_16.png" /><figcaption>**Fig. 14** Profile likelihood scan results for the sources Vela X and RX J1713.7-3946. The plots on the left hand side display the results for a purely leptonic (IC) input model ($f_{\text{in}} = 0$), whereas the plots on</figcaption></figure>

<figure><img src="image_17.png" /></figure>

RX J1713.7-3946

<figure><img src="image_18.png" /></figure>

the right hand side show those for a purely hadronic (PD) input model ($f_{\text{in}} = 1$). For more details, please refer to the caption of Fig. 11
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 15** Same as in Fig. 14, but for the sources Westerlund 1 and eHWC J1907+063</figcaption></figure>

**Appendix C: Example of a PDF for the hadronic fraction $f$**

In Fig. 16, we show an example of the PDF $\Phi(f)$ constructed from the likelihood scan (defined in Eq.) for an individual pseudo experiment in the hadronic scenario of Vela X. The pseudo experiment is hand-picked to show one of the few cases where the upper limit $f_{\text{max}}^{90}$ does not coincide with 1.0. We construct the limits such that all likelihood values included in the shaded interval are larger than those outside. For this reason, the intersections of the limits with the PDF are at the same height $y^{90}$, which results in the smallest 90% credible interval possible.

<figure><img src="image_23.png" /><figcaption>**Fig. 16** PDF for the hadronic fraction $f$ for a single pseudo experiment, constructed from the corresponding likelihood scan. The vertical line at the peak of the PDF marks the best-fit value $\hat{f}$, while the lines at $f_{\text{min}}^{90}$ and $f_{\text{max}}^{90}$ show the lower and upper bound of the 90% credible interval, respectively. The shaded area between the bounds contains 90% of the total area between 0 and 1. The horizontal dashed line at $y^{90}$ shows that the intersections of the PDF with the bounds occur at the same height</figcaption></figure>
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-01010-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.J. 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. C **76**, 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. Tachi-hara, 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. Alves Garre 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. Workman et al., Prog. Theor. Exp. Phys. **2022**, 083C01 (2022). https://doi.org/10.1093/ptep/ptac097
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. Adams et al., 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őcz, 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. Unbehaun¹,*, L. Mohrmann²,*, S. Funk¹

## KM3NeT Collaboration*

S. Aiello³, A. Albert⁴,⁵⁸, S. Alves Garre⁵, Z. Aly⁶, A. Ambrosone⁸,⁷, F. Ameli⁹, M. Andre¹⁰, E. Androutsou¹¹, M. Anghinolfi¹², M. Anguita¹³, L. Aphecetche¹⁴, M. Ardid¹⁵, S. Ardid¹⁵, H. Atmani¹⁶, J. Aublin¹⁷, C. Bagatelas¹¹, L. Bailly-Salins¹⁸, B. Baret¹⁷, Z. Bardačová⁴⁰, S. Basegmez du Pree¹⁹, Y. Becherini¹⁷, M. Bendahman¹⁶,¹⁷, F. Benfenati²⁰,²¹, M. Benhassi⁷,²², D. M. Benoit²³, E. Berbee¹⁹, V. Bertin⁶, S. Biagi²⁴, M. Boettcher²⁵, M. Bou Cabo²⁶, J. Boumaaza¹⁶, M. Bouta²⁷, M. Bouwhuis¹⁹, C. Bozza⁷,²⁸, R. M. Bozza⁷,⁸, H. Brânzaş²⁹, F. Bretaudieu¹⁴, R. Bruijn¹⁹,³⁰, J. Brunner⁶, R. Bruno³, E. Buis¹⁹,³¹, R. Buompane⁷,²², J. Busto⁶, B. Caiffi¹², D. Calvo⁵, S. Campion⁹,³², A. Capone⁹,³², F. Carenini²⁰,²¹, V. Carretero⁵, T. Cartraud¹⁷, P. Castaldi²⁰,³³, V. Cecchini⁵, S. Celli⁹,³², L. Cerisy⁶, M. Chabab³⁴, M. Chadolias¹, A. Chen³⁵, S. Cherubini²⁴,³⁶, T. Chiarusi²⁰, M. Circella³⁷, R. Cocimano²⁴, J. A. B. Coelho¹⁷, A. Coleiro¹⁷, R. Coniglione²⁴, P. Coyle⁶, A. Creusot¹⁷, A. Cruz³⁸, G. Cuttone²⁴, R. Dallier¹⁴, Y. Darras¹, A. De Benedittis⁷, B. De Martino⁶, V. Decoene¹⁴, R. Del Burgo⁷, L. S. Di Mauro²⁴, I. Di Palma⁹,³², A. F. Díaz¹³, D. Diego-Tortosa²⁴, C. Distefano²⁴, A. Domi¹⁹,³⁰, C. Donzaud¹⁷, D. Dornic⁶, M. Dörr³⁹, E. Drakopoulou¹¹, D. Drouhin⁴,⁵⁸, R. Dvornický⁴⁰, T. Eber¹, E. Eckerová⁵⁵, A. Eddymaoui¹⁶, T. van Eeden¹⁹, M. Eff¹, D. van Eijk¹⁹, I. El Bojaddaini²⁷, S. El Hedri¹⁷, A. Enzenhöfer⁶, G. Ferrara²⁴,³⁶, M. D. Filipović⁴¹, F. Filippini²⁰,²¹, L. A. Fusco²⁸, J. Gabriel⁴², T. Gal¹, J. García Méndez¹⁵, A. Garcia Soto⁵, C. Gatius Oliver¹⁹, N. Geißelbrecht¹, H. Ghaddari²⁷, L. Gialanella⁷,²², B. K. Gibson²³, E. Giorgio²⁴, A. Girardi⁹, I. Goos¹⁷, S. R. Gozzini⁵, R. Gracia¹, K. Graf¹, D. Guderian⁵⁹, C. Guidi¹²,⁴³, B. Guillon¹⁸, M. Gutiérrez⁴⁴, H. van Haren⁴⁵, A. Heijboer¹⁹, A. Hekalo³⁹, L. Hennig¹, J. J. Hernández-Rey⁵, F. Huang⁶, W. Idrissi Ibsnalih⁷, G. Illuminati²⁰,²¹, C. W. James³⁸, M. de Jong¹⁹,⁴⁶, P. de Jong¹⁹,³⁰, B. J. Jung¹⁹, P. Kalaczyński⁴⁷, O. Kalekin¹, U. F. Katz¹, N. R. Khan Chowdhury⁵, A. Khatun⁴⁰, G. Kistauri⁴⁸,⁴⁹, F. van der Knaap³¹, A. Kouchner¹⁷,⁵⁰, V. Kulikovskiy¹², R. Kvatadze⁴⁹, M. Labalme¹⁸, R. Lahmann¹, G. Larosa²⁴, C. Lastoria⁶, A. Lazo⁵, S. Le Stum⁶, G. Lehaut¹⁸, E. Leonora³, N. Lessing⁵, G. Levi²⁰,²¹, M. Lindsey Clark¹⁷, F. Longhitano³, J. Majumdar¹⁹, L. Malerba¹², J. Mańczak⁵, A. Manfreda⁷, M. Marconi¹²,⁴³, A. Margiotta²⁰,²¹, A. Marinelli⁷,⁸, C. Markou¹¹, L. Martin¹⁴, F. Marzaioli⁷,²², M. Mastrodocasa⁹,³², S. Mastroianni⁷, S. Miccichè²⁴, G. Miele⁷,⁸, P. Migliozzi⁷, E. Migneco²⁴, P. Mijakowski⁴⁷, M. L. Mitsou⁷, C. M. Mollo⁷, L. Morales-Gallegos⁷, C. Morley-Wong³⁸, A. Mosbrugger¹, A. Moussa²⁷, I. Mozun Mateo⁵¹,⁵², R. Muller¹⁹, M. R. Musone⁷,²², M. Musumeci²⁴, L. Nauta¹⁹, S. Navas⁴⁴, A. Nayerhoda³⁷, C. A. Nicolau⁹, B. Nkosi³⁵, B. Ó Fearraigh¹⁹,³⁰, V. Oliviero⁷,⁸, A. Orlando²⁴, E. Oukacha¹⁷, J. Palacios González⁵, G. Papalashvili⁴⁸, E. J. Pastor Gomez⁵, A. M. Păun²⁹, G. E. Păvălaş²⁹, S. Peña Martínez¹⁷, M. Perrin-Terrin⁶, J. Perronnel¹⁸, V. Pestel⁵², R. Pestes¹⁷, P. Piattelli²⁴, C. Poirè²⁸, V. Popa²⁹, T. Pradier⁴, S. Pulvirenti²⁴, G. Quéméner¹⁸, C. Quiroz¹⁵, U. Rahaman⁵, N. Randazzo³, S. Razzaque⁵³, I. C. Rea⁷, D. Real⁵, S. Reck¹, G. Riccobene²⁴, J. Robinson²⁵, A. Romanov¹²,⁴³, L. Roscilli⁷, A. Saina⁵, F. Salesa Greus⁵, D. F. E. Samtleben¹⁹,⁴⁶, A. Sánchez Losa⁵,³⁷, M. Sanguineti¹²,⁴³, C. Santonastaso⁷,⁵⁴, D. Santonocito²⁴, P. Sapienza²⁴, J. Schnabel¹, M. F. Schneider¹, J. Schumann¹, H. M. Schutte²⁵, J. Seneca¹⁹, N. Sennan²⁷, B. Setter¹, I. Sgura³⁷, R. Shanidze⁴⁸, Y. Shitov⁵⁵, F. Šimkovic⁴⁰, A. Simonelli⁷, A. Sinopoulou³, M. V. Smirnov¹, B. Spisso⁷, M. Spurio²⁰,²¹, D. Stavropoulos¹¹, I. Štekl⁵⁵, M. Taiuti¹²,⁴³, Y. Tayalati¹⁶, H. Tedjditi¹², H. Thiersen²⁵, I. Tosta e Melo³,³⁶, B. Trocme¹⁷, S. Tsagkli¹¹, V. Tsourapis¹¹, E. Tzamariudaki¹¹, A. Vacheret¹⁸, V. Valsecchi²⁴, V. Van Elewyck¹⁷,⁵⁰, G. Vannoye⁶, G. Vasileiadis⁵⁶, F. Vazquez de Sola¹⁹, C. Verilhac¹⁷, A. Veutro⁹,³², S. Viola²⁴, D. Vivolo⁷,²², H. Warnhofer¹, J. Wilms⁵⁷, E. de Wolf¹⁹,³⁰, T. Yousfi²⁷, G. Zarpapis¹¹, S. Zavatarelli¹², A. Zegarelli⁹,³², D. Zito²⁴, J. D. Zornoza⁵, J. Zúñiga⁵, N. Zywucka²⁵
¹ Erlangen Centre for Astroparticle Physics, University of Erlangen-Nuremberg, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany

² Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, Heidelberg 69117, Germany

³ INFN, Sezione di Catania, Via Santa Sofia 64, 95123 Catania, Italy

⁴ CNRS, IPHC UMR 7178, Université de Strasbourg, 67000 Strasbourg, France

⁵ IFIC-Instituto de Física Corpuscular (CSIC-Universitat de València), c/Catedrático José Beltrán 2, 46980 Paterna, Valencia, Spain

⁶ CNRS/IN2P3, CPPM, Aix Marseille Univ, Marseille, France

⁷ INFN, Sezione di Napoli, Complesso Universitario di Monte S. Angelo, Via Cintia ed. G, 80126 Naples, Italy

⁸ Università di Napoli “Federico II”, Dip. Scienze Fisiche “E. Pancini”, Complesso Universitario di Monte S. Angelo, Via Cintia ed. G, 80126 Naples, Italy

⁹ INFN, Sezione di Roma, Piazzale Aldo Moro 2, 00185 Rome, Italy

¹⁰ 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

¹¹ Institute of Nuclear and Particle Physics, NCSR Demokritos, Ag. Paraskevi Attikis, 15310 Athens, Greece

¹² INFN, Sezione di Genova, Via Dodecaneso 33, 16146 Genoa, Italy

¹³ Department of Computer Architecture and Technology/CITIC, University of Granada, 18071 Granada, Spain

¹⁴ Subatech, IMT Atlantique, IN2P3-CNRS, Université de Nantes, 4 rue Alfred Kastler-La Chantrerie, BP 20722, 44307 Nantes, France

¹⁵ 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

¹⁶ Faculty of Sciences, University Mohammed V in Rabat, 4 av. Ibn Battouta, B.P. 1014, R.P. 10000 Rabat, Morocco

¹⁷ CNRS, Astroparticule et Cosmologie, Université Paris Cité, 75013 Paris, France

¹⁸ LPC CAEN, Normandie Univ, ENSICAEN, UNICAEN, CNRS/IN2P3, 6 Boulevard Maréchal Juin, 14050 Caen, France

¹⁹ Nikhef, National Institute for Subatomic Physics, PO Box 41882, 1009 DB Amsterdam, The Netherlands

²⁰ INFN, Sezione di Bologna, v.le C. Berti-Pichat 6/2, 40127 Bologna, Italy

²¹ Dipartimento di Fisica e Astronomia, Università di Bologna, v.le C. Berti-Pichat 6/2, Bologna 40127, Italy

²² Dipartimento di Matematica e Fisica, Università degli Studi della Campania “Luigi Vanvitelli”, Viale Lincoln 5, 81100 Caserta, Italy

²³ E. A. Milne Centre for Astrophysics, University of Hull, Hull HU6 7RX, UK

²⁴ INFN, Laboratori Nazionali del Sud, Via S. Sofia 62, Catania 95123, Italy

²⁵ Centre for Space Research, North-West University, Private Bag X6001, Potchefstroom 2520, South Africa

²⁶ Instituto Español de Oceanografía, Unidad Mixta IEO-UPV, C/ Paranimf 1, 46730 Gandia, Spain

²⁷ Faculty of Sciences, University Mohammed I, BV Mohammed VI, B.P. 717, R.P. 60000 Oujda, Morocco

²⁸ Dipartimento di Fisica, Università di Salerno e INFN Gruppo Collegato di Salerno, Via Giovanni Paolo II 132, 84084 Fisciano, Italy

²⁹ ISS, Atomistilor 409, 077125 Măgurele, Romania

³⁰ Institute of Physics/IHEF, University of Amsterdam, PO Box 94216, 1090 GE Amsterdam, The Netherlands

³¹ TNO, Technical Sciences, PO Box 155, 2600 AD Delft, The Netherlands

³² Dipartimento di Fisica, Università La Sapienza, Piazzale Aldo Moro 2, 00185 Rome, Italy

³³ Dipartimento di Ingegneria dell'Energia Elettrica e dell'Informazione "Guglielmo Marconi", Università di Bologna, Via dell'Università 50, 47521 Cesena, Italy

³⁴ Physics Department, Faculty of Science Semlalia, Cadi Ayyad University, Av. My Abdellah, P.O.B. 2390, 40000 Marrakech, Morocco

³⁵ School of Physics, University of the Witwatersrand, Private Bag 3, 2050 Johannesburg, Wits, South Africa

³⁶ Dipartimento di Fisica e Astronomia "Ettore Majorana", Università di Catania, Via Santa Sofia 64, Catania 95123, Italy

³⁷ INFN, Sezione di Bari, Via Orabona 4, 70125 Bari, Italy

³⁸ International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia

³⁹ University Würzburg, Emil-Fischer-Straße 31, 97074 Würzburg, Germany

⁴⁰ Comenius University in Bratislava, Department of Nuclear Physics and Biophysics, Mlynska dolina F1, Bratislava 842 48, Slovak Republic
⁴¹ School of Computing, Engineering and Mathematics, Western Sydney University, Locked Bag 1797, Penrith, NSW 2751, Australia

⁴² IN2P3, LPC, Campus des Cézeaux 24, Avenue des Landais, BP 80026, 63171 Aubière Cedex, France

⁴³ Università di Genova, Via Dodecaneso 33, 16146 Genoa, Italy

⁴⁴ Dpto. de Física Teórica y del Cosmos & C.A.F.P.E., University of Granada, 18071 Granada, Spain

⁴⁵ NIOZ (Royal Netherlands Institute for Sea Research), PO Box 59, Den Burg, 1790 AB Texel, The Netherlands

⁴⁶ Leiden Institute of Physics, Leiden University, PO Box 9504, 2300 RA Leiden, The Netherlands

⁴⁷ National Centre for Nuclear Research, 02-093, Warsaw, Poland

⁴⁸ Department of Physics, Tbilisi State University, 3, Chavchavadze Ave., 0179 Tbilisi, Georgia

⁴⁹ Institute of Physics, The University of Georgia, Kostava str. 77, 0171 Tbilisi, Georgia

⁵⁰ Institut Universitaire de France, 1 Rue Descartes, 75005 Paris, France

⁵¹ IN2P3, 3, Rue Michel-Ange, 75794 Paris 16, France

⁵² LPC, Campus des Cézeaux 24, Avenue des Landais, BP 80026, 63171 Aubière Cedex, France

⁵³ Department Physics, University of Johannesburg, PO Box 524, 2006 Auckland Park, South Africa

⁵⁴ 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

⁵⁵ Czech Technical University in Prague, Institute of Experimental and Applied Physics, Husova 240/5, Prague 110 00, Czech Republic

⁵⁶ Laboratoire Univers et Particules de Montpellier, Place Eugène Bataillon-CC 72, 34095 Montpellier Cédex 05, France

⁵⁷ Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Remeis Sternwarte, Sternwartstraße 7, 96049 Bamberg, Germany

⁵⁸ Université de Haute Alsace, Rue des Frères Lumière, 68093 Mulhouse Cedex, France

⁵⁹ Institut für Kernphysik, University of Münster, Wilhelm-Klemm-Str. 9, 48149 Münster, Germany