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

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 G AMMAPY , 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 G AMMAPY package to generate pseudo data sets, where we assume 200 hours 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%.

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 hours of CTA observations and 10 years of KM3NeT detector operation.We then apply a threedimensional 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%.

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][7][8][9][10][11][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, 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 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 currentgeneration 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, which targets the detection of high-energy astrophysical neutrinos with TeV-PeV energies.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][18][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 reconstruction of the direction of the incoming neutrino.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 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 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 format1 used for the publicly available CTA IRFs.Then we have employed the GAMMAPY2 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 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].
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], 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 100 TeV 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 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 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 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 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, 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.

Visibility fraction
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.

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

Generation of KM3NeT IRFs
The KM3NeT IRFs are based on extensive simulations of neutrinos and anti-neutrinos 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 multiple 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][52][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 θ 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 (IACTcentred) 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 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.
subsequently averaged, assuming equipartition of the source flux between the two.

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.Fig. 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 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).
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 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.

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 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] 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, where A denotes the amplitude, E 0 the reference energy, Γ the spectral index, E cut the cut-off energy, and β the cutoff 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 χ 2 fit, keeping the reference energy and cut-off strength fixed (at E 0 = 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 scenarios based on γ-ray data alone.We note that the addition of lower-energy radio or X-ray data would further constrain the fit, but regard this as beyond the scope of this work.

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 hours, 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 Poissondistributed random counts are drawn based on those values.
As an example, in Figs.7 and 8, we show projections of one 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 200 hours.The black data points visualise one random Poisson realisation, drawn from the model predictions.
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  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.
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 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.

Likelihood analysis
The analysis is performed using the binned likelihood formalism implemented in the GAMMAPY package 7 .Leptonic and hadronic models are fitted to the generated pseudo data sets8 by minimising the 'Cash statistic' [58] C where denotes the total likelihood to observe the generated data, and is the Poisson probability to measure n i events in pixel i, given a model prediction ν i (ξ ) that depends on the model parameters ξ 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 x, while the other model parameters are optimised in each step.A confidence interval can then be derived from the difference in 'test statistic', where ξ are the parameter values for which L is 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.

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 ξ p describing the proton population and the parameters ξ e for the electron distribution (cf.Eq. 1).Our parameter of interest is then the 'hadronic fraction' where I had and I lep 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 P to the Cash statistic, where and f 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 ξ ) while maintaining a hadronic fraction close to f scan ± ∆ 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. 9 Values of A p = 0.1 and ∆ 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 ξscan the best-fit parameter values for a given f scan , 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 χ 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 ∆ TS by generating and fitting a large number (≫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) Φ( f ) from the likelihood ratio L ( ξscan )/L ( ξ ).By definition (cf.Eq. 5), which we use to derive the PDF as where f ≡ f ( ξscan ) and c is a normalisation constant that can be determined by requiring 1 0 Φ( f ) d f = 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 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 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 Appendix C.

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

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 ξ p and ξ 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 γ-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 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 − p) 2 /∆ p2 , where p and ∆ 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 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.

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, ∆ 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 ∆ 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 individual 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 .

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,  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.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 γ-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 γ-ray flux), whereas the predicted γ-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 γ-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 γ-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 γ-ray flux in the 10-100 TeV range.The large ∆ TS values obtained with the CTA data sets for the leptonic input model are a result of the strongly curved measured γ-ray spectrum, which is easier to reproduce with an IC model.While it appears clear that a purely hadronic origin of the γ-ray emission is already excluded by studies in the γ-ray domain (e.g.[28,62]), our results indicate that, in the case of a purely leptonic origin, a combined analysis of γ-ray and neutrino data may help in ruling out a potential small contribution to the γ-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−3946exhibits less curvature, which reflects in the flat ∆ 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 leadsin 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 yr of KM3NeT data, a differentiation between a leptonic and a hadronic scenario is challenging.

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. 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 yr 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 multiwavelength 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.    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 [39].For every source, both a PD and an IC model are fitted to the flux points.The goodness-of-fit of the χ 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.a 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.Vela X ( f in = 1.0)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 f , while the lines at f 90 min and f 90 max 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.

Fig. 2
Fig.2Background 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.θ 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. 4
Fig.4Comparison 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.

Fig. 5
Fig.5Comparison 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.

Fig. 6
Fig.6Fit 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.

Fig. 7
Fig. 7 Counts map of a pseudo CTA data set based on the PD model for Vela X with 200 hours 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. 8 have been extracted.

Fig. 9
Fig.9Counts 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.10have been extracted (same region as in Fig.7).

3 W e s t e r l u n d 1 fin = 0 fin = 1 Fig. 12
Fig.12Summary 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.

Fig. 13
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[39].For every source, both a PD and an IC model are fitted to the flux points.The goodness-of-fit of the χ 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 1
Galactic γ-ray sources investigated in this work.

Table 2
Best-fit parameter values of input models.