# Nonparametric signal separation in very-high-energy gamma ray observations
with probabilistic neural networks

M. Ullmo  E. Moulin

###### Abstract

An intriguing challenge in observational astronomy is the separation signals
in areas where multiple signals intersect. A typical instance of this in very-
high-energy (VHE, E$\gtrsim$100 GeV) gamma-ray astronomy is the issue of
residual background in observations. This background arises when cosmic-ray
protons are mistakenly identified as gamma-rays from sources of interest,
thereby blending with signals from astrophysical sources of interest. We
introduce a deep ensemble approach to determine a non-parametric estimation of
source and background signals in VHE gamma observations, as well as a
likelihood-derived epistemic uncertainty on these estimations. We rely on
minimal assumptions, exploiting the separability of space and energy
components in the signals, and defining a small region in coordinate space
where the source signal is assumed to be negligible compared to background
signal. The model is applied both on mock observations, including a simple toy
case and a realistic simulation of dark matter annihilation in the Galactic
center, as well as true observations from the public H.E.S.S. data release,
specifically datasets of the Crab nebula and the pulsar wind nebula MSH 15-52.
Our method performs well in mock cases, where the ground truth is known, and
compares favorably against conventional physical analysis approaches when
applied to true observations. In the case of the mock dark matter signal in
the Galactic center, our work opens new avenues for component separation in
this complex region of the VHE sky.

##  1 Introduction

Very-high-energy (VHE, E$\gtrsim$100 GeV) gamma rays are a versatile tool for
probing the most violent phenomena in the Universe. They can reveal
information about the regions they traverse and the mechanisms that produce
them, such as the acceleration of relativistic particles through plasma shocks
and magnetic fields, or processes like dark matter annihilation. Identifying
and characterizing such mechanisms requires a careful analysis of the gamma
ray sources, notably through the examination of their spatial, spectral, and
temporal features, prompting the need for precise analytical tools to
accurately capture them. However, at these very high energies the Earth’s
atmosphere is fully opaque. Thus ground-based telescopes must observe VHE
gamma rays indirectly, one event at a time, by detecting the Cherenkov light
emitted by the extensive air shower produced when a ray interacts with the
atmosphere. These Imaging Atmospheric Cherenkov Telescopes (IACTs) include
instruments like the High Energy Stereoscopic System (H.E.S.S. [1]), Major
Atmospheric Gamma Imaging Cherenkov (MAGIC [2]) telescopes, the Very Energetic
Radiation Imaging Telescope Array System (VERITAS [3]), and the upcoming
Cherenkov Telescope Array (CTA [4]), which promises to revolutionize the field
with unprecedented sensitivity and resolution.

Figure 1: The training setup. The networks output an estimation of
$S_{s},E_{s},S_{b}$ and $E_{b}$ over the region of interest in the coordinate
space. These are combined to make $F_{tot}$, which is then compared to the
observation through the loss given in Eq. (2.4). Optionally, constraints are
imposed on $S_{s}$ and $S_{b}$.

A critical challenge in ground-based gamma-ray astronomy is the background
noise from high-energy cosmic rays (CR) coming from all directions of the sky.
These CR create showers of secondary particles in the Earth’s atmosphere that
produce Cherenkov radiation similar to that of gamma-ray induced showers.
Preventing the telescopes from interpreting these CR-induced events as genuine
gamma-ray events requires advanced data analysis methods to discriminate
between the two when analyzing camera images of the showers. Several
techniques have been developed, from standard analytical methods like the
Hillas parameters method [5, 6] to machine-learning techniques like random
forest algorithms [7] or graph networks [8]. While these methods prove highly
effective and considerably reduce the amount of CR misidentified as gamma
rays, none are as of yet $100\%$ effective. Consequently, processed
observations still contain a residual background that is spatially isotropic
and often stronger than the sources of interest.

To address the problem of residual background, and depending on whether the
position the source is a priori known or not, different statistical methods
have been developed to estimate and remove background signal in the detected
event list. One common technique is the Ring Background method [9], where a
ring-shaped region around the on-source area (region of the observation
containing the source signal) is used to estimate the background. This
approach assumes that the background is isotropic and can be accurately
represented by the measurements in the surrounding region provided all VHE
emissions there are properly excluded for the background estimate within the
on region. Adaptative Ring background methods have been developed when
exclusion regions cover a significant part the field of view [10]. Another
method is the multiple OFF technique [9], which uses multiple OFF-source
regions to obtain a more reliable background estimate by averaging out spatial
fluctuations. Advanced statistical methods, such as the template-background
method [11], involve creating detailed models of the background based on the
distribution of cosmic-ray events across the field of view. Despite these
sophisticated techniques, limitations persist. The assumption of isotropy in
the ring background method may not always hold true, particularly in regions
with complex background structures. The multiple OFF technique, while reducing
statistical fluctuations, requires careful selection of OFF-source regions to
avoid biases. The template background method demands extensive computational
resources and accurate modeling, which can be challenging given the
variability of cosmic-ray events. Additionally, all these methods rely on the
ability to precisely reconstruct the direction and energy of incoming
particles, which is inherently uncertain due to the nature of the Cherenkov
emission process and atmospheric effects. Consequently, while these methods
significantly mitigate background contamination, they do not completely
eliminate it, prompting the need to seek alternative solutions.

Neural networks (NNs [12]) have demonstrated their effectiveness in such
signal separation tasks across diverse fields, including acoustics [13],
seismic monitoring [14], medical imaging [15], and even gamma-ray astronomy
[16].

In this article we introduce a proof of concept for a NN-based framework to
address the issue of CR background contamination, and as a first step for the
more general problem of multiple component separation in VHE astrophysics. Our
approach consists in jointly estimating the mean fields of the signals
corresponding to CR background and source of interest to maximize the
likelihood of a given observation, while providing a statistically-meaningful
uncertainty to account for model limitations as well as statistical noise in
the observations. Our method offers several significant advantages. It avoids
making any prior parametric assumptions about component shapes, making it
ideal for sources with complex shapes, and does not require a training
database, relying solely on the observation of interest for signal and
background estimation; conversely, this means that training is required for
each individual case, contrary to common machine learning approaches.
Additionally, it is flexible and can easily incorporate constraints based on
prior knowledge to improve the estimation process. The paper is organized as
follows: Section 2 addresses the problem setup and the description of our
neural network model. Section 3 presents the data on which the model is
applied. The results are presented in section 4. In section 5 we further
interpret results, address limitations and discuss outlooks of the present
work.

##  2 Method and Neural Network framework

###  2.1 Problem setup and formalism

We consider a given observation, a 3D histogram of events $(e,\vec{x})$ binned
in energy and space, as a Poisson realization of an underlying mean field
$F_{tot}$. We assume $F_{tot}$ is made up of two components, one for source
emission of interest ($F_{s}$) and the other due to the background noise
($F_{b}$). We additionally assume that both $F_{s}$ and $F_{b}$ can be further
decomposed into space and energy components $(S_{s},E_{s})$ and
$(S_{b},E_{b})$, respectively, such that:

$\begin{cases}F_{s}(\vec{x},e)=S_{s}(\vec{x})\times E_{s}(e),\\\
F_{b}(\vec{x},e)=S_{b}(\vec{x})\times E_{b}(e)\\\ \end{cases}$ (2.1)

This gives the signal mixture expressed as :

$F_{tot}(\vec{x},e)=S_{s}(\vec{x})\times E_{s}(e)+S_{b}(\vec{x})\times
E_{b}(e)\,.$ (2.2)

As a convention, we impose that $\sum_{obs}E_{s,b}=1$, such that $E_{s,b}$ are
dimensionless and physical quantities are expressed within $S_{s,b}$.

Networks | $S_{s},S_{b}$ | $E_{s},E_{b}$
---|---|---
Input dimension and coordinate | 2, $\vec{x}$ | 1, $e$
Number of layers | $4$ | $4$
Neurons per layer | $10$ | $5$
Layer Activation function | Leaky ReLU | Leaky ReLU
Final Activation function | ELU | ELU
Output dimension | 1 | 1
Figure 2: Hyperparameter settings in the network architecture used in this
work. All Neural Networks are trained using the Adam optimizer [17] with a
learning rate $lr$=10-3.

We aim to estimate the spatial and spectral shapes $S_{s},E_{s},S_{b}$ and
$E_{b}$ from the 2-component mixture in our observed field of view. To do so,
we build four distinct NNs to emulate the four components. The networks are
taken as simple dense nets (see architectures in table 2) and are trained in
the following manner (see Fig.1): they are provided coordinates as input
($\vec{x}$ for $S_{s}$ and $S_{b}$, $e$ for $E_{s}$ and $E_{b}$) and their
outputs (estimates of $S_{s},E_{s},S_{b}$ and $E_{b}$ over the observed
region) are combined as in Eq. (2.2) to form $F_{tot}^{est}$, which is then
compared to the observation.

Given that we aim to reconstruct $F_{tot}$ such that our observation is
maximally probable with respect to $F_{tot}$, our loss, which we aim to
minimize, is thus chosen to be:

$L=-\langle\log{\mathcal{P}(Obs|F_{tot})}\rangle_{\rm bins}$ (2.3)

Considering probabilities over individual bins of coordinates $(\vec{x},e)$
this becomes:

$L=-\frac{1}{n_{\rm
bins}}\sum_{\vec{x},e}\log{\mathcal{P}(Obs(\vec{x},e)|F_{tot}(\vec{x},e))}\,,$
(2.4)

where $\mathcal{P}(n|m)=\frac{e^{-m}m^{n}}{n!}$ is the Poisson probability of
observing $n$ counts given an average occurrence $m$, consistently with the
expected statistics of an observation.

The network architectures are deliberately designed to be relatively shallow
and low in parameters. This low parametrization imposes a certain smoothness
on the output functions and helps prevent overfitting, which would occur if
the network were complex enough to fit the Poisson noise in the observations.
We anticipate the background to be smoother in space compared to the source,
justifying the use of a lower-parameter architecture. Future research should
explore different architectures, weight initialization methods, learning rates
and gradient descent algorithms, as the current choices have been made
empirically and may benefit from optimization.

###  2.2 Network Training

####  2.2.1 Deep ensemble approach

Our approach essentially equates to a nonparametric fitting of source and
background, with a solution that maximizes likelihood. However, given the high
dimension of our problem we can expect a few complications. First, whether
parameter space has a single global maximum; given the redundancy in the NNs
we can already expect a manifold of parameter configurations that output the
same functions/predictions, however depending on several factors (noise in the
observation, overlapping functions, NN constraints, etc), we can additionally
expect a manifold of equally-likely solutions that correspond to the global
maxima. Additionally, depending on initial conditions we can expect networks
to find and get stuck into local maxima during training.

To account for local maxima cases and to better capture the manifold of
equiprobable maxima, we opt for a deep ensemble approach; instead of training
a single model (made up of four NNs), we train a set of models and observe
their combined results. In our case, we train 10 separate models. For each we
only change initial weights, which are randomized from a Gaussian
distribution. A more comprehensive approach would be to vary NN architectures
for each separate model as well, and experiment with varying set size; it is
let for further study.

To get a final estimate for all functions, we compute a weighted average of
all estimations. We additionally obtain a model uncertainty for all functions
by computing a weighted standard deviation. The weights for each estimation
$w_{i}$ are defined relative to estimation likelihood in the following way:

$w_{i}=\exp{\Bigg{(}\frac{1}{c_{\rm
tot}}\sum_{b}c_{b}\log{\mathcal{P}\Big{(}obs(b)|F_{tot}^{i}(b)\Big{)}}\Bigg{)}}\,,$
(2.5)

where $c_{\rm tot}$ is the total particle count, $c_{b}$ is the particle count
in bin $b$, and $F_{\rm tot}^{i}(b)$ is the total estimated field at bin $b$
for estimation $i$. In other words, every weight $w_{i}$ is itself a weighted
log-average probability per bin of the observation with respect to $F_{\rm
tot}$, where every bin is weighted by its particle count, such that regions of
the coordinate space with higher statistics are conferred more weight.

####  2.2.2 Adding constraints

Additionally we need to impose specific constraints on the model to
differentiate source and background signal (i.e., to break the symmetry in the
model as made explicit by Eq. (2.2)) and to leverage all available knowledge
for more precise estimations. We opt in this work for the following two
constraints.

##### Hat constraint

A first simple way to do this is to define an OFF region, where source signal
is expected to be negligible relative to background. Having defined this
region, we can impose our source signal to be zero within it simply by
multiplying its space component with a "hat" function $H$ such that
$H(\vec{x})=0$ in OFF region, $H(\vec{x})=1$ elsewhere.

For all cases we will be using this hat constraint, and will define the OFF
region as all bins beyond distance $d$ from the center of the observation. For
simplicity’s sake we set $d$ as two thirds of an observation’s field-of-view
radius, but further experiments should attempt varying $d$ or changing OFF
region shape altogether to see how results are affected. Furthermore, to
accommodate cases of extended sources where signal is not expected to
completely disappear within the field of view, we can instead opt for a soft
hat that imposes an upper limit (typically a fixed fraction of the background)
for the source in the OFF region, rather than making it zero.

##### Background threshold

In most realistic cases, source and background are difficult to separate,
typically because of overlapping energy spectra, leading to source signal
contaminating background estimation, and vice-versa. In extreme cases this can
lead to a one-component estimation where source or background is estimated to
be zero. Since using the hat constraint prevents background from being zero, a
complementary constraint is to impose an upper limit on the background’s
spatial estimation so that brighter points are assumed to be either due to
source signal or Poisson noise; thus limiting the chances of a null source
estimation. Since the instrumental background is expected to be spatially
isotropic, this is a hypothesis worth using if we know the observation to be
relatively spatially homogeneous. Supposing homogeneous background, we set the
mean count per spatial bin as threshold:

$T=\big{\langle}\sum_{\rm ebins}obs\big{\rangle}_{\rm sbins}$ (2.6)

Figure 3: Mock observation for the simple toy model case (left panels) and a
simulated observation of dark matter annihilation signals in the Galactic
center (right panels). Observations over spatial and energy bins are plotted
in the top an bottom panels, respectively. Energy and position are in
arbitrary units for the simple mock case.

###  2.3 Rejection Criterion

As certain cases are bound to be too complex for our model, leading to
erroneous estimations, we can perform a simple sanity check to ensure that the
final estimation is plausible. Indeed, one need only compare the likelihood of
the final total estimated field $F_{\rm tot}^{\rm est}$ to that of the one-
component solution $F_{\rm tot}^{\rm 1C}$, which corresponds simply to:

$\displaystyle\begin{cases}F_{tot}^{1C}=S^{1C}\times E^{1C},\\\
S^{1C}=\sum_{\rm ebins}obs,\\\ E^{1C}=\frac{1}{\sum obs}\sum_{\rm
sbins}obs,\\\ \end{cases}$ (2.7)

where the space and energy components simply correspond to the observation
summed over energy and space, respectively (normalized for energy). If
$\mathcal{P}(obs|F_{\rm tot}^{\rm 1C})\geq\mathcal{P}(obs|F_{\rm tot}^{\rm
est})$, the estimation should be rejected as erroneous.

This being said, in cases with low statistics, the one-component solution as
defined above can turn out to have a relatively high likelihood, possibly
higher than the ground truth solution, as it overfits the Poisson noise. In
such cases comparison against a one-component solution with similar
constraints on the smoothness of functions as our two-component models should
be preferred.

A more general approach hypothesizing more than one source should similarly
consider an iterative likelihood approach wherein the likelihood of an
$n$-component solution is compared against that of an optimized
$(n-1)$-component solution.

##  3 Datasets

###  3.1 Data format

A typical VHE gamma-ray observation consists of a list of gamma-like events
with reconstructed coordinates corresponding to positions in the sky, energy
and time. In what follows we will consider the search for steady emissions and
ignore the temporal dimension, but as it is of prime importance in the case of
VHE transient physics, it should be considered for subsequent work. For this
study, we work with binned events rather than an event list. Thus our training
data, which we henceforth refer to as "observation", corresponds to a 3D
histogram with two dimensions in space and one in energy.

###  3.2 Mock data

####  3.2.1 Toy model mock data

We start by building a simple mock observation with optimal conditions. This
allows us to establish that the models and method are sound for a simple case
and draw comparisons with messier realistic scenarios.

We build a 3D histogram with $50\times 50$ bins in "space" over $[0,1]^{2}$
and 20 bins in "energy" over $[0,1]$. We set $S_{s},S_{b},E_{s}$ and $E_{b}$
as follows:

$\displaystyle\begin{cases}S_{s}(\vec{x})=C_{s}\exp{-(\vec{x}-\vec{x_{0}})^{2}},\\\
S_{b}(\vec{x})=C_{b},\\\ E_{s}(e)=\alpha_{s}e^{2},\\\
E_{b}(e)=\alpha_{b}(1-e),\\\ \end{cases}$ (3.1)

where $C_{s}$ and $C_{b}$ are constant fixed such that
$\int{S_{s}}=\int{S_{b}}$, and such that average particle count per bin is 10,
and $\alpha_{s}$ and $\alpha_{b}$ are normalizing constants such that
$\int{E_{s}}=\int{E_{b}}=1$ on the considered interval.

Figure 4: Observations of the Crab nebula (left panels) and the pulsar wind
nebula MSH 12-52 (right panels) from the H.E.S.S. public release [18]. Top and
bottom panels show the spatial and energy count distributions of the
observations, respectively.

To construct an observation shown in the left-hand side of Fig. 3, we compute
a Poisson realization of $F_{tot}(\vec{x},e)$ for each bin $(\vec{x},e)$.
Looking at the energy spectra as shown on the left-hand side of Fig. 5, we can
note that they are easy to differentiate. The source is also easily
identifiable in space, owing to the high source-to-background
($\rm{StB}=\frac{S}{B}$) ratio (up to 60 in the brightest spatial bin) and,
consequently, the high signal-to-noise ratio ($\rm{StN}=\frac{S}{\sqrt{B}}$)
where the source is present. All these factors (distinctness of shapes, high
StB and StN) are decisive in the ground truth corresponding to the global
maximum of the likelihood (given our model constraints) and the models
consistently reaching this global maximum.

####  3.2.2 Mock dark matter signal in the Galactic Center

We produce mock data of H.E.S.S. observations of the GC region with the five-
telescope array [19]. These observations obtained in the H.E.S.S. Inner Galaxy
Survey provide the highest exposure in the GC region VHE gamma-rays amounting
to 546 hours distributed over the inner few degrees of the Milky Way halo. We
consider a 2-component mixture of a VHE gamma-ray signal induced by the
annihilation of Majorana dark matter particles populating the DM Milky-Way
halo and a residual background component. The energy-differential residual
background flux measured by H.E.S.S. is extracted from Ref. [19], and is
assumed to be spatially isotropic. Such an approximation is justified since no
significant spatial dependence was found down to 1% level [19]. For more
details on the mock data simulation framework, see Ref. [20]. The signal
component will be given by the energy-differential flux of gamma rays in a
solid angle $\Delta\Omega$ expected from the pair-annihilation of Majorana DM
particles of mass $m_{\rm MD}$ and in a DM halo of density $\rho$, expressed
as:

$\displaystyle\frac{d\Phi}{dE}$ $\displaystyle=\frac{\langle\sigma
v\rangle}{8\pi m_{\rm DM}^{2}}\sum_{f}{\rm
BR}_{f}\frac{dN_{f}}{dE}\,J(\Delta\Omega)\quad{\rm with}\quad J(\Delta\Omega)$
$\displaystyle=\int_{\Delta\Omega}d\Omega\int_{\rm
los}ds\,\rho^{2}(s[r,\theta])\,.$ (3.2)

Following the methodology given in Ref. [20], we compute the expected number
of signal events for annihilation of DM of 1 TeV mass in the W+W- channel with
an annihilation cross section of 4$\times$10-25 cm3s-1, and a DM density
profile following the Einasto parametrization with profile parameters
extracted from Ref. [19] in any energy bin $i$ and in any spatial Galactic
longitude bin $j$, and spatial Galactic longitude bin $k$. The time exposure
map is assumed to have a constant value of 500 hours. The energy resolution is
assumed to follow a Gaussian with width given by $\sigma$/E = 10%. The StB in
the brightest spatial bin is equal to 1. The events are binned in 30
logarithmically-spaced bins between 0.1 and 100 TeV (i.e. 0.1 dex.), and 40
bins in Galactic longitude in and latitude between $-10^{\circ}$ and
$+10^{\circ}$.

Figure 5: Estimations for the simple mock case for 10 models run with
different initial conditions. Upper left: Normalized energy counts for source
(red) and background(blue). Upper right: Total counts for bins over longitude
axis. Full lines show ground truth, dotted lines show model results. The
lines, ranging from darkest to lightest, represent estimations from highest to
lowest likelihood. Energy and position are in arbitrary units.Bottom left:
Loss evolution for the 10 runs. Bottom right: Spatial source estimation for
the 10 runs, arranged from highest to lowest likelihood from left to right and
top to bottom.

###  3.3 Observational public data

Finally we aim to test our model on true observational data. We choose two
sources of interest, the Crab Nebula and the pulsar wind nebula MSH 15-52. We
construct our observations by stacking observational runs (i.e. individual
observations of the sources over approximately 28 minutes) obtained from
public H.E.S.S. data [18] and counting the resulting event list in a 3D
histogram. The Crab and MSH 15-52 datasets correspond to observation times of
1.9 (4 runs) and 9.1 (20 runs) hours, respectively. We select a field of view
of 1.5 degrees centered on the source nominal position and consider events in
the energy range from 0.1 to 10 TeV. Within this region, we count events in a
histogram made over 50×50 equally spaced spatial bins and 20 logarithhmically-
spaced energy bins. In this configuration, the StB in the brightest spatial
bin is estimated to be of around 20 for Crab and 0.8 for MSH 15-52.

##  4 Results

###  4.1 Simple mock data

We first show results for the simple mock case (Fig. 3, left). As we can
control important parameters such as total event count and signal-to-noise
ratio when constructing the data, and know the ground truth distributions for
source and background signal, we can accurately determine the performance and
limitations of our network. The model is run ten times with varying randomized
initial weights, for n = 104 gradient updates. We apply only the hat
constraint.

In Fig. 5, we can note that results are both precise and consistent for all
runs, with near-perfect overlap of GT and estimated power spectra (center) and
space functions (right), and loss evolution following a similar pattern for
every run.

Looking closely at the spatial source estimations (bottom), we can note the
networks’ limitation in approximating a 2D gaussian, exhibiting more of a
branched star-like pattern; this is to be expected given our networks’ linear
architecture. However, averaging over a sufficient number of runs should
mitigate this effect.

Figure 6: Weighted average estimations for the simple mock case for hat
constraint. Left: average source space estimation. Middle: averaged energy
functions with standard deviation. Right: averaged spatial functions over
longitude axis with model standard deviation. Energy and position are in
arbitrary units.

We compute a weighted average of the recovered space and energy fields for the
ten runs. For all fields $F$ (corresponding to $S_{s},S_{b},E_{s}$ and
$E_{b}$), $F_{avg}$ is defined as $F_{\rm avg}=\sum_{\rm i}w_{\rm i}F_{\rm i}$
where $w_{\rm i}$ is the weight defined in Eq. (2.5). Correspondingly for all
$F$ we compute a weighted uncertainty $\nu_{F}=\sum_{\rm i}w_{\rm i}|F_{\rm
avg}-F_{\rm i}|$, which can be likened to our model’s epistemic uncertainty.
Looking at the recovered average fields from Fig. 6, we find that the space
gaussian is more closely approximated, while estimations both in space and
energy show near perfect overlap with ground truth, with negligible model
uncertainty.

Defining the weighted relative difference in the following manner:

$D_{R}=\frac{\int|F_{\rm{NN}}-F_{\rm{GT}}|}{\int|F_{\rm{GT}}|},$ (4.1)

the models recover the energy spectra (middle) within $0.4\%$ for source and
$0.3\%$ for background, and the space functions (left/right) within $6\%$ for
source and $1.4\%$ for background.

###  4.2 H.E.S.S. public data

We now compute results for real observational data from the H.E.S.S. public
release. Since we do not have a ground truth for comparison, we compare our
field estimations with those obtained using standard methods (SM); for energy
spectra estimation this corresponds to the MultipleOFF method [9], whereas for
spatial estimation we use the ring background method [9].

####  4.2.1 The point-like Crab nebula

We first observe results for 10 individual runs as shown in Fig. 7, using only
the hat constraint. Looking at the energy spectra (upper-left panel), we can
see that estimations for background spectra show little variation and are
consistent both with one another and with SM. This is to be expected given our
hat constraint, which defines a background-only OFF region, thereby
constraining the background spectrum. Comparatively, we can see more variation
in estimations of the source spectrum; however they remain consistent with SM
measurements.

Figure 7: Estimations for the Crab nebula for 10 models run with different
initial conditions. Upper left: Normalized energy counts in log-bins for
source (red) and background(blue). Upper right: Total counts for bins over
longitude axis. Full lines show standard model results, dotted lines show NN
model results. The lines, ranging from darkest to lightest, represent
estimations from highest to lowest likelihood. Bottom left: Loss evolution for
the 10 runs. Bottom right: Source space estimation for the 10 runs, arranged
from highest to lowest likelihood from left to right and top to bottom.

Looking at spatial source estimations (bottom panel), results seem somewhat
consistent with one another, recovering a Gaussian of slightly varying
amplitude and width, along with surrounding low-amplitude excess. A one-
dimensional spatial observation of source and background (upper-middle panel)
reveals the reason for the varying amplitude of the estimated source, as we
can see the source signal more or less "contaminating" the estimated
background, decreasing the corresponding source estimation in the process.
This can be mitigated by thresholding the background. The spatial NN
estimations are somewhat consistent with SM, though overly negative estimation
of the excess (corresponding to source signal + Poisson fluctuation of the
background) suggests that here SM is overestimating the background in the OFF
region.

The loss functions shown in the upper-right panel, are rescaled
($L^{\prime}=L-c$) to better show their spread (here $c=0.23$). There is very
little difference between final losses, and therefore between the likelihoods
of the final estimations.

As before we now observe the likelihood-weighted average for all fields as
shown in Fig. 8. Looking first at the hat-only constraint case (top), we can
see that the average spatial source estimation (left panel) is smooth in
shape, with the low-amplitude excesses of individual estimations mostly
vanishing. Looking over one spatial dimension (middle panel), we can see that
the average estimated background signal remains apparently contaminated by
source signal, although a flat estimation is within our method’s region of
uncertainty (as defined by weighted standard deviation). This in turn reduces
average estimated source signal, creating a visible inconsistency with SM.
Looking at the averaged energy spectra (right panel) we can see that
estimations are consistent with SM, with a larger uncertainty and disparity
with SM for the source.

Figure 8: Weighted average estimations for the Crab Nebula for the hat-only
(top panels) and hat + bck threshold (bottom panels) constraints,
respectively. Left panel: average source space estimation. Middle panel:
averaged energy functions with standard deviation. Right panel: averaged
spatial functions over longitude axis with Poisson and model standard
deviation.

Upon adding a threshold constraint on the background (bottom panel), we find
that uncertainty is largely decreased for all estimations, with near perfect
agreement with SM for the energy spectra (right panel), and in space (middle)
SM’s estimated excess mostly falling within one Poisson deviation (shown as
yellow-shaded region) of the NN’s average estimated source. This improved
agreement makes sense, given that we are applying similar constraints to that
of the SM.

####  4.2.2 The extended source MSH 15-52

We first observe results for 10 individual runs as shown in Fig. 9, using only
the hat constraint. For the energy spectra (top-left panel), estimations for
background spectra barely show any variation and perfectly overlap with SM
estimation. This is due to the long observation time which provides increased
statistics, allowing for a precise estimation of the background spectrum.

Figure 9: Estimations for the pulsar wind nebula MSH15-52 for 10 models run
with different initial conditions. Upper left: Normalized energy counts in
log-bins for source (red) and background (blue). Upper right: Total counts for
bins over longitude axis. Full lines show standard model results, dotted lines
show NN model results. The lines, ranging from darkest to lightest, represent
estimations from highest to lowest likelihood. Bottom left: Loss evolution for
the 10 runs. Bottom right: Source space estimation for the 10 runs, arranged
from highest to lowest likelihood from left to right and top to bottom.

Here again, while still consistent with SM, the source estimations show more
variation. We can note two spurious cases (yellow and light orange curves),
inconsistent with the other estimations. These correspond to the cases with
highest final losses (upper-right panel), i.e. with lowest final likelihoods.
Looking at the corresponding spatial estimations, it becomes clear that these
correspond to a one-component solution; the estimated source in space is zero
(bottom panel, left-most images), and looking at the one-dimensional case
(upper-middle panel) we can see that the equivalent background estimations
correspond to the total field.

Figure 10: Weighted average estimations for the MSH 15-52 for the hat-only
(top panels) and hat + background threshold (bottom panels) constraints,
respectively. Left panel: average source space estimation. Middle panel:
averaged energy functions with standard deviation. Right: averaged space
functions over longitude axis with Poisson and model standard deviation.

We now observe the likelihood-weighted average for all fields as shown in Fig.
10. In the hat-only constraint case (top panels), the average spatial source
estimation (left panel) is somewhat smooth in shape, with some low-amplitude
excess remaining. More runs might be warranted to smooth these out. Looking
over one spatial dimension (middle panel), the average estimated background
signal is once more contaminated by source signal, especially due to the
previously mentioned one-component estimations affecting the average. A flat
estimation is within our method’s region of uncertainty as defined by weighted
standard deviation. As in the Crab case this reduces average estimated source
signal, creating a slight inconsistency with SM. Still, SM-estimated excess is
essentially within a Poisson deviation of NN-estimated source as shown by the
yellow field.

In the averaged energy spectra (right panel) we can see that background
estimation overlaps perfectly with SM, while source estimation shows more
discrepancy with SM due to one-component estimations. Putting a threshold on
the background (bottom panel) forces a two-component solution, removing the
issues brought up in the previous case; the estimated background is flatter,
all NN estimations are more consistent with SM. However unlike for the Crab
case, source estimations maintain a high uncertainty even with this additional
constraint, likely due to the complex shape and fainter emission of the
source.

###  4.3 Mock diffuse dark matter signal

Figure 11: Estimations for the DM Mock data for 10 models run with different
initial conditions. Upper left: Normalized energy counts in log-bins for
source (red) and background (blue). Upper middle: Total counts for bins over
longitude axis. Full lines show ground truth, dotted lines show model results.
The lines, ranging from darkest to lightest, represent estimations from
highest to lowest likelihood. Bottom left: Loss evolution for the 10 runs.
Bottom right: Source space estimation for the 10 runs, arranged from highest
to lowest likelihood from left to right and top to bottom.

We examine here results for the more complex case of the mock DM signal. Given
the generally low signal-to-background ratio, it is necessary to use both hat
and background threshold constraints to obtain satisfactory results. Since the
source is extended such that it covers the entire field of view, we opt for a
soft hat, where spatial source signal is thresholded to $5\%$ of spatial
background in the OFF region.

We first observe results for 10 individual runs as shown in Fig. 11. For the
energy spectra (top-left panel), estimations for background spectra show once
again little variation and overlap with GT, due to the hat constraint. On the
other hand, source estimations show significantly more variation, with all but
the highest-probability estimations showing significant discrepancy with GT.
These aberrations can be explained by the relatively low total source signal,
which implies low impact on total likelihood estimation even in very erroneous
cases. Inspecting the spatial estimations (upper-middle and bottom panels), we
can note that background estimations are generally identical and consistent
with GT. Conversely, source estimations can be divided into two categories;
essentially flat estimations where source value corresponds to the imposed hat
threshold, and pointed estimations that better capture GT signal shape. We can
note that the second category corresponds to higher likelihood scores, which
will be relevant when computing a weighted average. Observing the losses, we
can see that most first-category (lower-likelihood) cases stagnate throughout
training; this is likely due to the early highly erroneous estimations of
energy spectra preventing a gradual convergence of spatial estimations towards
GT. These cases correspond to local minima.

Now looking at the likelihood-weighted average for all fields shown in Fig.
12, the importance of weighting takes its full measure, as the averaged
estimated source fields suddenly become near-perfectly consistent with GT.
Indeed, we can see that the average spatial source estimation (left panel) is
smooth in shape, that energy functions (middle panel) overlap closely (with an
average relative difference of $7.6\%$ for source and $1.6\%$ for background)
and with little variation for both source and background, and that 1D space
estimations, while showing more significant uncertainty for source, overlap
very well on average with GT, with a relative difference in 2D space of $10\%$
for source and $0.4\%$ for background. A significant finding is that the
model’s prediction of the spectral shape of the signal component accurately
captures the bump-like gamma-ray feature of the $W^{+}W^{-}$ annihilation
channel near the DM mass. This advancement allows for the reconstruction of
DM-induced features in a fully non-parametric manner. It would benefit DM
searches, particularly within realistic DM models where multiple annihilation
channels, weighted by their branching ratios, contribute to the final-state
gamma-ray spectrum, in particular at energies close to the DM mass.

Finally, inspecting the weights (i.e. likelihoods) of individual estimations,
we find that only the top three cases contribute significantly to the weighted
average, while more abhorrent estimations are conferred negligible weight,
allowing for excellent results despite many "failed" runs. This in turn brings
up the necessity of doing a sufficient amount of runs. Indeed, the more
difficult a problem (low StB, low statistics, overlapping spectra, lighter
constraints) the more runs it will take to find and better cover global maxima
for the likelihood.

Figure 12: Weighted average estimations for the DM Mock data, for hat + bkg
threshold constraints. Left: average source space estimation. Middle: averaged
energy functions with standard deviation. Right: averaged space functions over
longitude axis with model (red) standard deviation.

##  5 Discussion

In this paper we developed a NN-based deep ensemble approach for source-
background separation in noisy very-high-energy gamma-ray observations. This
method consists in training $n$ sets of networks to estimate the spatial and
energy components for two fields corresponding to the average source and
background, so as to maximize the likelihood for the observation to be a
realization of their sum (the total estimated field). Constraints are imposed
on the estimations to better differentiate source and background, as well as
to render the task easier. One such constraint (hat) consists in defining and
imposing an OFF region in space where the source is made to be negligible, and
another (bkg threshold) consists in imposing a threshold on the spatial
background intensity, such that brighter regions have to be considered as due
to source or Poisson fluctuation of the background by the NN. The $n$ runs are
then weighted according to their likelihood and averaged to obtain a final
estimation for all components, as well as an uncertainty corresponding to
their weighted standard deviation. We tested this approach both on simple and
more complex realistic mock observations, as well as real observations from
the H.E.S.S. public data release. We compared results to ground truth in mock
cases and results derived from standard methods in real cases. We found that
our approach worked remarkably well in all considered cases, with average
component estimations showing good to near-perfect agreement with ground truth
and standard method, and with performance clearly improving when combining
constraints.

Going forward, it is important to keep in mind a few limitations in our
approach. For one, the model relies on strong hypotheses: that there are only
two signals in the observation, that space and energy components in the
signals are separable, and, depending on the constraints we choose, that we
can confidently define an OFF region in an observation, and that the
observation is "flat" enough (i.e. relatively homogeneous acceptance rate with
respect to space in considered observation) that we can posit a threshold for
the background. These assumptions can reasonably be made for all the cases
considered in this article, however in larger grid-surveys, we can expect more
inhomogeneity in observations to arise due to varying conditions such as
weather, changing azimuthal angles and varying energy acceptance with respect
to observation angle causing space/energy mixing. We can also expect more than
one source per observation, prompting the need for a more generalized
structure that could allow for multiple sources. Additionally, all cases
considered have relatively high ($>0.1$) StB throughout space, with sufficient
exposition that StN is $>1$ throughout the observation’s considered space.
More difficult cases, such as dividing source signal strength by ten for the
mock DM data, did not provide satisfactory results within similar training
time frames and number of runs. Indeed, in such cases, we can expect models to
be more likely to fall into local maxima, while reaching global maxima could
take more time, requiring longer training and more attempts. As this is time-
and energy-consuming, methods to optimize training and better tackle the issue
of local minima (adapting training rate, simulated annealing [21], detecting
stagnation to stop training early, sequential training, etc) should be
considered before addressing such cases.

Further work on this subject could take several paths. Most obvious would
include steps towards generalization, such as adding a temporal dimension, or
allowing for a varying number of sources. Conversely, another path would be to
further make use of physical knowledge for more precise estimations. This
could include taking into account additional instrumental input like space and
energy uncertainty in events, which could help to smooth out Poisson noise, or
information such as instrument response function. Alternatively this approach
could be used in tandem with other analysis methods to better define
constraints or a-priori hypotheses like number of sources or assumed OFF
regions. Finally, experimenting with various architectures, within a set of
deep ensemble runs or for all sets, should be instrumental in optimizing
results. Dense NNs are interesting in the context of a proof of concept, but
one could expect a convolutional neural network (CNN, [22]) taking the
observation as input as well as a base of comparison with the output to
compute the loss (as in Ref. [23]), our model only doing the latter, to be
more suited to the task. One might also explore more direct non-parametric
methods for fitting the observation histogram, such as using wavelets with
smoothness constraints [24].

When aiming to solve the problem of source-background separation in noisy
observations, we sought a method to find and display a solution, or rather a
manifold of solutions, that maximizes the likelihood of the observation while
optimally exploring function space and displaying model uncertainty, as could
be expected in noisy observations. We found that our approach, using a deep-
ensemble set of low-parameter dense nets, proved greatly effective in this
regard. Critical to obtaining optimal and accurate final estimations were two
key factors. The first was the use of informative constraints on source and
background shapes. These constraints helped guide the model towards more
realistic and reliable solutions. The second, and a key finding of this work,
was the use of an adequate weighting method when combining individual results
of the deep ensemble. This method allowed us to automatically discard more
erroneous estimations while enabling similarly plausible estimations to
contribute to the final average. In conclusion, our methodology demonstrates a
robust framework for source-background separation in the presence of noise,
leveraging the strengths of deep ensembles and effective weighting techniques.
This approach not only enhances the accuracy of the results but also provides
a clear representation of model uncertainty, paving the way for more reliable
and insightful analyses in noisy observational data. Future work can build on
these findings to further refine and expand the applicability of our approach
to a broader range of astrophysical problems.

## Acknowledgments

The authors would like to thank Alexander Belikov and Alessandro Montanari, as
well as Mathieu de Bony de Lavergne for interesting and informative
conversations, and Jérôme Bobin for his helpful advice. This project was
funded by the Programme Transversal de Compétences Simulation Numérique of
CEA.

## References

  * [1] https://www.mpi-hd.mpg.de/HESS/.
  * [2] http://www.magic.iac.es.
  * [3] https://veritas.sao.arizona.edu.
  * [4] CTA collaboration, _Science with the Cherenkov Telescope Array_ , WSP (11, 2018), 10.1142/10986, [1709.07997].
  * [5] A. Hillas, _Evolution of ground-based gamma-ray astronomy from the early days to the cherenkov telescope arrays_ , _Astroparticle Physics_ 43 (2013) 19.
  * [6] T.C. Weekes, M.F. Cawley, D. Fegan, K. Gibbs, A. Hillas, P. Kowk et al., _Observation of tev gamma-rays from the crab nebula using the atmospheric cherenkov imaging technique_ , _Astrophysical Journal_ 342 (1989) 379.
  * [7] J. Albert et al., _Implementation of the Random Forest Method for the Imaging Atmospheric Cherenkov Telescope MAGIC_ , _Nucl. Instrum. Meth. A_ 588 (2008) 424 [0709.3719].
  * [8] MAGIC collaboration, _The performance of the MAGIC telescopes using deep convolutional neural networks with CTLearn_ , _PoS_ Gamma2022 (2024) 220 [2211.16009].
  * [9] D. Berge, S. Funk and J. Hinton, _Background Modelling in Very-High-Energy gamma-ray Astronomy_ , _Astron. Astrophys._ 466 (2007) 1219 [astro-ph/0610959].
  * [10] H.E.S.S. collaboration, _The H.E.S.S. Galactic plane survey_ , _Astron. Astrophys._ 612 (2018) A1 [1804.02432].
  * [11] G.P. Rowell, _A new template background estimate for source searching in TeV gamma -ray astronomy_ , _Astrom. Astrophys._ 410 (2003) 389 [astro-ph/0310025].
  * [12] J. Schmidhuber, _Deep learning in neural networks: An overview_ , _Neural networks_ 61 (2015) 85.
  * [13] D. Wang and J. Chen, _Supervised speech separation based on deep learning: An overview_ , _IEEE/ACM transactions on audio, speech, and language processing_ 26 (2018) 1702.
  * [14] W. Zhu, S.M. Mousavi and G.C. Beroza, _Seismic signal denoising and decomposition using deep neural networks_ , _IEEE Transactions on Geoscience and Remote Sensing_ 57 (2019) 9476.
  * [15] J. Xu and H. Liu, _Deep-learning-based separation of a mixture of dual-tracer single-acquisition pet signals with equal half-lives: A simulation study_ , _IEEE Transactions on Radiation and Plasma Medical Sciences_ 3 (2019) 649.
  * [16] J. Aschersleben, T. Arnesen, R. Peletier, M. Vecchi, C. Vlasakidis and M. Wilkinson, _Signal-background separation and energy reconstruction of gamma rays using pattern spectra and convolutional neural networks for the small-sized telescopes of the cherenkov telescope array_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 1059 (2024) 168942\.
  * [17] D.P. Kingma and J. Ba, _Adam: A Method for Stochastic Optimization_ , _arXiv e-prints_ (2014) arXiv:1412.6980 [1412.6980].
  * [18] https://www.mpi-hd.mpg.de/hfm/HESS/.
  * [19] H.E.S.S. collaboration, _Search for Dark Matter Annihilation Signals in the H.E.S.S. Inner Galaxy Survey_ , _Phys. Rev. Lett._ 129 (2022) 111101 [2207.10471].
  * [20] A. Montanari, E. Moulin and N.L. Rodd, _Toward the ultimate reach of current imaging atmospheric Cherenkov telescopes and their sensitivity to TeV dark matter_ , _Phys. Rev. D_ 107 (2023) 043028 [2210.03140].
  * [21] D. Bertsimas and J. Tsitsiklis, _Simulated annealing_ , _Statistical science_ 8 (1993) 10.
  * [22] K. O’shea and R. Nash, _An introduction to convolutional neural networks_ , _arXiv preprint arXiv:1511.08458_ (2015) .
  * [23] F. List, N.L. Rodd, G.F. Lewis and I. Bhat, _Galactic center excess in a new light: Disentangling the $\gamma$-ray sky with bayesian graph convolutional neural networks_, _Phys. Rev. Lett._ 125 (2020) 241102.
  * [24] O. Rioul and M. Vetterli, _Wavelets and signal processing_ , _IEEE signal processing magazine_ 8 (1991) 14.