# A hierarchical Bayesian approach to point source analysis in high-energy
neutrino telescopes

Francesca Capel Max Planck Institute for Physics
Boltzmannstraße 8, 85748 Garching  Julian Kuhlmann Technical University of
Munich
James-Franck-Str. 1, 85748 Garching  Max Planck Institute for Physics
Boltzmannstraße 8, 85748 Garching  Christian Haack Erlangen Centre for
Astroparticle Physics
Nikolaus-Fiebiger-Str. 2, 91058 Erlangen  Martin Ha Minh Technical University
of Munich
James-Franck-Str. 1, 85748 Garching  Hans Niederhausen Michigan State
University
288 Farm Lane, 48824 East Lansing MI  Lisa Schumacher Erlangen Centre for
Astroparticle Physics
Nikolaus-Fiebiger-Str. 2, 91058 Erlangen

###### Abstract

We propose a novel approach to the detection of point-like sources of high-
energy neutrinos. Motivated by evidence for emerging sources in existing data,
we focus on the characterisation and interpretation of these sources. The
hierarchical Bayesian model is implemented in the Stan platform, enabling
computation of the posterior distribution with Hamiltonian Monte Carlo. We
simulate a population of weak neutrino sources detected by the IceCube
experiment and use the resulting data set to demonstrate and validate our
framework. We show that even for the challenging case of sources at the
threshold of detection and using limited prior information, it is possible to
correctly infer the source properties. Additionally, we demonstrate how
modelling flexible connections between similar sources can be used to recover
the contribution of sources that would not be detectable individually, going
beyond what is possible with existing stacking methods.

High energy astrophysics (739) — Astronomical methods (1043) — Neutrino
astronomy (1100) — Astrostatistics (1882) — Bayesian statistics (1900) —
Hierarchical models (1925)

††journal: ApJ††facilities: Computations were performed on the HPC system
Raven at the Max Planck Computing and Data Facility.††software: stan and
cmdstanpy (Stan Development Team, 2024), astropy (Astropy Collaboration et
al., 2013, 2018), numpy (Harris et al., 2020), matplotlib (Hunter, 2007),
scipy (Virtanen et al., 2020), h5py (Collette, 2013; Collette et al., 2023),
arviz (Kumar et al., 2019)

(LABEL:)

##  1 Introduction

Neutrino astronomy is in an exciting period, with the discovery of
astrophysical neutrinos confirmed, but the search for their sources still
ongoing (Kurahashi et al., 2022). Recent results from the IceCube
Collaboration present evidence for the association of high-energy neutrinos
with the blazar TXS 0506+056 (Aartsen et al., 2018a, b), the Seyfert galaxy
NGC 1068 (Abbasi et al., 2022a) and the Galactic plane (Abbasi et al., 2023).
Independent analyses making use of public information also claim significant
associations of neutrino events with blazars (Giommi & Padovani, 2021; Buson
et al., 2022, 2023), tidal disruption events (TDEs, Stein et al. 2021; van
Velzen et al. 2021; Reusch et al. 2022), the Seyfert galaxies NGC 4151 and NGC
3079 (Neronov et al., 2024), and the Cygnus region (Neronov et al., 2023).
However, few reports have crossed the $5\sigma$ significance threshold
typically used to define detections and the physical interpretation of these
results remains challenging.

The approach to searching for point sources in neutrino data used by the
IceCube Collaboration makes use of hypothesis testing techniques in a
likelihood-based frequentist framework, as described in Braun et al. (2008,
2010). The reconstructed event directions, energies, and angular uncertainties
are used to distinguish source and background through a likelihood ratio test,
comparing null and signal hypotheses. The significance of a potential source
location is then calculated as a $p$-value by comparing the observed test
statistic value in experimental data to the test statistic distribution
expected under the null model.

We present an alternative approach to point source searches within the
framework of Bayesian hierarchical modelling. The motivation is to make the
most of existing data. Our framework can handle fits of complex models with
large numbers of free parameters. As such, more information from both theory
and experiment can be brought together resulting in more interpretable
statistical analyses. With several large-scale neutrino observatories either
in operation (IceCube: Aartsen et al. 2017, Antares: Ageron et al. 2011,
Baikal-NT: Belolaptikov et al. 1997), in development (KM3Net: Adrián-Martínez
et al. 2016; Margiotta 2022; Baikal-GVD: Avrorin et al. 2021; Dvornicky 2023,
P-ONE: Agostini et al. 2020) or planned in the future (IceCube-Gen2: Aartsen
et al. 2021, TRIDENT: Ye et al. 2023), we also focus on developing methods
that can adapt and scale as we learn more about neutrino sources. As such, we
focus on the characterisation of sources in addition to their discovery.

In this work, we introduce our method in Section 2 and demonstrate its
performance on simulated data in Section 3. We then discuss the performance of
our framework in the context of existing methods in Section 4 before
concluding in Section 5. The code used in this work and relevant examples are
available in the
hierarchical_nu111https://github.com/cescalara/hierarchical_nu (to be made
available upon publication) python package.

##  2 Methods

Our approach is centred on the derivation of a hierarchical or multi-level
likelihood function that captures the key phenomenology of astrophysical
neutrino production and detection. We start with the high-level model
parameters and connect to the observables: the neutrino energies and arrival
directions, along with their respective reconstruction uncertainties. We
describe the key aspects of the physical model and the statistical
implementation here.

###  2.1 Physical model

We consider point-like sources of astrophysical neutrinos with a power-law
spectrum defined between source frame energies $E^{\prime}_{\mathrm{min}}$ and
$E^{\prime}_{\mathrm{max}}$

$\frac{\differential[2]{n{}}}{\differential{E}\differential{t}}\propto
E^{-\gamma_{\mathrm{s}}{}},$ (1)

where $\gamma_{\mathrm{s}}$ is the spectral index. The isotropic source
luminosity with this energy range leads to an energy flux at Earth

$F{}=\frac{L}{4\pi D_{L}{}^{2}(z)}$ (2)

in the redshifted energy range $E_{\text{min, max}}=E^{\prime}{}_{\text{min,
max}}/(1+z)$, where $D_{L}{}(z)$ is the luminosity distance at redshift $z$.
It then follows that the differential flux from a single source at Earth is
given by

$\frac{\differential[3]{n{}}}{\differential{E}\differential{t}\differential{A}}=\frac{Lk_{\gamma}E^{-\gamma_{\mathrm{s}}{}}}{4\pi
D_{L}(z)^{2}}=\phi_{\mathrm{s}}\Bigg{(}\frac{E}{E_{0}}\Bigg{)}^{-\gamma_{\mathrm{s}}{}},$
(3)

where $k_{\gamma}$ is defined such that the spectrum is normalised to $L$ in
the source frame and in the final step we introduce the normalisation energy,
$E_{0}$, and summarise the differential flux at $E_{0}$ as
$\phi_{\mathrm{s}}$. Throughout this work, we assume a flat $\Lambda$CDM
cosmology with $H_{0}=70$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$,
$\Omega_{m}=0.3$ and $\Omega_{\Lambda}=0.7$. We choose a standard power-law
spectral model for point sources here to allow for straightforward comparison
to other methods. However, more complex modelling possibilities are
implemented in our framework and we plan to explore their application in
future work.

In addition to point sources, we also consider two diffuse components as
possible sources of neutrino emission: a diffuse astrophysical background and
an atmospheric background.

Individual sources of interest typically belong to a population of objects
with similar astrophysical properties. In point-source searches, we typically
only hope to resolve some fraction of the total population of sources.
Additionally, individual sources or lists of sources may not be solely
responsible for the observed astrophysical flux. With these factors in mind,
we model a diffuse astrophysical component that accounts for any astrophysical
flux that cannot be associated with resolved point sources. To avoid
assumptions on the cosmological distribution or luminosities of these unknown
sources, we simply model this contribution as an isotropic spectrum over the
whole sky described as

$\frac{\differential[4]{n{}}}{\differential{E}\differential{t}\differential{A}\differential{\omega{}}}=\phi_{\mathrm{d}}{}\Bigg{(}\frac{E}{E_{0}}\Bigg{)}^{-\gamma_{\mathrm{d}}{}},$
(4)

where $\phi_{\mathrm{d}}$ is the differential flux normalisation at $E_{0}$,
and $\gamma_{\mathrm{d}}$ is the spectral index of the bounded power law
spectrum that is defined between $E_{\text{min}}$ and $E_{\text{max}}$. This
spectral model is standard in diffuse astrophysical neutrino analyses, but
could be adapted as new information becomes available (see e.g. Naab et al.
2023).

Another important diffuse background, particularly at energies $E<100$
$\mathrm{TeV}$, is that due to atmospheric neutrinos, produced by the
interactions of cosmic rays in our Earth’s atmosphere. The atmospheric
neutrino flux depends on the zenith angle and the spectrum is not well-
described by a single power law. Therefore, we use MCEq (Fedynitch et al.,
2015) to model the atmospheric arrival direction distribution and spectra. We
use the H4a cosmic ray flux model described in Gaisser (2012), the atmospheric
density profile implemented in NRLMSISE-00 (Picone et al., 2002) and SIBYLL
2.3c to model hadronic interactions (Riehn et al., 2017). Furthermore, we
assume that the normalisation of the atmospheric flux is not exactly known and
parameterise it as $\Phi_{\mathrm{a}}$.

High-energy neutrino telescopes measure the secondary Cherenkov radiation
produced when incoming neutrinos interact in a large instrumented volume of
water or ice. The energies and arrival directions of the primary neutrinos are
not observable, but the reconstructed energies and directions of secondaries
serve as a proxy. The reconstruction of these observables yields an associated
uncertainty. Here, we consider the IceCube neutrino observatory to demonstrate
our approach, but it is straightforward to extend our framework to other
experiments. In particular, we consider the publicly available data of track-
like events from muon neutrino interactions (IceCube Collaboration, 2021) that
is used in Aartsen et al. (2020a), but with limited information on the
provided instrument response functions (IRFs). The IRFs consist of the
effective area, $A_{\mathrm{eff}}$, as function of neutrino energy, $E$, and
declination, $\delta$. Further, the energy resolution and point spread
function are available as a tabulated mapping from $E$ and $\delta$ to
$\hat{E}$ and $\hat{\omega}$. The corresponding public data set lists
reconstructed muon events with energy $\hat{E}$, and direction $\hat{\omega}$.
Further details are discussed in Appendix A.

###  2.2 Statistical formalism

The aim of our statistical framework is to quantify the association of
neutrinos with possible sources, whilst simultaneously inferring the physical
properties of these sources. We expect sources to share similar
characteristics across their classes or populations, making a hierarchical
framework a natural approach.

Our method builds on a hierarchical mixture model formalism for cross-
identification (see e.g. Budavári & Loredo 2015), extended to account for
reconstructed event energies as well as directions, important selection
effects in the observed samples, and relevant source properties (Capel &
Mortlock, 2019). The likelihood has the form of an inhomogeneous Poisson point
process, the rate of which is a mixture model over the different possible
source contributions. This can be written as (Streit, 2010)

$\mathcal{L}(E{},\Theta)=e^{-\bar{n}_{\nu}^{\mathrm{tot}}{}}\prod_{i{}=1}^{n_{\nu}}\sum_{j{}=1}^{N_{\mathrm{c}}}r_{j}{}(\hat{E}{}_{i}{},\hat{\omega}{}_{i}{}|E{}_{i}{},\Theta_{j}),$
(5)

for $n_{\nu}$ detected neutrino events and $N_{\mathrm{c}}$ model components
(where $N_{\mathrm{c}}$ = $N_{\mathrm{s}}$ \+ 2 diffuse components). The
expected number of detected neutrino events, $\bar{n}_{\nu}^{\mathrm{tot}}$,
is given by the sum over all model components such that
${\bar{n}_{\nu}^{\mathrm{tot}}{}=\sum_{j{}=0}^{N_{\mathrm{c}}{}}\bar{n}_{j}}$.
Similarly, the expected number of events from point sources is given by
${\bar{n}_{\nu}{}=\sum_{j{=0}}^{N_{\mathrm{s}}{}}\bar{n}_{j}}$. We summarise
the vector of all high-level parameters as $\Theta$, and the subset of those
relevant for individual model components as $\Theta_{j}{}$. All energies are
summarised as $E{}$. For $N_{\mathrm{s}}$ point sources,
${\Theta_{j}{}=\\{L_{j}{},\gamma_{j}{},\omega{}_{j}{},D_{j}{}\\}}$, whereas
for the diffuse astrophysical and atmospheric components
${\Theta_{j}{}=\\{\Phi_{\mathrm{d}}{},\gamma_{\mathrm{d}}{}\\}}$ and
${\Theta_{j}{}=\\{\Phi_{\mathrm{a}}{}\\}}$, respectively. $r_{j}{}$ is the
$j$-th component’s rate parameter, relating to the number of expected events,
$\bar{n}_{\nu}^{\mathrm{tot}}$, through integration as shown in Eq.
eq:app_nex. The relationship between parameters in Eq. eq:likelihood is
summarised graphically in Fig. 1 for the case of point source components.

$\omega{}_{j}{}$$D_{j}{}$$\gamma_{j}{}$$L_{j}{}$$E{}_{i}{}$$\lambda_{i}{}$$\Phi_{\mathrm{s}}{}$$\hat{\omega}{}_{i}{}$$\hat{E}{}_{i}{}$$\bar{n}_{\nu}{}$

$N_{\mathrm{s}}{}$

Figure 1: Summary of the likelihood for point source components. Open and
shaded circles are used to show model parameters and fixed observables
respectively, with the arrows showing the connections between them. The boxes
are used to separate the $N_{\mathrm{s}}$ source compenents and $n_{\nu}$
events. The parameters are explained in the text. The levels from upper to
lower show the source properties, latent true parameters and quantities that
depend on the IRF. The latent source labels, $\lambda_{i}$, are discrete
parameters marking the mixture model component that each event belongs to.

All model parameters are left free, but the highest-level parameters also have
associated prior distributions. As a default, we make use of weakly
informative priors for these parameters, as shown in Table 1. Such priors
allow us to include our knowlegde on reasonable physical values, while not
driving the results of the eventual inference. In Section 3.3, we also
demonstrate the application of more informative priors, as could be motivated
by, e.g., multi-messenger information or theoretical predictions. Where
relevant, we verify that the results are robust to reasonable variations in
the choice of priors.

| Distribution | Units
---|---|---
$L/L_{j}{}$ | $\mathrm{Lognorm}(8.0\times 10^{43},4.0)$ | $\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$
$\gamma_{\mathrm{s}}{}/\gamma_{j}{}$ | $\mathrm{Norm}(2.0,0.25)$ | –
$\Phi_{\mathrm{d}}{}$ | $\mathrm{Lognorm}(5.4\times 10^{-8},0.3)$ | ${\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$
$\gamma_{\mathrm{d}}{}$ | $\mathrm{Norm}(2.5,0.04)$ | –
$\Phi_{\mathrm{a}}{}$ | $\mathrm{Lognorm}(3.0\times 10^{-5},0.08)$ | ${\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$
Table 1: Default prior assumptions for model hyperparameters, given as the
distributions used and the corresponding $\mu$ and $\sigma$ values.
$L/L_{j}{}$ is bounded between $[0,10^{52}]$
$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$, $\Phi_{\mathrm{d}}$ and
$\Phi_{\mathrm{a}}$ are also bounded between $[10^{-10},10^{-7}]$ and
$[0,3\times 10^{-4}]$ ${\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$. Both
$\gamma_{\mathrm{s}}{}/\gamma_{j}{}$ and $\gamma_{\mathrm{d}}$ are bounded
between 1 and 4.

For a given data set of $n_{\nu}$ neutrino events, we perform inference on the
model parameters in a Bayesian framework. The posterior distribution is
proportional to the likelihood given in Eq. eq:likelihood multiplied by the
joint prior summarised in Table 1. We generate samples from this posterior
using Stan (Stan Development Team, 2024), which implements an efficient
variant of a Markov chain Monte Carlo algorithm called Hamiltonian Monte
Carlo. This class of algorithm guarantees convergence to the target
distribution in the limit of infinite samples. Here, we use a set of
diagnostics to assess convergence in the case of finite samples. In
particular, we run 4 separate chains of 2000 samples each (1000 warm-up
samples to tune the algorithm and find the target distribution, 1000 actual
samples of the target distribution) and require an effective sample size,
$n_{\mathrm{eff}}>1000$ and a Gelman-Rubin statistic, $\hat{R}<1.1$ with no
divergent transitions for all model parameters in our analyses (Gelman et al.,
2013).

###  2.3 Interpretation of results

The results of the fits using Stan take the form of samples representing the
joint posterior distribution over all model parameters, also including latent
parameters. Useful summaries of the fit parameters can be derived from these
samples, such as the most probable value and highest posterior density
interval (HDI), representing the “best-fit” value and its associated
uncertainty. We evaluate the “goodness-of-fit” using posterior predictive
checks (PPCs), which involve generating simulated data under the assumptions
of the fitted model and comparing to the observed data to check the inferences
are reasonable (Gelman et al., 2013, Chapter 6.3).

As shown in Fig. 1, each neutrino event has a latent true energy parameter,
$E{}_{i}{}$. The marginal posterior of this parameter provides additional
information useful to the interpretation of the origin of this event,
considering all possible model components and uncertainties, which is non-
trivial to infer from the reconstructed event energies, $\hat{E}{}_{i}{}$,
alone.

To directly address the possible association of neutrino events with source
components of the model, we introduce additional derived parameters that
capture the relevant information. Each neutrino event has a discrete label
parameter, $\lambda_{i}{}$, which identifies its source component. While the
$\lambda_{i}{}$ are not explicity sampled during inference, we can compute the
marginal posterior distribution for $\lambda_{i}{}$,
$\Pr(\lambda_{i}{}=j{}|\hat{E}{},\hat{\omega}{})$, as shown in Appendix B.
This marginal posterior gives the association probability for each
event–source component pair given the available data, allowing for a more
direct and insightful intepretation of the results.

We can also quantify discovery and sensitivity within this Bayesian framework.
For individual event–source associations, we can confidently claim an
association if

$\Pr(\lambda_{i}{}=j{}|\hat{E}{},\hat{\omega}{})>\alpha,$ (6)

where $\alpha$ is a threshold probability to be defined. However, the
information in individual events is naturally limited, and therefore if using
small event sample sizes it is important to study the prior dependence of the
association probability (see Capel et al. 2023 for an example). A useful
parameter that summarises our expectation over the whole sample of possible
event–source associations is the expected number of detected neutrinos from
the $j$-th point source, $\bar{n}_{j}$. We define the detection of a source as

$\Pr(\bar{n}_{j}\geq 1|\hat{E}{},\hat{\omega}{})>\alpha.$ (7)

In the case of a non-detection or when quantifying the sensitivity, we can
find the HDI of $\Phi_{j}$ or $L_{j}$ for a given probability level, $\alpha$,
and report the corresponding upper limit.

The probability thresholds, $\alpha$, defined above are somewhat arbitrary and
are not directly comparable to a $p$-value from a frequentist analysis. It is
possible to calibrate $\alpha$ via repeated simulations to give some desired
coverage, true detection rate or false detection rate (Betancourt, 2018).
However, within a Bayesian framework, the resulting posterior probability for
event–source associations or expected number of events is itself the main
result. For the purpose of this work we choose $\alpha=0.95$ (similar
definitions can be found in, e.g., Aggarwal et al. 2021; Sottosanti et al.
2021). This is not equivalent to a $p$-value of 0.05, and is importantly a
statement about the alternate hypothesis (point-source model components)
rather than the null hypothesis (diffuse background model components).

##  3 Application

In this work, we focus the application of our method to a population of weak
sources that are below the detection threshold of existing methods.

###  3.1 Simulated data set

We simulate a population of neutrino sources using popsynth (Burgess & Capel,
2021), characterised by the physical properties introduced in Section 2.1: a
redshift, $z$, a luminosity, $L$, a spectral index, $\gamma$, and a direction,
$\omega$. We choose a source redshift evolution based on the shape of the star
formation rate given in Madau & Dickinson (2014) and a luminosity function
that is described by a broken power law. Source spectral indices follow a
Gaussian distribution and sources are disributed isotropically on the sky. The
choices regarding the source population properties are detailed in Table 2.

| Distribution
---|---
$z$ |  $\frac{\differential{N}}{\differential{V}}=100\frac{1+4.8z}{1+(z/2.7)^{3.9}}$ ${\mathrm{Gpc}}^{-3}\text{\,}\mathrm{sr}$
L | $\frac{\differential{N}}{\differential{L}}=\begin{cases}CL^{-2}&\mbox{if }L\leq L_{\mathrm{br}}\\\ CL^{-3}L_{\mathrm{br}}&\mbox{if }L>L_{\mathrm{br}}\end{cases}$ ${\mathrm{erg}}^{-1}\text{\,}\mathrm{s}$
$\gamma$ | $\frac{\differential{N}}{\differential{\gamma}}=\mathrm{Norm}(2.0,0.25)$
$\omega$ |  $\frac{\differential{N}}{\differential{\omega}}=\frac{1}{4\pi}$ ${\mathrm{sr}}^{-1}$
Table 2: Source population properties and their distributions. We use
$L_{\mathrm{br}}=5\times 10^{44}$ $\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$ and
$C$ is defined such that the distribution is normalised to 1. The the
following ranges are also used: $z\in[0,6]$, $L\in[5\times 10^{43},5\times
10^{48}]$ $\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$ and $\gamma\in[1,4]$.

Our source population is set up to provide a generic but relevant test
scenario that is not connected to a particular class of astrophysical sources.
We ensure that our simulated population is consistent with general population
constraints (Murase & Waxman, 2016; Capel et al., 2020) in that it does not
exceed the measured diffuse astrophysical neutrino flux and sources are
neither too rare nor too bright. Beyond $z\sim 0.5$, the neutrino signal from
the population is effectively diffuse ($\bar{n}_{\nu}\lesssim 1$ for
individual sources).

To generate our simulated neutrino data set, we select sources in the Northern
hemisphere ($\delta>-5^{\circ}$) with $F\geq 5\times 10^{12}$
$\mathrm{erg}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$, giving
$N_{\mathrm{s}}=6$. The properties of these selected sources are given in
Table 3.

| $(RA,\delta)_{j}$ | $z$ | $L_{j}$ | $\gamma_{\mathrm{s}}{j}$ | $n_{j{}}$
---|---|---|---|---|---
| [$\mathrm{\SIUnitSymbolDegree}$] |  | [$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$] |  |  $\geq 300$ $\mathrm{GeV}$ |  $\geq 1$ $\mathrm{TeV}$
1 | (183.6, 5.5) | 0.04 | $5.8\times 10^{43}$ | 2.05 | 14 | 11
2 | (56.0, 29.3) | 0.43 | $5.8\times 10^{45}$ | 2.30 | 7 | 6
3 | (273.2, 28.5) | 0.13 | $5.3\times 10^{44}$ | 2.52 | 5 | 5
4 | (23.8, 15.4) | 0.06 | $6.6\times 10^{43}$ | 1.99 | 4 | 4
5 | (179.6, 20.8) | 0.17 | $1.6\times 10^{45}$ | 1.62 | 2 | 2
6 | (325.0, 62.3) | 0.12 | $1.81\times 10^{44}$ | 1.81 | 1 | 1
Table 3: Properties of the simulated neutrino point sources. Sources are
listed in order of decreasing number of simulated events, $n_{j{}}$, that are
included in the event selection used here.

Fig. 2 shows the selected sources in comparison with the sensitivity and
discovery potential of skyllh, which implements a frequentist approach that is
typically used in neutrino point source searches (Wolf, 2019; Bellenghi et
al., 2023a). The sensitivity flux is calculated as the median 90% CL upper
limit on the flux normalisation over repeated simulations of only background.
The discovery potential is defined as the flux that leads to a $5\sigma$
deviation from the background expectation in 50% of repeated simulations. The
same curves from Aartsen et al. (2020b) are also shown, but as these have been
calculated using more detailed detector Monte Carlo that is not publicly
available, they do not represent a fair comparison in this work.

All sources are below the estimated discovery potential, but it is important
to note that the sensitivity and discovery potential fluxes are calculated for
a fixed power-law spectrum with index $\gamma_{\mathrm{s}}{}=2$, while the
point sources simulated here have $\gamma_{j}$ in the range $[1.62,2.52]$. In
Fig. 2, we also show the corresponding number of events in our simulation.
This visualiation better accounts for the different spectral index
assumptions, but the number of events required for a detection in a given
simulation will depend on their energies. In this work, we only investigate
one possible realisation of our simulation and Fig. 2 is only intended to
provide a rough comparison with existing work.

Figure 2: Upper panel: The differentital flux normalisation at 1
$\mathrm{TeV}$ for simulated point sources. The SkyLLH sensitivity assuming an
$E^{-2}$ power-law spectrum as reported in Bellenghi et al. (2023a) is shown
for comparison, as this calculation makes use of the same publicly available
IRF that is used here. We also plot an estimation of the discovery potential
by rescaling the sensitivity as a function of declination. This rescaling is
calculated as the ratio of the discovery potential to the sensitivity,
according to the results presented in Aartsen et al. (2020b), which are also
shown for comparison. Lower panel: The total expected number of neutrino
events for simulated point sources is shown in comparison to that required by
the sensitivity and discovery potential curves introduced above.

All other sources from the population are treated as a sub-dominant
contribution to the diffuse flux and are not simulated explicitly. We include
a diffuse astrophysical flux with $\phi_{\mathrm{d}}=1.8\times 10^{-18}$
${\mathrm{GeV}}^{-1}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{sr}}^{-1}$
and $\gamma_{\mathrm{d}}=2.5$ based on recently reported observations (Abbasi
et al., 2022b; Naab et al., 2023; Abbasi et al., 2024). The atmospheric flux
is also included using MCEq as described in Section 2.1. All source components
are simulated over a wide energy range from $E_{\text{min}}=100$
$\mathrm{GeV}$ to $E_{\text{max}}=100$ $\mathrm{PeV}$. We simulate the
detection of neutrinos from these source components using the provided IRF and
place a cut on the reconstructed muon energy of $\hat{E}_{\mathrm{min}}=300$
$\mathrm{GeV}$. This choice is to include thresholding effects due to events
with true energies $E{}<300$ $\mathrm{G}\mathrm{e}\mathrm{V}$ being mis-
reconstructed above $\hat{E}_{\mathrm{min}}$ and vice versa. For this
simulated data set, we simulate $n_{j{}}$ as the most likley integer value
based on $\bar{n}_{j{}}$, for a clearer interpretation and discussion of the
different source cases in the next sections. A summary of the resulting data
set is shown in Fig. 3.

Figure 3: The simulated data set. The left panel shows a skymap of the event
directions with the selected source positions overlaid. The right panel shows
the distribution of events in $\hat{E}$ and $\sin(\delta)$, cf. Fig. 1 in
IceCube Collaboration (2021).

###  3.2 Single source analysis

To demonstrate the analysis of individual sources in our framework we focus on
the point source that has the most detected events in our simulated data set,
source #1 from Table 3. We make a selection on our data set, defining a
circular region of interest (ROI) on the sky centred on the source location
with a radius of $5^{\circ}$222We verified that for the source modelling
assumptions used here and the angular resolution provided by the IRF, this
radius is large enough to avoid impacting the results..

We use our method to fit the data in this ROI, using the priors detailed in
Table 1. We show the results for the marginal posterior of the source
parameters in Fig. 4, demonstrating the correct reconstruction of the known
true parameter values.

Figure 4: The joint posterior density of $L$ and $\gamma_{\mathrm{s}}$ with
contours showing the 30%, 60% and 90% HDIs. The true parameter values are
indicated by the horizontal and vertical solid lines. The upper and right
panels show the marginal distributions in comparison with the prior.

The marginal posterior of the expected number of neutrino events from source
#1, $\bar{n}_{1}$, is given in Fig. 5. To validate the correct reconstruction
of $\bar{n}_{1}$, we also fit data sets where we add in source events one by
one. Fig. 5 also shows the resulting $\bar{n}_{1}$ posterior for
$n_{\nu}{}=[0,5,10,14]$. We can see that in the case of $n_{\nu}{}=0$, the
$\bar{n}_{1}$ posterior is strongly peaked at zero and gradually moves away in
a consistent manner as more point source events are added to the simulation.
Using the definition of source detection given in Eq. eq:discovery_nex with
$\alpha=0.95$, we can identify this source above the detection threshold in
the simulated dataset that includes all 14 source events. We also performed
fits to 100 simulations of background-only events, verifying that none of
these result in false point-source detections according to our detection
criterion.

Our framework further enables us to give posterior association probabilities
of single events to each source component, as detailed in Section 2.3. We
summarise this information in Fig. 6 in terms of both the event positions and
energies. We see that the higher energy events that are closer to the source
have a larger association probability, as expected. Point source events that
have a low association probability are mostly mis-reconstructed as belonging
to the atmospheric background component, which dominates at lower energies.
Using the event-based definition of association given in Eq.
eq:discovery_event, only the event with the second highest $\hat{E}$ value
passes the threshold with an association proability of 0.99. The event with
the highest $\hat{E}$ value is close behind, with an association probability
of 0.95.

In Fig. 6, we see the complementary information that is available in terms of
the individual event–source association probabilities. For example, Fig. 5
shows that the overall $\bar{n}_{1}$ posterior is consistent with 14 events
for source #1, but Fig. 6 shows that we have probable individual associations
for only about half of these events, typically at higher energies. The
remaining contribution to $\bar{n}_{1}$ comes from many events with smaller
association probabilities. We can also see that for events that have a $\sim
50\%$ association probability, the posterior for $E$ can have two peaks,
corresponding to the possible latent $E$ that would best fit the source or
background components.

As means of verifying goodness-of-fit we perform PPCs. Histograms of $\hat{E}$
and $\hat{\omega}$ for 100 generated data sets from the posterior predictive
distribution are shown in Fig. 7, indicating no obvious signs of mismodelling
from visual inspection. We also quantified the discrepancy in each bin by
calculating a $p$-value, and found only relatively large $p$-values that
followed a uniform distribution between 0 and 1.

Figure 5: Marginal posterior density of the expected number of point source
events for source #1, $\bar{n}_{1}$. The number of selected source events in
each fit is colour coded, with $n_{1}$ = [0, 5, 10, 14] shown. Figure 6:
Colour-coded posterior association probabilities of events to source #1. Left
panel: Spatial distribution of events. The marker size does not represent the
angular resolution and the point source events are highlighted in red borders.
Right panel: Marginal posteriors for the true neutrino energies of events
included in the ROI. Markers on the top row indicate $\hat{E}$ for all events
included in the analysis, with dashed connecting lines to the posteriors for
the events from the tested point source. Figure 7: Distribution of observed
and posterior predictive data. Shaded bands show 30%, 60%, and 90% quantiles
from darkest to lightest.

###  3.3 Joint analysis of multiple sources

We now perform joint fits of all six selected point sources. Motivated by
computational considerations, we make a further cut of $\hat{E}{}\geq 1$
$\mathrm{TeV}$ on our simulated dataset, as only 4 of the total 33 point
source events are lost in this case (see Table 3). Additionally, based on the
results in Section 3.2, we do not expect the lower energy events to contribute
to the detectability of these sources. Our framework is designed to allow for
flexible definitions of possible connections between individual sources, and
we explore three cases here as examples: 1) All sources share the same
luminosity, $L$, and spectral index, $\gamma_{\mathrm{s}}$; 2) All sources
have individual and unknown $L_{j}{}$ and $\gamma_{j}{}$; and 3) Informative
priors on the individual $L_{j}$ and $\gamma_{j}{}$. We note that our
framework can also handle the case between 1) and 2), including some balance
between global and individual source parameters, and this is discussed further
in Section 4.

For the joint fits, all information shown in Section 3.2 is available and can
be analysed independently for all sources. The difference between fitting all
sources separately is that the joint fit allows the results of one source to
influence the others, according to the details of the model for their
population. For brevity, we focus here on few key results that highlight the
implications of the three example cases introduced above. Fig. 8 shows the
marginal posterior distributions for the point source luminosity, $L$ and the
expected number of events from each source, $\bar{n}_{\nu}$, for each of the
three studied cases.

Figure 8: The marginal posteriors for $L$, $\gamma_{\mathrm{s}}{}$ and
$\bar{n}_{j}$ for all 6 point sources, shown by the different colours. The
upper, middle and lower rows show the results for cases 1), 2) and 3),
respectively. The true values are indicated as solid vertical lines and the
priors are shown in grey, where relevant.

For the first case of shared $L$ and $\gamma_{\mathrm{s}}$, we see that the
source with most events, source #1, drives the values of these parameters. As
this source has a relatively low $L$333Despite having a low $L$, as source #1
is the closest source, it still has the highest number of detected events in
the simulated data set., $\bar{n}_{j}$ is reduced for all sources other than
source #4, which has a similar $L$. For this case, $\gamma_{\mathrm{s}}$ is
relatively unconstrained due to the low event numbers, so there is only a
slight information gain of the posterior relative to the prior. The overall
impact is the underestimation of $\bar{n}_{j}$ for harder sources, such as
sources #3, #5 and #6.

For the case of individual free source parameters, $L_{j}$ and $\gamma_{j}$,
the variety of the sources is better captured. While the $L_{j}{}$ and $n_{j}$
posteriors are all consistent with the corresponding true values, the
distributions remain relatively unconstrained for sources #2, #4 and #6. The
$\gamma_{j}$ posteriors show some shift away from the prior for sources #5 and
#1, but are otherwise rather similar to the posterior from the shared
parameter case above.

Finally, we investigate using informative priors for the individual $L_{j}{}$
and $\gamma_{j}$. Cases 1) and 2) both use shared priors for the source
parameters, as described in Table 1. We now suppose that further information
is available on the $L_{j}$ and $\gamma_{j}$ for the individual sources. This
could be in the form of, e.g., neutrino observations of other sources, multi-
messenger information, or theoretical predictions. We use log-normal priors
for $L_{j}$, centred on the true values and with a width of 2.0. For the
$\gamma_{j}$, we use normal priors centred on the true values with
$\sigma_{\gamma}=0.1$. While informative, these priors still allow for the
possibility that $\bar{n}_{j{}}=0$ for all sources, and the allowed parameter
ranges leave room for the data to drive the posterior results. We see that in
this case, we have better reconstruction of both $L_{j}$ and $\bar{n}_{j{}}$
for all sources, but the most impact is seen for sources #2, #4, and #6, which
were the least constrained in case 2). The $\gamma_{j}$ posteriors mostly
follow the priors, with no strong additional constraints from the data.

We further summarise the above results for the three example cases by
comparing the results for the total number of expected point source events
across all sources, $\bar{n}_{\nu}$, as shown in Fig. 9. We see that assuming
the $L$ and $\gamma_{\mathrm{s}}$ are shared by all sources as in case 1), we
underestimate the total contribution of these point sources as $\sim 1/2$ of
the true value. For case 2) we improve the result, but still underestimate the
contribution as only $\sim 3/4$ of the truth and the extra information in case
3) is necessary to recover the contribution from all point sources. As
$\Pr(\bar{n}_{\nu}{}\geq 1|\hat{E}{},\hat{\omega}{})=1.0$ for all cases, it
makes more sense to state the $\bar{n}_{\nu}$ threshold for which this
expression is $>0.95$. For cases 1), 2) and 3) this is 7, 9 and 14,
respectively. We can interpret this is a confident detection of this many
events from all point sources included in the fit.

Figure 9: The marginal posterior distributions for $\bar{n}_{\nu}$ from all 6
point sources for the three example cases discussed in Section 3.3.

We give the $\Pr(\bar{n}_{j{}}\geq 1|\hat{E}{},\hat{\omega}{})$ values for
each source considered in Table 4. Interestingly, for case 1), both sources #1
and #4 are detectable according to the definition in Eq. eq:discovery_nex with
$\alpha=0.95$, but only source #1 remains detectable for cases 2) and 3). This
can be understood as sources #1 and #4 have very similar $L_{j}$ and
$\gamma_{j}$, so it is beneficial to assume shared parameters in this case,
especially given the strength of source #1. However, the overall detectability
of all sources is better served by considering their individual properties,
and the value of $\Pr(\bar{n}_{j}\geq 1|\hat{E}{},\hat{\omega}{})$ increases
significantly for all sources other than #4.

Source | $\Pr(\bar{n}_{j{}}\geq 1|\hat{E}{},\hat{\omega}{})$
---|---
Case 1 | Case 2 | Case 3
#1 | 0.998 | 0.998 | 1.000
#2 | 0.000 | 0.098 | 0.566
#3 | 0.158 | 0.777 | 0.859
#4 | 0.989 | 0.305 | 0.568
#5 | 0.007 | 0.661 | 0.769
#6 | 0.110 | 0.316 | 0.436
Table 4: Impact of the three different modelling assumptions described in the
text on the detectability of all sources considered.

The event–source association probabilities generally reflect the results
discussed above. It is worth mentioning that source #5 produces an event with
$\hat{E}{}=1.3$ $\mathrm{PeV}$ in our simulation, which has an association
probability of 0.92, 0.95 and 0.99 for cases 1), 2) and 3), respectively. As
this is a faint source with only 2 simulated events and more extreme
properties, it makes sense that the association probability increases as we
allow the sources to have independent parameters and include more prior
information.

##  4 Discussion

The results shown in Section 3.2 and 3.3 demonstrate the validity of our
method and introduce the possible analyses that can be performed within this
framework. In addition to parameter estimation, we address how probabilities
of interest can be quantified, such as the event–source association
probability and the probability that a source or population contributes at
least $\bar{n}_{\nu}$ events to the data. The most relevant output will
typically depend on the details of the application to different physical
scenarios. For example, the event–source association probabilities are likely
most useful for studying high-energy neutrino alert events or analysing the
impact of different source spectral models on possible associations. On the
other hand, the overall expected event contribution is a direct way to
investigate the detectability of sources and their populations that naturally
includes the relevant uncertainties. Example applications could include
characterising a population of hard-spectrum sources that TXS 0506+056-like
sources to (Buson et al., 2022, 2023; Bellenghi et al., 2023b), or
constraining the contribution from a population of NGC 1068-like sources
(Glauch et al., 2023; Saurenhaus & Capel, 2023). With our framework, we aim to
provide a consistent setting for application to these different cases.

In Section 3.3, we also illustrate the power that including more information
into the analysis can have when trying to detect a number of weak sources in
the data. Even for events with a relatively poor energy resolution, we see in
Fig. 9 that including prior information on the spectral index makes a
significant difference when estimating the contribution of these point
sources. For a non-detection, this equates to stronger constraints on the
proposed model, allowing us to make the most of the available data in either
case.

Several recent works in the field of multi-messenger astrophysics explore
Bayesian approaches for individual source–event associations (see e.g. Ashton
et al. 2018, Bartos et al. 2019 and Veske et al. 2021, Kowalski 2021). These
methods still frame the problem as a hypothesis test, preferring Bayes factors
and odds ratios to $p$-values, or using the Bayes factor as a test statistic
in order to draw conclusions. Here, we focus more on addressing the questions
of interest via parameter estimation rather than model comparison. The
foreseen workflow is to be able to develop and refine models that are
consistent with the data, while exploring their implications. In this
workflow, model rejection is a subset of the possible outcomes.

Frequentist hypothesis testing methods are typically used in searches for
neutrino point sources. In particular, likelihood ratio methods as introduced
in Braun et al. (2008) and implemented in SkyLLH (Wolf, 2019; Bellenghi et
al., 2023a) have been used in IceCube analyses to find evidence for TXS
0506+056, NGC 1068 and the Galactic plane as neutrino sources (Aartsen et al.,
2018a, b; Abbasi et al., 2022a, 2023). We see our Bayesian approach as
complemetary to the standard methods in that the focus is on the evaluation
and charaterisation of source models rather than the rejection of the
background-only hypothesis. Due to the different definition of probability in
frequentist and Bayesian statistics and the different goals of these two
methods, it is non-trivial to directly compare their performance. As such, we
only make qualitative comparisons in this work and highlight the complementary
features below.

With our definition of detection and the realisation of the simulated dataset
studied, we can detect source #1 with only weakly informative priors, as well
as source #4 when assuming shared source parameters as in case 1) of Section
3.3. Source #5 has one very high energy event and is therefore detectable for
case 3). Sources #2, #3 and #6 cannot be detected independently as they are
softer (#2 and #3) or happen to only produce lower energy events in this
simulated data set (#6), but their contribution to the population can be
recovered when including additional information, as shown in Fig. 9.

We recognise that our joint source fits are conceptually similar to the
“stacking” technique that has been used in previous point source searches to
increase sensitivity to particular source modelling assumptions (e.g. Glauch
et al. 2023). Case 1) is most similar to a distance-weighted stacking analysis
and case 3) is similar to a flux-weighted stacking analysis. However, an
important difference is that the priors and modelling that we use are set up
in a way such that relevant uncertainties are included and the data can
overrule the prior in more informative cases. This modelling, together with
the interpretation offered by the Bayesian approach results in a more flexible
analysis.

Cases 1) and 2) explore the possibility that all sources are the same
(“complete pooling”) and that all sources and independent (“no pooling”),
respectively. Neither of these assumptions are completely realistic, and in
practice we expect some balance between global and individual source
properties (“partial pooling”), as modelled in our simulation. We see for case
2) that the neutrino data alone does not contain enough information to
significantly constrain the population hyperparameters (the shape of the
luminosity function and spectral index distribution). However, the future
increase in performance expected from planned experiments mentioned in Section
1 will provide the data sets necessary for the “partial pooling” case to be
leveraged, and we plan to explore the impact of these new possibilities in
upcoming work.

Thanks to the implementation of our statistical model with Stan’s Hamiltonian
Monte Carlo algorithm, we are able to perform fits with over $\sim 7000$ free
parameters in Section 3.3, where the majority of the free parameters are the
latent $E$ of the events. In principle, these latent $E$ parameters could be
marginalised over in the likelihood to speed up the fits, but as discussed in
Section 2.3 and demonstrated in Section 3.2, these parameters add to the
interpretability of the results. By including these free parameters here, we
also demonstrate that it is relatively straightforward to add further model
complexity in terms of parameters for both the source modelling and detector
modelling. Large numbers of free parameters can be challenging for optimisers
currently implemented in SkyLLH, as it was not designed to fit the latent $E$
parameters or more complex source models. Markov chain Monte Carlo methods,
such as those used here, could also be implemented in a frequentist setting to
address these challenges.

Another way in which we include more complexity in this work is that we model
the atmospheric and diffuse background components, while keeping the total
number of expected events conserved. In this way, if events are fit to point
sources, the diffuse astrophysical and/or atmospheric components are reduced
to compensate. This approach is interesting when considering the bigger
picture of possible neutrino sources, with competing source populations
expected (e.g., Bartos et al. 2021), in addition to contributions from the
Galactic plane.

A challenge of the frequentist hypothesis testing approach arises when testing
multiple hypotheses, necessitating a trial correction factor. Consequently,
source lists for studies have to be limited or the discovery threshold raised,
leading to less use of the data and decreased sensitivity. Furthermore, it is
non-trivial to keep track of trial factors across independent analyses of the
same data. While the Bayesian methods used here do not guarantee a certain
false positive rate, the structure included through the priors and hierarchy
of our model naturally tends to mitigate the effects of outliers (Gelman et
al., 2009). Should a certain false positive rate be desired, it is possible to
calibrate the probability thresholds introduced in Section 2.3 through
repeated simulation and fits, as detailed in Betancourt (2018).

Along with guaranteed coverage and false positive rates, one complementary
aspect of the standard approach is that by focusing on rejecting background,
unexpected signals can be identified even if the alternate hypothesis is not
well-specified. In our approach, mismodelling can be also be studied with PPCs
and used to improve models for a better match with the data. Additionally,
computational challenges mean that our approach is currently better suited to
building and testing specific point source models rather than performing an
uninformed scan across the whole sky. Our implementation in Stan includes
within-chain parallelisation of the likelihood evaluation that can be scaled
according to the available computational resources. As a benchmark, for point
sources near the Equator where event rates are the largest, fits using 10
years of data selected in an ROI of $None$ radius may take up to 2 hours on
144 threads.

##  5 Conclusions

We present a hierarchical Bayesian approach to searching for point sources of
astrophysical neutrinos. Our method is an alternative to existing frequentist
approaches, with a focus on the characterisation of sources and
interpretability of results.

We demonstrate our approach through application to a simulated data set
containing 6 weak point sources hidden in typical backgrounds. Even for low
event numbers of $n_{j{}}\leq 14$, we are able to recover the constribution of
the strongest source and provide constraints on its luminosity and spectral
index. The contribution of the remaining sources can also be inferred by
leveraging similar source properties and more informative priors for their
spectral shape. These results show the potential gain of more complex and
flexible modelling when studying weak sources that is relevant to the expected
signals from current neutrino source candidates, such as blazars or Seyfert
galaxies.

We plan to apply our framework to the existing public data sets, investigate
the impact of more complex source modelling, and study the information gain
from expected future data sets. Our framework is written in a modular way,
such that it can also be extended to other event types and detectors,
providing a useful open-source tool to the community for the evaluation of
different astrophysical models and detector configurations.

## acknowledgments

The authors express their thanks to Chiara Bellenghi and Martin Wolf for very
helpful discussions about the used data set, Lena Saurenhaus for contributing
to the production of Fig. 3. The authors are also grateful to Allen Caldwell
and Daniel Mortlock for feedback on an early version of this manuscript.

Julian Kuhlmann acknowledges support from the DFG through the
Sonderforschungsbereich SFB 1258 “Neutrinos and Dark Matter in Astro- and
Particle Physics” (NDM). Francesca Capel acknowledges the financial support
from the Excellence Cluster ORIGINS, which is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094-390783311 at the early stage of this work.

##  Appendix A Detector model

Using the information provided in the publicly available IRFs (IceCube
Collaboration, 2021), we can calculate the number of expected events in a
sample as

$\bar{n}_{\nu}^{\mathrm{tot}}{}=T\int_{E_{\text{min}}{}}^{E_{\text{max}}{}}\differential{E{}}\int_{\hat{E}_{\mathrm{min}}{}}^{\hat{E}_{\mathrm{max}}{}}\differential{\hat{E}{}}\int_{S_{2}}\differential{\omega{}}\int_{\mathrm{ROI}}\differential{\hat{\omega}{}}\frac{\differential[4]{n{}}}{\differential{E{}}\differential{t}\differential{A}\differential{\omega{}}}A_{\mathrm{eff}}{}(E{},\omega{})\,\mathrm{IRF}(\hat{E}{}|E{},\omega{})\,\mathrm{IRF}(\hat{\omega}{}|\omega{},E{},\hat{E}{}),$
(A1)

where $T$ is the total observation time, ROI is the region of interest and
$S_{2}$ the surface of a sphere. The integrals over solid angle reflect the
modelling of the sources and the cuts on reconstructed direction,
respectively. In fits, the angular resolution is modelled using a Rayleigh
distribution embedded on a sphere, for a derivation see Glauch (2021). In
simulations, $\sigma_{\hat{\omega}}{}$ is sampled from the provided IRFs.

Due to the required differentiability of the model likelihood with respect to
the neutrino energy, we have to interpolate the provided energy resolution,
which is provided as histograms in reconstructed energies, covering half a
decade of neutrino energy, split over three declination bands. For each
histogram $\Pr(\hat{E}{}|E{})$ we fill zero-entries in between non-zero
entries by interpolation and extend the empty flanks by a steep power-law with
index $\pm 15$. After renormalising to unity, we stack the histograms’
logarithm along a new axis of logarithmic neutrino energy and create a two-
dimensional spline representation. Evaluations of this spline at each
$\hat{E}$ over a dense grid of $\log_{10}(E{})$ is handed over to Stan for
interpolation. Similarly, we interpolate the logarithm of effective area over
logarithmic energies.

As the provided instrument response is only valid for events originating as
neutrinos we are restricted to analysis in the Northern hemisphere, as in the
Southern hemisphere the event rate is dominated by atmospheric muons.

##  Appendix B Event–source association

We can construct an association probability of events to source components by
ratios of rate parameters. For every event, we define a categorical variable
$\lambda_{i}{}$, taking on possible values $1,\dots,N_{\mathrm{c}}{}$:

$\Pr(\lambda_{i}{}=j{}|\hat{E}{},\hat{\omega}{},\Theta)=\frac{r_{j}{}(\hat{E}{},\hat{\omega}{}|\Theta_{j}{})}{\sum_{l=1}^{N_{\mathrm{c}}{}}r_{l}(\hat{E}{},\hat{\omega}{}|\Theta_{l})},$
(B1)

This expression evaluated for the point source and marginalised over the joint
posterior is used as colour scale in Fig. 6. We refer for more details to
Streit (2010, Chapter 3.2) and the Stan documentation444https://mc-
stan.org/docs/stan-users-guide/finite-mixtures.html#recovering-posterior-
mixture-proportions.

## References

  * Aartsen et al. (2018a) Aartsen, The IceCube Collaboration, Fermi-LAT, et al. 2018a, Science, 361, eaat1378, doi: 10.1126/science.aat1378
  * Aartsen et al. (2018b) Aartsen, M., Ackermann, M., Adams, J., et al. 2018b, Science, 825, eaat2890 , doi: 10.1126/science.aat2890
  * Aartsen et al. (2020a) —. 2020a, Physical Review Letters, 124, doi: 10.1103/physrevlett.124.051103
  * Aartsen et al. (2017) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017, Journal of Instrumentation, 12, P03012, doi: 10.1088/1748-0221/12/03/p03012
  * Aartsen et al. (2020b) —. 2020b, Physical Review Letters, 124, 051103, doi: 10.1103/physrevlett.124.051103
  * Aartsen et al. (2021) Aartsen, M. G., Abbasi, R., Ackermann, M., et al. 2021, Journal of Physics G: Nuclear and Particle Physics, 48, 060501, doi: 10.1088/1361-6471/abbd48
  * Abbasi et al. (2022a) Abbasi, R., Ackermann, M., Adams, J., et al. 2022a, Science, 378, 538, doi: 10.1126/science.abg3395
  * Abbasi et al. (2022b) —. 2022b, The Astrophysical Journal, 928, 50, doi: 10.3847/1538-4357/ac4d29
  * Abbasi et al. (2023) —. 2023, Science, 380, 1338, doi: 10.1126/science.adc9818
  * Abbasi et al. (2024) —. 2024, Characterization of the Astrophysical Diffuse Neutrino Flux using Starting Track Events in IceCube.  https://arxiv.org/abs/2402.18026
  * Adrián-Martínez et al. (2016) Adrián-Martínez, S., Ageron, M., Aharonian, F., et al. 2016, Journal of Physics G: Nuclear and Particle Physics, 43, 084001, doi: 10.1088/0954-3899/43/8/084001
  * Ageron et al. (2011) Ageron, M., Aguilar, J. A., Samarai, I. A., et al. 2011, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 656, 11 , doi: 10.1016/j.nima.2011.06.103
  * Aggarwal et al. (2021) Aggarwal, K., Budavári, T., Deller, A. T., et al. 2021, The Astrophysical Journal, 911, 95, doi: 10.3847/1538-4357/abe8d2
  * Agostini et al. (2020) Agostini, M., Böhmer, M., Bosma, J., et al. 2020, Nature Astronomy, 4, 913, doi: 10.1038/s41550-020-1182-4
  * Ashton et al. (2018) Ashton, G., Burns, E., Canton, T. D., et al. 2018, The Astrophysical Journal, 860, 6, doi: 10.3847/1538-4357/aabfd2
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Avrorin et al. (2021) Avrorin, A. D., Avrorin, A. V., Aynutdinov, V. M., et al. 2021, Physics of Atomic Nuclei, 84, 513, doi: 10.1134/s1063778821040062
  * Bartos et al. (2019) Bartos, I., Veske, D., Keivani, A., et al. 2019, Physical Review D, 100, 083017, doi: 10.1103/physrevd.100.083017
  * Bartos et al. (2021) Bartos, I., Veske, D., Kowalski, M., Márka, Z., & Márka, S. 2021, The Astrophysical Journal, 921, 45, doi: 10.3847/1538-4357/ac1c7b
  * Bellenghi et al. (2023a) Bellenghi, C., Karl, M., Wolf, M., et al. 2023a, in Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), Vol. 444, 1061, doi: 10.22323/1.444.1061
  * Bellenghi et al. (2023b) Bellenghi, C., Padovani, P., Resconi, E., & Giommi, P. 2023b, The Astrophysical Journal Letters, 955, L32, doi: 10.3847/2041-8213/acf711
  * Belolaptikov et al. (1997) Belolaptikov, I., Bezrukov, L., Borisovets, B., et al. 1997, Astroparticle Physics, 7, 263, doi: 10.1016/s0927-6505(97)00022-4
  * Betancourt (2018) Betancourt, M. 2018, Calibrating Model-Based Inferences and Decisions.  https://arxiv.org/abs/1803.08393
  * Braun et al. (2010) Braun, J., Baker, M., Dumm, J., et al. 2010, Astroparticle Physics, 33, 175, doi: 10.1016/j.astropartphys.2010.01.005
  * Braun et al. (2008) Braun, J., Dumm, J., Palma, F. D., et al. 2008, Astroparticle Physics, 299 , doi: 10.1016/j.astropartphys.2008.02.007
  * Budavári & Loredo (2015) Budavári, T., & Loredo, T. J. 2015, Annual Review of Statistics and Its Application, 2, 113 , doi: 10.1146/annurev-statistics-010814-020231
  * Burgess & Capel (2021) Burgess, J., & Capel, F. 2021, Journal of Open Source Software, 6, 3257, doi: 10.21105/joss.03257
  * Buson et al. (2022) Buson, S., Tramacere, A., Pfeiffer, L., et al. 2022, The Astrophysical Journal Letters, 933, L43, doi: 10.3847/2041-8213/ac7d5b
  * Buson et al. (2023) Buson, S., Tramacere, A., Oswald, L., et al. 2023, Extragalactic neutrino factories.  https://arxiv.org/abs/2305.11263
  * Capel et al. (2023) Capel, F., Kuhlmann, J., Haack, C., et al. 2023, in ICRC 2023, Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), 1576, doi: 10.22323/1.444.1576
  * Capel & Mortlock (2019) Capel, F., & Mortlock, D. J. 2019, Monthly Notices of the Royal Astronomical Society, 484, 2324–2340, doi: 10.1093/mnras/stz081
  * Capel et al. (2020) Capel, F., Mortlock, D. J., & Finley, C. 2020, Physical Review D, 101, 123017, doi: 10.1103/physrevd.101.123017
  * Collette (2013) Collette, A. 2013, Python and HDF5 (O’Reilly)
  * Collette et al. (2023) Collette, A., Kluyver, T., Caswell, T. A., et al. 2023, h5py/h5py: 3.8.0, 3.8.0, Zenodo, doi: 10.5281/zenodo.7560547
  * Dvornicky (2023) Dvornicky, R. 2023, in Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), Vol. 444, 976, doi: 10.22323/1.444.0976
  * Fedynitch et al. (2015) Fedynitch, A., Engel, R., Gaisser, T. K., Riehn, F., & Stanev, T. 2015, Calculation of conventional and prompt lepton fluxes at very high energy.  https://arxiv.org/abs/1503.00544
  * Gaisser (2012) Gaisser, T. K. 2012, Astroparticle Physics, 35, 801, doi: 10.1016/j.astropartphys.2012.02.010
  * Gelman et al. (2013) Gelman, A., Carlin, J. B., Stern, H. S., et al. 2013, Bayesian Data Analysis, 3rd edn. (Chapman and Hall/CRC), doi: 10.1201/b16018
  * Gelman et al. (2009) Gelman, A., Hill, J., & Yajima, M. 2009, Why we (usually) don’t have to worry about multiple comparisons.  https://arxiv.org/abs/0907.2478
  * Giommi & Padovani (2021) Giommi, P., & Padovani, P. 2021, Universe, 7, 492, doi: 10.3390/universe7120492
  * Glauch (2021) Glauch, T. 2021, PhD thesis, Technische Universität München
  * Glauch et al. (2023) Glauch, T., Kheirandish, A., Kontrimas, T., et al. 2023, in Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), Vol. 444, 1052, doi: 10.22323/1.444.1052
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * IceCube Collaboration (2021) IceCube Collaboration. 2021, IceCube Data for Neutrino Point-Source Searches Years 2008-2018, IceCube Neutrino Observatory, doi: 10.21234/CPKQ-K003
  * Kowalski (2021) Kowalski, M. 2021, Nature Astronomy, 5, 732, doi: 10.1038/s41550-021-01431-y
  * Kumar et al. (2019) Kumar, R., Carroll, C., Hartikainen, A., & Osvaldo, M. 2019, Journal of Open Source Software, doi: 10.21105/joss.01143
  * Kurahashi et al. (2022) Kurahashi, N., Murase, K., & Santander, M. 2022, Annual Review of Nuclear and Particle Science, 72, 365
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, Annu. Rev. Astron. Astrophys., 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Margiotta (2022) Margiotta, A. 2022, The KM3NeT infrastructure: status and first results.  https://arxiv.org/abs/2208.07370
  * Murase & Waxman (2016) Murase, K., & Waxman, E. 2016, Physical Review D, D94, 103006, doi: 10.1103/physrevd.94.103006
  * Naab et al. (2023) Naab, R., et al. 2023, in Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), Vol. 444, 1064, doi: 10.22323/1.444.1064
  * Neronov et al. (2024) Neronov, A., Savchenko, D., & Semikoz, D. V. 2024, Physical Review Letters, 132, 101002, doi: 10.1103/physrevlett.132.101002
  * Neronov et al. (2023) Neronov, A., Semikoz, D., & Savchenko, D. 2023, Neutrino signal from Cygnus region of the Milky Way.  https://arxiv.org/abs/2311.13711
  * Picone et al. (2002) Picone, J. M., Hedin, A. E., Drob, D. P., & Aikin, A. C. 2002, Journal of Geophysical Research: Space Physics, 107, SIA 15, doi: 10.1029/2002ja009430
  * Reusch et al. (2022) Reusch, S., Stein, R., Kowalski, M., et al. 2022, Physical Review Letters, 128, 221101, doi: 10.1103/physrevlett.128.221101
  * Riehn et al. (2017) Riehn, F., Dembinski, H., Engel, R., et al. 2017, in ICRC 2017, Proceedings of 35th International Cosmic Ray Conference — PoS(ICRC2017), 301, doi: 10.22323/1.301.0301
  * Saurenhaus & Capel (2023) Saurenhaus, L., & Capel, F. 2023, in Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), Vol. 444, 1509, doi: 10.22323/1.444.1509
  * Sottosanti et al. (2021) Sottosanti, A., Bernardi, M., Brazzale, A. R., et al. 2021, Identification of high-energy astrophysical point sources via hierarchical Bayesian nonparametric clustering.  https://arxiv.org/abs/2104.11492
  * Stan Development Team (2024) Stan Development Team. 2024, Stan Modeling Language Users Guide and Reference Manual.  http://mc-stan.org/
  * Stein et al. (2021) Stein, R., Velzen, S. v., Kowalski, M., et al. 2021, Nature Astronomy, 1, doi: 10.1038/s41550-020-01295-8
  * Streit (2010) Streit, R. L. 2010, Poisson Point Processes: Imaging, Tracking, and Sensing, 1st edn. (Springer New York, NY), doi: 10.1007/978-1-4419-6923-1
  * van Velzen et al. (2021) van Velzen, S., Stein, R., Gilfanov, M., et al. 2021, Establishing accretion flares from massive black holes as a major source of high-energy neutrinos.  https://arxiv.org/abs/2111.09391
  * Veske et al. (2021) Veske, D., Márka, Z., Bartos, I., & Márka, S. 2021, The Astrophysical Journal, 908, 216, doi: 10.3847/1538-4357/abd542
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wolf (2019) Wolf, M. 2019, in Proceedings of 36th International Cosmic Ray Conference — PoS(ICRC2019), Vol. 358, 1035, doi: 10.22323/1.358.1035
  * Ye et al. (2023) Ye, Z. P., Hu, F., Tian, W., et al. 2023, Nature Astronomy, 7, 1497, doi: 10.1038/s41550-023-02087-6