††institutetext: Department of Physics & National Center for Supercomputing
Applications (NCSA)
University of Illinois at Urbana-Champaign

# A Detailed Study of Interpretability of Deep Neural Network based Top
Taggers

Ayush Khot  Mark S. Neubauer  and Avik Roy111Corresponding author.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Recent developments in the methods of explainable AI (xAI) methods allow us to
explore the inner workings of deep neural networks (DNNs), revealing crucial
information about input-output relationships and realizing how data connects
with machine learning models. In this paper we explore interpretability of DNN
models designed for identifying jets coming from top quark decay in the high
energy proton-proton collisions at the Large Hadron Collider (LHC). We review
a subset of existing such top tagger models and explore different quantitative
methods to identify which features play the most important roles in
identifying the top jets. We also investigate how and why feature importance
varies across different xAI metrics, how feature correlations impact their
explainability, and how latent space representations encode information as
well as correlate with physically meaningful quantities. Our studies uncover
some major pitfalls of existing xAI methods and illustrate how they can be
overcome to obtain consistent and meaningful interpretation of these models.
We additionally illustrate the activity of hidden layers as Neural Activation
Pattern (NAP) diagrams and demonstrate how they can be used to understand how
DNNs relay information across the layers and how this understanding can help
us to make such models significantly simpler by allowing effective model
reoptimization and hyperparameter tuning. While the primary focus of this work
remains a detailed study of interpretability of DNN-based top tagger models,
it also features state-of-the art performance obtained from modified
implementation of existing networks.

##  1 Introduction

Machine learning models are ubiquitous in experimental high energy physics.
With an ever increasing volume of data coupled with complex detector
phenomenology, these models have been extremely useful to find meaningful
information from these very large datasets. Over time, machine learning models
have grown in complexity and simpler regression and classification models have
been replaced by intricate and deep neural networks. Owing to their
intractably large number of trainable parameters and arbitrarily complex non-
linear nature, deep neural networks have been treated as black boxes. It has
always been difficult to understand how different input features contribute to
the network’s computational process and how the inter-connected neural
pathways convey information. In recent years, advances in explainabale
Artificial Intelligence (xAI) MILLER20191  have made it possible to build
intelligible relationship between an AI model’s inputs, architecture, and
predictions xAI-intro ; xAI-review ; xAI-review-2 . While many methods remain
model agnostic, a substantial subset of these methods have been developed to
infer interpretability of computer vision models where an intuitive reasoning
can be extracted from human-annotated datasets to validate xAI techniques.
However, in other data structures, like large tabular data or relational data
constructs like graph, use of xAI methods are still quite novel
sahakyan2021explainable ; xAI-GNN-survey . In recent times, xAI has been
successful in learning the underlying physics of a number of problems in high
energy detectors 9302535 , including parton showers at the LHC LAI2022137055
and jet reconstruction using particle flow algorithms mokhtar2021explaining .

One of the major applications of ML in the field of HEP is classification i.e.
tagging of jets. Jets represent hadronic showers observed as a conical spray
of particles originating from quarks and gluons produced in the high energy
collisions at high energy colliders like the LHC. Identifying jets that
originate from decay products of a certain particle like the top quark ($t$)
and being able to separate them from other jet categories, like jets
originating from the QCD background is an important challenge in many physics
analyses. Traditional top tagging based on kinematic cuts (see Refs.
kaplan2008top ; almeida2009top ; almeida2010template ; plehn2012top  for
example) have been used in particle phenomenology research as well as by ATLAS
and CMS experiments and their predecessors. In Run 1 physics analyses, top
taggers based on simple kinematic cuts and low-complexity statistical models
like decision trees provided satisfactory performance in terms of tagging
accuracy and efficiency aad2016identification ; CMS-PAS-JME-09-001 ; CMS-PAS-
JME-13-007 . However, models based on deep neural networks (DNNs) started
becoming popular in Run 2 at a higher center of mass energy of 13 TeV
atlastopodnn ; cms2020identification . For top quarks produced with large
momenta, the decay products can be confined within a smaller physical cone and
be reconstructed as a single jet. For such boosted jets, top tagging can be
particularly challenging and require a better analysis of jet substructures, a
collection of constituents and their derivative properties that can offer
better discrimination properties. DNNs have been proved to be widely useful to
utilize the jet substructure properties in performing jet classification. A
wide variety of DNN models have been developed to optimize top tagging
pearkes2017jet ; liam2019reports ; datta2017how ; louppe2019qcd ;
butter2018deep ; komiske2019energy ; qu2020jet ; macaluso2018pulling ;
erdmann2019lorentz ; egan2017long ; moreno2020jedi ; gong2022efficient . An
elaborate review and comparison of many of these models has been given in Ref.
kasieczka2019machine . While some of these models have exploited DNN’s
capacity to approximate arbitrary non-linear function hornik1989multilayer
and their huge success with problems in the field of computer vision, other
models have been inspired by the underlying physics and have taken advantage
of physical symmetries and physics-inspired feature engineering. Nevertheless,
exploring interpretability of jet tagger models is quite novel in literature.
Development of ML models for top tagging has been mostly interpretation-
agnostic– exploration of tractable input-output relationship or information
propagation through neural architecture has been largely ignored. Attempts to
develop interpretable, physics-inspired architectures has mostly been
concentrated within clever feature engineering, in creating or augmenting
input feature spaces with physically meaningful quantities erdmann2019lorentz
; chakraborty2019interpretable ; agarwal2021explainable . However, identifying
the relative importance of different input features to establish tractable
input-output relationships or standardize model training paradigm has been
only minimally addressed. As it has been pointed out in Ref.
shanahan2022snowmass , explainability of DNN-based models is important for
validation of such models and to understand how differences in physical
distributions between simulation and real data can lead to differences in a
network’s behavior. The broader scope of uncertainty quantification in
association with ML models relies on developing robust explanations
seuss2021bridging  and in the field of HEP for problems like top tagging will
require dedicated understanding of how robust as well as interpretable these
models are grojean2022lessons .

Yet another remarkable application of interpretability is to understand how
the model conveys information and in doing so, which parts of a DNN most
actively engage in forward propagation of information. Such studies could be
useful to understand model complexity and eventually lead to novel techniques
of model reoptimization. Given DNNs have shown remarkable success in jet and
event classification, recent work has placed emphasis on developing DNN-
enabled FPGAs for trigger-level applications at the LHC duarte2018fast ;
iiyama2021distance ; heintz2020accelerated . As resource consumption and
latency of FPGAs directly depend on the size of the network to be implemented,
it is definitely easier to embed simpler networks on these devices. Hence,
methods that allow interpreting a network’s response patterns as well as
provide critical insights about model optimization without compromising its
performance can greatly benefit these new budding field of ML applications,
especially for online event selection and jet tagging at current and future
high energy colliders.

Application of state-of-the art explainability techniques for interpreting jet
tagger models is receiving more attention recently agarwal2021explainable ;
roy2022interpretability ; neubauer2022explainable  and has been demonstrated
to be successful in identifying feature importance for models like the
Interaction Network IN .

In this paper, we study the interpretability of a subset of existing ML-based
top tagging models. The models we have chosen use multi-layer perceptrons
(MLPs) as underlying neural architecture in order to elucidate the
applicability and limitations of existing xAI methods without convoluting them
with the complexity of larger models with unorthodox data structures. To
compare our results for different models as well as with existing benchmarks
in published literature, we use the dataset developed by the authors of Ref.
butter2018deep  and later used in the top tagger model review in Ref.
kasieczka2019machine . The models explored in this paper along with the
dataset have been reviewed in section 2. The model hyperparameters explained
in this section will constitute the baseline model in each category. Variants
of each model are studied to better understand their interpretability where
the underlying architecture remains same but model hyperparameters, input
features, or data preprocessing might be changed. Section 3 reviews state-of-
the-art xAI methods that we will be using in investigating explainability of
top tagger models. In section 4, we analyze the results of applying xAI
methods on different top tagger models.

##  2 Review of Top Tagging Dataset and Models

The dataset used in this paper has been used to perform model benchmarking
studies in Ref. kasieczka2019machine  and publicly available at Ref. topdata .
This dataset consists of 1 million top (signal) jets and 1 million QCD
(background) jets genetated with Pythia8 pythia8  with its default tune at 14
TeV center of mass energy for proton-proton collisions. The detector
simulation was performed with Delphes delphes  and jets were reconstructed
using the $anti-k_{t}$ algorithm cacciari2008anti  with a jet radius of
$R=0.8$ using FastJet cacciari2012fastjet . Only jets with transverse momenta
within the range of $550$ and $650$ GeV are considered. For each jet, the
dataset contains the four momenta of up to 200 constituents with zero-padded
entries for missing constituents. The dataset is split into traning,
validation, and testing sets with a 6:2:2 split. Some characteristic jet
features from a random subsample of the training data are shown in Figure 1.

(a)

(b)

(c)

Figure 1: Distribution of LABEL:sub@fig:Nconst number of constituent
particles, LABEL:sub@fig:jet-pt jet transverse momentum ($p_{T}$), and
LABEL:sub@fig:jet-m jet mass for background and signal jets

In this paper we consider three different NN-based models for top tagging.
Given the tagger distinguishes between two jet classes, minimizing the
standard binary cross-entropy (BCE) loss has been used as the training
objective for all models. The training is done using the Adam optimizer with
minibatches. All networks showed comparable performance with different
batchsizes. The architecture, hyperparameters, and data preprocessing for each
of the baseline models is summarized below-

  * TopoDNN pearkes2017jet ; atlastopodnn : The simplest top tagging model we consider is a fully connected multi-layer perceptron (MLP) trained with transverse momentum ($p_{T}$), azimuthal angle ($\phi$), and pseudorapidity ($\eta$) of the 30 most energetic particles. Usually referred to as TopoDNN, this model is trained on preprocessed data where (i) the jet is rotated on the $\eta-\phi$ plane to have the most energetic component aligned along the central coordinate (0,0), (ii) the second most energetic component falls along the negative-$\phi$ axis, and (iii) all momenta are scaled with an arbitrarily chosen factor of 1/1700. The transformations (i) and (ii) take advantage of the underlying Lorentz invariance of collider physics and (iii) converts the momenta into unitless quantities and scales them down to a numerical range comparable to those of $\eta,\phi$ quantities. The baseline model is constructed with 4 hidden layers with 300, 102, 12, and 6 nodes respectively and ReLU activation function. The output layer consists of a single node which is converted by the sigmoid function to represent the probability of the jet being classified as a signal jet.

  * Multi-body $N$-subjettiness (MB$N$S) liam2019reports ; datta2017how : Top-tagging with $N$-subjettiness variables uses an MLP as the underlying trainable architecture. However, the input to the network is different from the usual kinematic variables. It uses the multi-body $N$-subjettiness variables thaler2011identifying , defined as

$\tau_{n}^{(\beta)}=\frac{1}{p_{T,J}}\sum_{i}p_{T,i}\min\left\\{\Delta
R_{1i}^{\beta},\Delta R_{2i}^{\beta},...,\Delta R_{ni}^{\beta}\right\\}$ (1)

where $p_{T,J}$ and $p_{T,i}$ represent the transverse momenta of the jet and
its $i$-th constituent and $\Delta R_{ki}$ is the distance between the $k$-th
jet axis and the $i$-th particle constituent. The $n$ jet axes chosen for
claculating $\tau_{n}^{(\beta)}$ are obtained using the $k_{t}$ algorithm
ellis1993successive  with $E$-scheme recombination blazeya2000run . The input
to MB$N$S tagger is the set of subjettiness variables

$\left\\{\tau_{1}^{(0.5)},\tau_{1}^{(1)},\tau_{1}^{(2)},\tau_{2}^{(0.5)},\tau_{2}^{(1)},\tau_{2}^{(2)},...,\tau_{N-2}^{(0.5)},\tau_{N-2}^{(1)},\tau_{N-2}^{(2)},\tau_{N-1}^{(1)},\tau_{N-1}^{(2)}\right\\}\bigcup\left\\{p_{T,J},m_{J}\right\\}$
(2)

where, besides the subjettiness variables, the jet $p_{T}$ and jet mass
$(m_{J})$ variables are used as inputs to provide a kinematic scale for the
jet event. However, the latter inputs are scaled by a factor of 1/1000 to
mitigate the several orders of magnitude gap between their numerical range and
those of the $\tau$s.

(a)

(b)

(c)

(d)

Figure 2: Distribution of LABEL:sub@fig:tau_1_10 $\tau_{1}^{(1)}$,
LABEL:sub@fig:tau_2_10 $\tau_{2}^{(1)}$, LABEL:sub@fig:tau_3_05
$\tau_{3}^{(0.5)}$, and LABEL:sub@fig:tau_5_20 $\tau_{5}^{(2)}$ for background
and signal jets

In our work we consider the MB$8$S model where the MLP consists of 4 hidden
layers with (200, 200, 50, 50) nodes respectively. The ReLU activation
function is used for the hidden layers. The output layer consists of 2 nodes
transformed by the SoftMax function to respectively represent the
probabilities of the jet being classified as a background or signal jet.

  * Particle Flow Network (PFN) komiske2019energy : PFNs are built as infrared- and colinear- (IRC) safe networks based on the deep set zaheer2017deep  architecture. The deep set architecture inherently makes the network invariant under permutation of particle constituents. The model is implements the following relation

$\mathrm{PFN}=F\left(\sum_{i=0}^{N-1}\Phi\left(p_{i}\right)\right)$ (3)

where $F$ and $\Phi$ represent non-linear functions implemented as trainable
NNs, $N$ is the number of jet constituents and $p_{i}$ is the four momentum of
the $i$-th jet constituent. We train the network with the $p_{T},\eta,\phi$ of
jet constituents as input. As a part of data preprocessing, we standardized
the constituents’ $\eta$ and $\phi$ by subtracting the jet’s $\eta$ and
$\phi$. Also, the $p_{T}$ values of the jet constituents are scaled by a the
inverse of sum of constituent $p_{T}$s, i.e.
$\nicefrac{{1}}{{\sum_{i}p_{T,i}}}$. The $\Phi$ network is implemented as an
MLP with 3 layers of 100, 100, and 256 nodes respectively. Each layer is
followed by a ReLU activation layer. The output layer of $\Phi$ represents a
256 dimensional latent space of jet representation. The $F$ network consists
of 3 hidden layers with 100 nodes per layer with ReLU activations. The output
consists of two nodes transformed by the SoftMax operation to represent the
probabilities corresponding to each jet class.

##  3 Interpretability of Machine Learning Models: Tools and Methods

While a number of xAI techniques have been developed, how any one these
methods actually explain an ML model can actually be quite different from the
others wang2009feature ; van2021evaluating ; jesus2021can . Often we find the
xAI techniques producing diverging explanations, making it challenging to rely
on these methods. We investigate a number of these techniques and compare the
corresponding results and try to understand what may contribute to their
divergence. Here, we summarize the methods that we are going to use to explore
the interpretability of top tagger models.

  * The $\Delta\mathrm{AUC}$ method: Identifying feature importance has been an important part of studying classification models tang2014feature . In standard feature selection tasks, a reasonable subset of the features that excels in some model performance metric is chosen. Although it is conceptually different from feature ranking in post-hoc model interpretation, many interpretation metrics also rely on identifying feature importance by optimizing a simpler surrogate model by minimizing a model’s performance loss ribeiro2016should . One of the most useful model analysis tool of a binary classification is the Region Operator Characteristic (ROC) curve, and the Area Under the Curve (AUC) serves as a scalar metric for evaluating model performance. ROC-AUC based feature ranking has been widely promoted in ML literature chen2008fast ; wang2009feature ; serrano2010feature . We adapt those same principles for our model interpretation studies. One straightforward way of evaluating a feature’s contribution in making predictions is to investigate the model’s performance when a particular feature is masked from the input- by replacing it with a population-wide average value or a zero value, whichever is contextually relevant to the model’s relationship with the training dataset.

  * Shapely Additive Explanations (SHAP) lundberg2017unified : The SHAP scores represent a game theoretic approach in identifying the importance of difference features. For each instance of the dataset, the input features are assigned an additive score that determines to what extent much a particular feature contributes to the classifier prediction. An average model prediction is determined by replacing each feature by its population average and then individual features are added back to the model to find their impact on leading the prediction towards the optimal value. For each feature, the SHAP score is determined by evaluating the average contribution of adding the feature over all possible feature subsets defined without that feature. Given that evaluating exact SHAP scores require iterating over $2^{n}$ sets of feature combinations for each data instance with $n$ features, several simplifying assumptions are made to reduce the computational complexity. In our work, we use the kernel SHAP method- a model agnostic approach to obtain local explanations similar to the LIME framework ribeiro2016model  and obtained by generating random samples around the data point and performing a mean-squared-error-minimizing linear regression over the samples to evaluate the SHAP scores. In order to avoid any overfitting in obtaining the SHAP score, the number of samples in each model were chosen to be at least twice as many as the number of input features

  * Layerwise Relevance Propagation (LRP) LRP-NN ; LRP-overview : The LRP technique propagates the classification score predicted by the network backwards through the layers of the network and attributes a partial relevance score to each input. The backpropagation of LRP scores in an MLP network is obtained by the following relation-

$r_{j}^{(n)}=\sum_{k}\frac{a_{j}^{(n)}w_{jk:n}}{\sum_{j}a_{j}^{(n)}w_{jk:n}}r_{k}^{(n+1)}$
(4)

where $a_{j}^{(n)}$ and $r_{j}^{(n)}$ are the activation and relevance scores
of the $j$-th node in $n$-th layer and $w_{jk:n}$ is the weight that
determines the contribution of the $j$-th activation in the $n$-th layer to
the $k$-th node in layer $n+1$. The inputs to the network are identified as
the $0$-th layer and the relevance scores assigned to them are denoted as
$r_{j}^{(0)}$. The original LRP method has been developed for simple MLP
networks. Variants of this method have been explored to propagate relevance
across convolutional neural networks LRP-pixel  and graph neural networks
schnake2021higher . While the basic LRP rule in Eqn. 4 conserves the total
relevance score i.e. the classifier network’s output, based on the
distribution of weights and activations, relevance scores can become unbounded
when $\sum_{k}a_{j}^{(n)}w_{jk:n}\to 0$. To overcome this, the LRP rule is
modified to treat postive and negative weights asymmetrically. We use the so
called LRP-$\gamma$ rule defined as-

$r_{j}^{(n)}=\sum_{k}\frac{a_{j}^{(n)}\left(w_{jk:n}+\gamma
w_{jk:n}^{+}\right)}{\sum_{j}a_{j}^{(n)}\left(w_{jk:n}+\gamma
w_{jk:n}^{+}\right)}r_{k}^{(n+1)}$ (5)

where $w^{+}=w\cdot\Theta(w)$, $\Theta$ being the Heaviside step function and
$\gamma$ is a regularization parameter. In our representation, we always
present the normalized relevance scores so that $\sum_{j}r_{j}^{(0)}=1$.

  * Neural Activation Pattern (NAP) Diagrams roy2022interpretability ; neubauer2022explainable : While the aforementioned methods help identify the importance of features for a trained NN model, the NAP diagrmas visualize the information propagation pathways through the network’s architecture. The NAP diagram visualizes the Relative Neural Activation (RNA) score, defined as-

$\mathrm{RNA}(j,k;\mathcal{S})=\frac{\sum_{i=1}^{N}a_{j,k}(s_{i})}{\max_{j}\sum_{i=1}^{N}a_{j,k}(s_{i})}$
(6)

where $\mathcal{S}=\\{s_{1},s_{2}...s_{N}\\}$ represents a set of samples over
which the $\mathrm{RNA}$ score is evaluated. The quantity $a_{j,k}(s_{i})$ is
the activation of $j$-th neuron in the $k$-th layer when the input to the
network is $s_{i}$. When summed over all the samples in the evaluation set
$\mathcal{S}$, this represents the cumulative neural response of a node, which
is normalized with respect to the largest cumulative neural response in the
same layer to obtain the $\mathrm{RNA}$ score. Hence, in each layer, there
will be at least one node with an $\mathrm{RNA}$ score of 1. Since the neurons
are activated with ReLU activation in the models we consider, the
$\mathrm{RNA}$ score will be strictly non-negative, and $\leq 1$. In a
qualitative way, we are trying to see which neurons most actively engage to
obtain the predictions from our models. Since the MLPs in our models consist
of only Dense layers, each layer takes all the activations from the previous
layer as inputs. As all nodes within a given layer are subject to the same set
of inputs, we can reliably estimate how strongly they perceive and transfer
that information to the next layer by looking at their activation values. For
the same reason, we normalize the cumulative activation of a node with respect
to the largest aggregate in the same layer.

The NAP diagram is obtained by presenting the RNA scores for the different
layers of the model as a two dimensional heatmap where along the horizontal
axis lies the different activation layers of the network and the vertical axis
represents the different nodes in those activation layers. NAP diagrams
illustrate the relative activity level of different nodes within each layer
and hence can demonstrate the sparsity of the model’s activity.

##  4 Model Interpretability for Top Taggers

###  4.1 TopoDNN

Since TopoDNN is arguably the simplest NN-based model to perform top tagging,
it is perhaps the most ideal model to investigate different aspects of xAI.
Given that correlation among features has been demonstrated to be an important
aspect of identifying feature rankings in classical machine learning
tolocsi2011classification  as well as modern xAI methods, we start by
examining the pairwise Pearson correlation coefficient for a subset of the
input features for background and signal jets in Figure 3. The correlation
matrices for both jet categories are mostly sparse except for some large anti-
correlations between pt_0, the transverse momentum of the most energetic jet
constituent and that of some of the low energy constituents. Given the dataset
has been generated within a limited $p_{T}$ range, such anti-correlations are
expected- the higher the energy of the most energetic constituent, the lower
the energy of the remaining constituents. Given that the numerical range of
the $p_{T}$ of lower energy constituents is typically much smaller than that
of the highest energy constituent, we can expect the impact of their anti-
correlations with the pt_0 on the NN’s performance to be rather small. We can
indeed verify that in Figures 4a–4c where we identify the important features
for the TopoDNN model using the $\Delta\mathrm{AUC}$ (Figure 4a) and SHAP
scores (Figures 4b and 4c).

(a)

(b)

Figure 3: Feature correlation matrices for $p_{T},\eta,$ and $\phi$ of the 10
highest energy constituents of LABEL:sub@fig:featTopoDNNbkg QCD background and
LABEL:sub@fig:featTopoDNNbkg top jets.

The $\Delta\mathrm{AUC}$ score cannot independently identify the features that
contribute to identification of signal and background jets. However, by
evaluating the SHAP scores for subsets of the dataset that only contain one
kind of jets, one can identify the features that most dominantly contribute to
identification of the corresponding jet class. However, as shown in Figures 4b
and 4c, there is a significant overlap between the features that are
identified as the most important ones for both jet categories. Unlike computer
vision models that most deal with image or videos as input data and importance
distribution for different images can vary based on which pixels carry the
most relevant information, the same feature can contribute equally importantly
for different classes in models with tabular data. Why the network treats the
same set of variables as important becomes obvious upon inspecting the
distribution of some of these features as shown in Figures 4d–4g. phi_1,
phi_2, eta_2, all show strong classification characteristics and given these
variables are either loosely correlated or almost uncorrelated as shown in
Figure 3, they all can independently contribute to the network’s ability to
tell apart the different jet classes. On the other hand, pt_0 by itself is a
modest discriminator, and hence identified as having a modest impact on the
model’s performance. This has also been verified by training a variant of the
TopoDNN model that excludes the $p_{T,0}$ variable and as shown in Table A1,
performs almost equally as well as the baseline model.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 4: Feature rankings obtained from the LABEL:sub@fig:dAUCTopoDNN
$\Delta\mathrm{AUC}$ scores, LABEL:sub@fig:bkgSHAPTopoDNN SHAP scores for QCD
background, and LABEL:sub@fig:sigSHAPTopoDNN SHAP scores for top jets. Figures
LABEL:sub@fig:phi1TopoDNN, LABEL:sub@fig:phi2TopoDNN,
LABEL:sub@fig:eta2TopoDNN, LABEL:sub@fig:pt0TopoDNN show the distribution of
phi_1, phi_2, eta_2, and pt_0 for two jet categories

However, we see a stark difference in the distribution of the relevance scores
(Figure 5) among different features, obtained from the LRP method, when
compared to other feature ranking metrics we have considered so far. Unlike
SHAP or $\Delta\mathrm{AUC}$ scores, a subset of the $p_{T}$ variables have
the largest relevance scores. While the features mostly highly ranked in the
previous two methods show strong discriminating characterisitcs, some of the
highly ranked features from the LRP method show very little discriminating
capacity. This difference can be understood from the nature of these ranking
methods. Both $\Delta\mathrm{AUC}$ and SHAP calculate the model’s deviation
from the mean behavior, i.e. qualitatively, they both represent how much
information is obtained from inclusion of the true value of a feature instead
of using the population mean as a feature mask. On the other hand, LRP
calculates the feature’s cumulative relevance, which additively includes the
relevance scores attributed to each feature’s mean behavior. Assuming
$\vec{x}=\\{x_{i}\\}$ be a sample jet event taken from the set of events $X$
and $\vec{{x}}_{\backslash
k}=\vec{x}\backslash\\{x_{k}\\}\bigcup\\{\mathbf{E}\left(X_{k}\right)\\}$ be
the event set where we mask the $k$-th input feature by replacing it with its
mean value, the linear order behavior of the NN can be approximated using the
Deep Taylor Decomposition formalism montavon2017explaining -

$f(\vec{x})\approx f(\vec{{x}}_{\backslash k})+\frac{\partial f}{\partial
x_{k}}\left(x_{k}-\bar{x}_{k}\right)$ (7)

where $f(\vec{x})$ represents the output of the NN before the final Sigmoid
activation. While we recognize that Eqn. 7 is an approximation, especially
given that a random sample $\vec{x}$ can be quite far away from
$\vec{{x}}_{\backslash k}$ in the Euclidean vector space spanned by $X$ so
that nonlinear effects might become important, we only intend to use this
relation as a diagnostic tool to inspect the linear order effect observed on
relevance scores when the model’s inputs are varied away from the population
mean. Noting that the relevance scores additively distribute the functional
output among the different inputs, i.e. $f(\vec{x})=\sum_{i}r(x_{i})$ and
$f(\vec{{x}}_{\backslash k})=\sum_{i\neq k}r(x_{i})+r(\bar{x}_{k})$, we can
rewrite Eqn. 7 as,

$\sum_{i}r(x_{i})\approx\sum_{i\neq k}r(x_{i})+r(\bar{x}_{k})+\frac{\partial
f}{\partial x_{k}}\left(x_{k}-\bar{x}_{k}\right)$ (8)

Collecting terms with equivalent indices and ignoring higher order effects, we
can write $r(x_{k})\approx r(\bar{x}_{k})+\delta r_{k}$ where we denote
$r(\bar{x}_{k})$ as the mean-behavior relevance score and $\delta
r_{k}=f(\vec{x})-f(\vec{{x}}_{\backslash k})\approx\frac{\partial f}{\partial
x_{k}}\left(x_{k}-\bar{x}_{k}\right)$ is identified as the differential
relevance score. Figure 6a show the absolute mean behavior relevance scores of
different features and the relative size and distribution of the relevance
scores are very similar to what we can see in Figure 5. This explains that a
large contribution of the LRP scores actually come from the mean behavior
relevances, and has very little to do with the network’s ability to
distinguish different jet types. In fact, many of the angular variables that
are regarded as highly important by $\Delta\mathrm{AUC}$ and SHAP methods also
show large differential relevance, as shown in the distributions of Mean
Absolute Differential Relevance (MAD Relevance) scores (normalized with
respect to relevance scores from the baseline model) in Figures 6b and 6c.

(a)

(b)

Figure 5: Distribution of the relevance scores obtained from LRP for
LABEL:sub@fig:bkgLRPTopoDNN QCD background and LABEL:sub@fig:sigLRPTopoDNN top
jets.

(a)

(b)

(c)

Figure 6: Distribution of the LABEL:sub@fig:avgLRPTopoDNN absolute mean
behavior relevance scores, LABEL:sub@fig:bkgdLRPTopoDNN MAD relevance score
for QCD background, and LABEL:sub@fig:sigdLRPTopoDNN MAD relevance score top
jets.

Now we turn our focus to examine the behavior of the internal architecture of
the model with NAP diagrams. As discussed in Section 3, NAP diagrams plot the
RNA scores of different nodes of the activation layers within the network.
Figure 7 shows the 2D map of RNA scores for QCD and top jets where the RNA
scores of the former are plotted as negative values to allow simultaneous
visualization. It can be readily understood that that the network in Figure 7
is quite sparse. We can heuristically quantify sparsity of the network by the
fractional number of hidden activation nodes with an RNA score less than a
given threshold. We arbitrarily choose this threshold to be 0.2 and find that
the network’s sparsity measures for background and signal jet categories are
0.86 and 0.76 respectively, giving an overall sparsity measure of 0.70. This
implies that about 70% of the nodes show a cumulative activity level of less
than 20% compared to that of the most active node in the corresponding layer.
We also see in Figure 7 that the most active nodes for different jet
categories are almost completely disentangled by the time information
propagates to the third layer. Large sparsity of the network and early
disentanglement of jet categories indicate the network’s complexity can be
reduced without any noticeable compromise in its performance. To demonstrate
this, we have trained variants of the TopoDNN model with lesser complexity by
simultaneously reducing the depth and width of the MLP model. As shown in
Table A1, these simplified models perform almost equally well while the model
complexity is significantly reduced.

Figure 7: NAP diagram for TopoDNN model visualizing a 2D map of RNA score for
different nodes of the activation layers. To simultaneously visualize the
scores for QCD and top jets, we project the RNA scores of the former as
negative values.

###  4.2 MB$8$S

Although the underlying architecture of the MB$8$S model is an MLP, there is
stark difference among the input features. Unlike the input to the TopoDNN
model explored in Section 4.1, the inputs to the MB$8$S model are highly
correlated for both jet categories (Figure 8). Such large correlations among
features can make it hard to distinguish whether a feature truly conveys
independent discriminating characteristics or if a feature is deemed important
by a model simply because it is correlated with another feature. Training a NN
with correlated feature inputs can contribute to increased model complexity
ayinde2019regularizing , overfitting cogswell2015reducing , and obscure its
interpretability kaur2020interpreting . The MB$8$S network has been trained
with DropOut layers srivastava2014dropout  with a dropout rate of 0.2 (0.1)
for the first (final) two hidden layers to protect it from the problem of
overfitting. However, having large correlations among features can make it
harder to interpret and identify the relative ranking of the input features.
Figure 9a shows the distribution of $\Delta\mathrm{AUC}$ score for the top 10
features. When compared with the SHAP score distributions for background and
signal jets in Figures 9b and 9c, the sets of top ranking features for these
different methods show significant overlap but their raking shows some
noticeable difference, especially for the jet mass feature. As shown in Figure
1c, the distribution of jet mass is very different for the two types of jets
and naturally, it is expected to be a strong discriminator. This is also
reflected in the distribution of SHAP scores. But the $\Delta\mathrm{AUC}$
ranking places the jet mass variable at a lower ranking compared to
subjettiness variables $\tau_{1}^{(2)}$ and $\tau_{2}^{(2)}$. As shown in
Figure 8, these variables demonstrate almost 100% correlation with
$\tau_{1}^{(1)}$ and $\tau_{2}^{(1)}$ respectively, both of which are
identified as top ranking variables in both $\Delta\mathrm{AUC}$ and SHAP
rankings and have strong discriminative distributions as shown in Figures 2a
and 2b. We do not expect a one-to-one correspondence between the feature
rankings from $\Delta\mathrm{AUC}$ and SHAP scores since there is no
established, trivial translation between these scores. Nevertheless, having a
lower rank for jet mass compared to variables that display almost perfect
correlations with other strong discriminators naturally intrigues the question
whether it is overshadowing the importance of variables that actually adds new
information to the classifier. It has been observed that inclusion of jet mass
significantly improves the performance of the MB$N$S taggers liam2019reports ,
it is only expected that jet mass should rank higher.

(a)

(b)

Figure 8: Feature correlation matrices LABEL:sub@fig:featTopoDNNbkg QCD
background and LABEL:sub@fig:featTopoDNNbkg top jets for the different input
features for the MB$8$S network

(a)

(b)

(c)

Figure 9: Feature rankings obtained from the LABEL:sub@fig:dAUCMB8S
$\Delta\mathrm{AUC}$ scores, LABEL:sub@fig:bkgSHAPMB8S SHAP scores for QCD
background, and LABEL:sub@fig:sigSHAPMB8S SHAP scores for top jets for the
MB$8$S network. Only the top 10 features are shown.

In order to disentangle the effect of strong feature correlations, we train a
variant of the MB$8$S network where only the $\tau_{x}^{(1)}$ variables are
included, along with the mass and transverse momentum of jets. This choice is
inspired from the block-diagonal large concentration of pairwise correlations
in Figure 8. This alternate network has almost identical performance as
compared to the baseline MB$8$S network with an AUC score of 97.55% while the
former has an AUC score of 97.98%. The feature rankings via
$\Delta\mathrm{AUC}$ and SHAP scores for this network are shown in Figure 10.
We still see that the rank of jet mass is different between
$\Delta\mathrm{AUC}$ and SHAP score distributions. However, having removed
highly correlated features from the network’s inputs , we can safely assume
these differences are primarily methodological. The top ranking variables in
this network display relatively weaker correlations among themselves
(correlation coefficients $\leq 0.4$) and hence, correctly attributes high
importance to features that contribute new information to the classifier’s
decision making process. Moreover, since the two networks demonstrate
equivalent performance, the highly correlated subjettiness variables only
marginally impact the network’s performance and hence, the relatively high
$\Delta\mathrm{AUC}$ score attributed to $\tau_{1}^{(2)}$ and $\tau_{2}^{(2)}$
in Figure 9a must be an artifact of strong feature correlations.

(a)

(b)

(c)

Figure 10: Feature rankings obtained from the LABEL:sub@fig:dAUCMB8S_taux1
$\Delta\mathrm{AUC}$ scores, LABEL:sub@fig:bkgSHAPMB8S_taux1 SHAP scores for
QCD background, and LABEL:sub@fig:sigSHAPMB8S_taux1 SHAP scores for top jets
for the MB$8$S network trained with
$\\{\tau_{x}^{(1)}\\}\bigcup\\{p_{T,J},m_{J}\\}$ variables

Next we turn our attention to relevance scores attributed to different input
features by the LRP method. Figures 11a and 11c show the mean absolute
relevance scores attributed to different features for background and signal
jets. The feature ranking from relevance scores attributed for background jets
is very different from what we find from previously explored feature ranking
methods. Jet $p_{T}$ is assigned the largest relevance score for background
jets. However, this variable is one of the least expressive features for the
network. It has almost no correlation with other features, and has a very
similar distribution for both jet types. Hence, assigning this feature a very
large relevance score definitely raises some concerns about reliability of the
LRP interpretation of this network. However, distribution of MAD relevance
score in Figure 11b gives a more appropriate distribution for feature
importances. For the top jets, the relevance distribution identifies
$\tau_{1}^{(1)},\tau_{2}^{(1)},m_{J}$ as the most important features, similar
to what we have found from SHAP scores. However, in the distribution of MAD
relevance score (Figure 11d) $\tau_{1}^{(2)}$ ranks higher. Also,
$\tau_{2}^{(2)}$ and $\tau_{1}^{(0.5)}$ have larger MAD relevance scores
compared to $\tau_{2}^{(1)}$. This ambiguity can be attributed to the very
large correlations among these features.

(a)

(b)

(c)

(d)

Figure 11: Distribution of LABEL:sub@fig:LRPMB8Sbkg relevance scores for QCD
background, LABEL:sub@fig:dLRPMB8Sbkg MAD relevance scores for QCD background,
LABEL:sub@fig:LRPMB8Ssig relevance scores for signal top jets, and
LABEL:sub@fig:dLRPMB8Ssig MAD relevance scores for signal top jets for the
MB$8$S network

To demonstrate how the different hidden activation layers contribute to
information propagation for the two jet categories, we show the 2D map of the
RNA scores of different nodes in Figure 12. Although the network appears to be
relatively sparse for different jet categories, with a sparsity measure of
0.74 and 0.64 for background and signal jets measured with respect to a
threshold of 0.2 on RNA scores, different nodes are most strongly activated
for signal and jet categories. As a result, the network’s overall sparsity
becomes 0.44. However, we can already see that the nodes that are most
strongly activated by the two jet categories are almost completely
disentangled at Layer 3. This indicates that the network might be simplified
by choosing a shallower network and indeed verified in Table A1 where we see
that a model trained without the fourth layer has identical performance
metrics when compared with the baseline.

Figure 12: NAP diagram for MB$8$S model visualizing a 2D map of RNA score for
different nodes of the activation layers. To simultaneously visualize the
scores for QCD and top jets, we project the RNA scores of the former as
negative values.

###  4.3 PFN

While the previously analyzed TopoDNN and MB$8$S models both employed a single
MLP to perform the jet classification, PFN utilizes a deep set topology that
linearly combines particle-level neural embeddings to obtain a jet level
latent space, which eventually is used to train a second MLP to learn the
classification task. Befor going into details about interpretability of PFN,
we would like to point out that our implementation of PFN obtained
significantly better performance benchmark compared to what is reported in
Ref. kasieczka2019machine  where the different models reported accuracies
between 0.89 and 0.93 with ROC-AUC score ranging between 0.955 and 0.985. We
obtained an accuracy of 0.977 and an ROC-AUC score of 0.997, which also
outperforms the Interaction Network based model in Ref. moreno2020jedi  and
the recent Graph Neural Network based LorentzNet model gong2022efficient . The
major difference comes from the model’s inputs. While the PFN model trained in
Refs kasieczka2019machine  directly used the constituents’ four momenta as
input to the, we trained the network with $p_{T},\eta,\phi$ variables of the
constituents. The constituent $\eta$ and $\phi$ values were standardized
according to the prescription in Ref. komiske2019energy  where
$(\eta_{i},\phi_{i},p_{T,i})\to(\eta_{i}-\eta_{J},\phi_{i}-\phi_{J},\frac{p_{T,i}}{\sum_{i}p_{T,i}})$.

(a)

(b)

(c)

Figure 13: Feature rankings obtained from the LABEL:sub@fig:dAUCptetaphiPFN
$\Delta\mathrm{AUC}$ scores, LABEL:sub@fig:bkgdLRPptetaphiPFN MAD relevance
scores for QCD background, and LABEL:sub@fig:sigdLRPptetaphiPFN MAD relevance
scores for top jets for the PFN network. Only the top 20 features are shown.

We start by examining the feature ranking from the $\Delta\mathrm{AUC}$ and
MAD relevance scores in Figure 13. The rankings obtained from both methods
show a consistent distribution with the azimuthal distribution of the higher
energy constituents being identified as some of the most important features,
along with the transverse momenta of the two most energetic constituents.
Although the TopoDNN network is also trained on a similar set of input and
there is some overlap in the set of important features for these two networks,
there are two major differences that prohibit one-to-one comparison between
the feature rankings for these variables. Firstly, the input data
preprocessing is very different for these two networks which may result into
significant differences on how these variables are treated by the
corresponding models. And secondly, particle level information fed to PFN is
directly used to obtain particle level neural embeddings while they are
directly used to obtain jet classification for the TopoDNN network. Hence, the
importance metric attributed to each component in PFN is indicative of how
impactful that constituent is in obtaining the jet level embedding while the
relative ranking of the features of a given constituent suggests how important
those features are for the network to obtain the constituent level embedding.
For instance, the variables $\phi$ and $p_{T}$ for any given constituent
outrank the corresponding $\eta$ which suggests that the $\Phi$ network
assigns more importance to the former in order to obtaian its neural
embedding.

(a)

(b)

Figure 14: NAP diagrams for LABEL:sub@fig:NAPPHIPFN $\Phi$ and
LABEL:sub@fig:NAPFCPFN $F$ networks for the PFN model

The particle level embeddings obtained by the $\Phi$ network are summed over
to obtain a latent space representation of the jet. Characterization of latent
spaces has been a topic of general interest in many areas of machine learning
application. For instance, disentangling semantic features of images via
latent spaces in Variational AutoEncoders (VAEs) kingma2013auto  and its
variants burgess2018understanding ; hadjeres2017glsr ; bajaj2021invariance
has been widely studied in modern machine learning literature. In the context
of collider physics, how latent spaces embed information and can be used as
effective candidates for anomaly detection and bump hunting has been studied
moreno2020jedi ; komiske2019energy ; bortolato2022bump . The proponents of PFN
performed detailed studies showing how the latent space representation forms
discernible contours in the $(\eta,\phi)$ plane of jet image representations.
While such studies are useful to divulge geometric features of the latent
space configuration, interpreting how they actually contribute to the
network’s decision making process, especially for large latent spaces with
$\mathcal{O}(100)$-dimensions, remains unexplored.

Not unlike most NN-based implementations, the latent space of jet
representation for PFN is obtained by training a NN to perform a non-trivial,
nonlinear transformation of the input feature space. The network is trained to
obtain classification scores in two disentangled dimensions in the very last
layer of $F$ network, it is only expected that the information propagation
pathways for the different types of jets will show some level of
disentanglement within the hidden layers of the networks. This is indeed
verified by the NAP diagrams shown in Figure 14. These NAP diagrams reveal
some crucial insights. We clearly see that the activity level of different
nodes in the final layer of the $\Phi$ network for background and signal jets
is very similar. It implies that the network embeds the jet level information
in the same latent subspace. Also, the latent space appears to be very sparse
and we indeed found that many of the latent space variables are identically
zero for all events in both jet categories. We retrained variant networks with
latent space dimensions of 128 and 64 and still got comparable performance
(Table A1). A third observation is that the $F$ network effectively learns to
disentangle the representation of jet classes by the second hidden layer
(Figure 14b). This implies that the network can be simplified by reducing the
number of hidden layers. Building on top these observations, we trained
alternate models with modified number of nodes in different hidden layers of
the two networks and found that the network’s accuracy is only marginally
compromised even with much smaller networks.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

Figure 15: Latent space feature rankings obtained from the
LABEL:sub@fig:dAUCPFNz $\Delta\mathrm{AUC}$ scores, LABEL:sub@fig:bkgdLRPPFNz
MAD relevance scores for QCD background, and LABEL:sub@fig:sigdLRPPFNz MAD
relevance scores for top jets. Figures in the bottom two rows show the
distribution of some of the top ranking latent space features.

Given the sparsity of the latent space representation, we expect that only a
small subset of these latent features will actually have a strong contribution
towards the decision making process of $F$ in the baseline model. The ranking
of different features using the $\Delta\mathrm{AUC}$ score and MAD relevance
scores are shown in Figures 15a–15c. There are noticeable overlaps among the
features that rank high with these methods, though the actual sequence of
latent variables, understandably, show some differences. We show the actual
distributions of some of these latent space embeddings in Figures 15d–15k.
These distributions highlight a stark contrast between the the latent space
representation for PFN when compared to those studied in the context of
generative models like VAEs. The latter class of models is especially popular
in computer vision as their latent spaces often provide semantic
disentanglement. Such disentanglement creates a clear separation in latent
space dimensions that account for variabilities in input distributions. For
instance, VAEs trained on the popular celebrity protrait dataset CelebA
liu2015faceattributes  have successfully demonstrated disentanglement of
semantics like age, gender, skin tone, hairstyle etc. burgess2018understanding
. In the context of out top tagging dataset, it is the jet kinematics that is
embedded in the latent space and such embedding alone is often insufficent for
classification dillon2021better . In fact, for models like PFN that are
trained to learn the jet classes and not the distribution of training dataset,
the model has no additional incentive in disentangling the jet classes.
Rather, PFN learns to embed the information regarding jet classes in
correlations among latent features. This can be seen from the latent space
correlation matrices of the two jet classes shown in Figure 16. The pairwise
feature correlations are very different for the two jet classes, creating a
clearer context for the classifier network $F$ to obtain the desired jet
classification.

(a)

(b)

Figure 16: Correlation matrices for latent space features for
LABEL:sub@fig:FeatcorrbkgzPFN background QCD and LABEL:sub@fig:FeatcorrsigzPFN
top jets. Only the top 30 latent features obtained from the
$\Delta\mathrm{AUC}$ metric are shown.

We can further illustrate this by taking a deeper look into the correlations
among the latent features. Since the size of the latent space is large but
sparse, we select the top ranking subset of latent features so that
simultaneously masking each latent feature in the remaining subset causes at
most 1% drop in the AUC score from the test data. We found that for our
baseline PFN model, this requires choosing 85 of the 256 latent features. One
of the classical data analysis tools to identify how features contribute to
variability in data is Principal Component Analysis (PCA)
jolliffe2016principal . PCA performs a linear transformation on these features
to obtain a set of orthogonal feature space with no cross-correlation among
the transformed features. The reduced feature space preserves the sample
variance of the dataset but redistributes them along orthogonal basis vectors.
We found that 99% of sample variance of the test data was described by the top
30 principal components. We show the distribution of some of these components
in Figures 17a\- 17d. We can readily see how these PCA-transformed latent
features can visibly differentiate between the two jet classes in Figures
17e\- 17g.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 17: Distribution of some of the principal components for signal and
background jets (top row) and pairwise distribution between these components
for these jet classes (bottom row)

Strong disentanglement between the jet classes by the principal components of
the latent space demonstrates how the jet class information is actually
embedded in the correlation among different dimensions of the latent space. It
is also instructive to investigate the physical nature of the latent space
learned by the network. While it is neither trivial nor obvious for neural
networks to learn about features that bear resemblance with features that bear
meaning to humans, it has been seen that latent space networks can
occasionally learn about physical variables moreno2020jedi . In case of the
PFN, since the latent space dimensions are highly correlated with each other,
we chose to study the correlation between the principal components and jet
features like jet mass and the number of constituents which, as shown in
Figure 1, can have moderate to strong discriminative feature. We found that
the first principal component, $z_{pca,0}$ (Figure 17a) shows a strong
correlation with jet mass for both jet categories with correlation
coefficients being 0.91 and 0.62 for background and signal jets respectively.
This is particularly impressive considering that the network was trained with
only information regarding the 3-momentum of the constituents and information
about the energy or jet mass was absent from the set of inputs. $z_{pca,0}$
also shows strong correlations with the number of jet constituents, while
$z_{pca,5}$, $z_{pca,11}$ also appears to demonstrate moderate correlations
with constituent multiplicity.

##  5 Conclusion

This paper presents a comprehensive study of the interpretability of DNN based
top tagger models. Our work has unveiled a number of important aspects
regarding how these models connect with the corresponding datasets. We have
observed intriguing inconsistencies in features ranking from different ranking
methods, especially how the LRP method can produce a relevance distribution
that can lead to a misleading interpretation of feature importance. Modifying
LRP to obtain differential relevance scores has been found to be more
consistent with other approaches in xAI. Furthermore, explainability metrics
need to be carefully studied and understood for models trained with highly
correlated input features since, as we show with the MB$8$S model,
interpretability of AI models can be obscured in such cases. Our investigation
suggests that models learn to embed jet class information in correlations
among latent space embeddings and can learn to mimic distributions that
closely resemble physical jet features. On the other hand, RNA scores and NAP
diagrams can lead to effective understanding how information is propagated
through different layers of a network and can lead to efficient model
reoptimization strategies. This opens a possibility of incorporating these
methods to obtain in-situ model optimization during training.

## Acknowledgements

We would like to thank Volodymyr Kindratenko and the Center for Artificial
Intelligence Innovation at the NCSA for support. We are grateful to Gregor
Kasieczka and Huilin Qu for pointing us to important resources and code
repositories. A.K.,A.R.,M.S.N: This work was supported by the FAIR Data
program of the U.S. Department of Energy, Office of Science, Advanced
Scientific Computing Research, under contract number DE-SC0021258.

##  Appendix A Performance of Baseline and variant models

Model Architecture | Model Description | ROC-AUC score | Accuracy | Spartsity
---|---|---|---|---
TopoDNN | Baseline | 0.971 | 0.916 | 0.70
Trained without $p_{T,0}$ | 0.970 | 0.914 | 0.85
Hidden Layers: $(240,80,10)$ | 0.972 | 0.916 | 0.82
Hidden Layers: $(120,40,6)$ | 0.972 | 0.916 | 0.59
MB$8$S | Baseline | 0.980 | 0.928 | 0.44
Trained without jet $p_{T}$ | 0.980 | 0.928 | 0.45
Trained with $\\{\tau_{x}^{(1)}\\}\bigcup\\{p_{T,J},m_{J}\\}$ | 0.976 | 0.921 | 0.43
Hidden Layers: $(200,200,50)$ | 0.980 | 0.928 | 0.48
PFN | Baseline | 0.997 | 0.977 | 0.708 $(\Phi)$, 0.467 $(F)$
$\Phi:(100,100,128)$ | 0.997 | 0.976 | 0.628 $(\Phi)$, 0.463 $(F)$
$\Phi:(100,100,64)$ | 0.994 | 0.971 | 0.633 $(\Phi)$, 0.480 $(F)$
$\Phi:(100,64,64),F:(64,32)$ | 0.994 | 0.969 | 0.535 $(\Phi)$, 0.531 $(F)$
| $\Phi:(100,64,64),F:(32,32)$ | 0.991 | 0.963 | 0.601 $(\Phi)$, 0.500 $(F)$

Table A1: Performance and sparsity of the baseline and model variants for
different model architectures

## References

  * (1) T. Miller, _Explanation in artificial intelligence: Insights from the social sciences_ , _Artificial Intelligence_ 267 (2019) 1.
  * (2) D. Gunning, M. Stefik, J. Choi, T. Miller, S. Stumpf and G.-Z. Yang, _XAI—explainable artificial intelligence_ , _Science Robotics_ 4 (2019) eaay7120.
  * (3) P. Linardatos, V. Papastefanopoulos and S. Kotsiantis, _Explainable AI: a review of machine learning interpretability methods_ , _Entropy_ 23 (2020) 18.
  * (4) G. Vilone and L. Longo, _Explainable artificial intelligence: a systematic review_ , _arXiv preprint arXiv:2006.00093_ (2020) .
  * (5) M. Sahakyan, Z. Aung and T. Rahwan, _Explainable artificial intelligence for tabular data: A survey_ , _IEEE Access_ 9 (2021) 135392\.
  * (6) H. Yuan, H. Yu, S. Gui and S. Ji, _Explainability in graph neural networks: A taxonomic survey_ , _arXiv preprint arXiv:2012.15445_ (2020) .
  * (7) D. Turvill, L. Barnby, B. Yuan and A. Zahir, _A survey of interpretability of machine learning in accelerator-based high energy physics_ , in _2020 IEEE/ACM International Conference on Big Data Computing, Applications and Technologies (BDCAT)_ , pp. 77–86, 2020, DOI.
  * (8) Y.S. Lai, D. Neill, M. Płoskoń and F. Ringer, _Explainable machine learning of the underlying physics of high-energy particle collisions_ , _Physics Letters B_ 829 (2022) 137055.
  * (9) F. Mokhtar, R. Kansal, D. Diaz, J. Duarte, J. Pata, M. Pierini et al., _Explaining machine-learned particle-flow reconstruction_ , _arXiv preprint arXiv:2111.12840_ (2021) .
  * (10) D.E. Kaplan, K. Rehermann, M.D. Schwartz and B. Tweedie, _Top tagging: A method for identifying boosted hadronically decaying top quarks_ , _Phys. Rev. Lett._ 101 (2008) 142001.
  * (11) L.G. Almeida, S.J. Lee, G. Perez, I. Sung and J. Virzi, _Top quark jets at the LHC_ , _Phys. Rev. D_ 79 (2009) 074012.
  * (12) L.G. Almeida, S.J. Lee, G. Perez, G. Sterman and I. Sung, _Template overlap method for massive jets_ , _Physical Review D_ 82 (2010) 054034.
  * (13) T. Plehn and M. Spannowsky, _Top tagging_ , _Journal of Physics G: Nuclear and Particle Physics_ 39 (2012) 083001.
  * (14) The ATLAS Collaboration, _Identification of high transverse momentum top quarks in pp collisions at $\sqrt{s}=8$ tev with the ATLAS detector_, _Journal of high energy physics_ 2016 (2016) 1.
  * (15) The CMS Collaboration, _A Cambridge-Aachen (C-A) based Jet Algorithm for boosted top-jet tagging_ , Tech. Rep. CMS-PAS-JME-09-001, CERN, Geneva (2009).
  * (16) The CMS Collaboration, _Boosted Top Jet Tagging at CMS_ , Tech. Rep. CMS-PAS-JME-13-007, CERN, Geneva (2014).
  * (17) The ATLAS Collaboration, _Performance of top-quark and $w$-boson tagging with atlas in run 2 of the lhc_, _Eur. Phys. J. C_ 79 (2019) 1.
  * (18) CMS collaboration, _Identification of heavy, energetic, hadronically decaying particles using machine-learning techniques_ , _Journal of Instrumentation_ (2020) .
  * (19) J. Pearkes, W. Fedorko, A. Lister and C. Gay, _Jet constituents for deep neural network based top quark tagging_ , _arXiv preprint arXiv:1704.02124_ (2017) .
  * (20) L. Moore, K. Nordström, S. Varma and M. Fairbairn, _Reports of my demise are greatly exaggerated: $N$-subjettiness taggers take on jet images_, _SciPost Phys._ 7 (2019) 036.
  * (21) K. Datta and A. Larkoski, _How much information is in a jet?_ , _J. High Energy Phys._ 2017 (2017) 1.
  * (22) G. Louppe, K. Cho, C. Becot and K. Cranmer, _QCD-aware recursive neural networks for jet physics_ , _Journal of High Energy Physics_ 2019 (2019) 1.
  * (23) A. Butter, G. Kasieczka, T. Plehn and M. Russell, _Deep-learned top tagging with a Lorentz layer_ , _SciPost Physics_ 5 (2018) 028\.
  * (24) P.T. Komiske, E.M. Metodiev and J. Thaler, _Energy flow networks: deep sets for particle jets_ , _Journal of High Energy Physics_ 2019 (2019) 1.
  * (25) H. Qu and L. Gouskos, _Jet tagging via particle clouds_ , _Physical Review D_ 101 (2020) .
  * (26) S. Macaluso and D. Shih, _Pulling out all the tops with computer vision and deep learning_ , _Journal of High Energy Physics_ 2018 (2018) 1.
  * (27) M. Erdmann, E. Geiser, Y. Rath and M. Rieger, _Lorentz boost networks: autonomous physics-inspired feature engineering_ , _Journal of Instrumentation_ 14 (2019) P06006.
  * (28) S. Egan, W. Fedorko, A. Lister, J. Pearkes and C. Gay, _Long short-term memory (LSTM) networks with jet constituents for boosted top tagging at the lhc_ , _arXiv preprint arXiv:1711.09059_ (2017) .
  * (29) E.A. Moreno, O. Cerri, J.M. Duarte, H.B. Newman, T.Q. Nguyen, A. Periwal et al., _JEDI-net: a jet identification algorithm based on interaction networks_ , _The European Physical Journal C_ 80 (2020) 1.
  * (30) S. Gong, Q. Meng, J. Zhang, H. Qu, C. Li, S. Qian et al., _An efficient lorentz equivariant graph neural network for jet tagging_ , _arXiv preprint arXiv:2201.08187_ (2022) .
  * (31) G. Kasieczka, T. Plehn, A. Butter, K. Cranmer, D. Debnath, B.M. Dillon et al., _The machine learning landscape of top taggers_ , _SciPost Physics Proceedings_ 7 (2019) .
  * (32) K. Hornik, M. Stinchcombe and H. White, _Multilayer feedforward networks are universal approximators_ , _Neural networks_ 2 (1989) 359\.
  * (33) A. Chakraborty, S.H. Lim and M.M. Nojiri, _Interpretable deep learning for two-prong jet classification with jet spectra_ , _Journal of High Energy Physics_ 2019 (2019) 1.
  * (34) G. Agarwal, L. Hay, I. Iashvili, B. Mannix, C. McLean, M. Morris et al., _Explainable AI for ML jet taggers using expert variables and layerwise relevance propagation_ , _Journal of High Energy Physics_ 2021 (2021) 1.
  * (35) P. Shanahan, K. Terao and D. Whiteson, _Snowmass 2021 computational frontier CompF03 topical group report: Machine learning_ , _arXiv preprint arXiv:2209.07559_ (2022) .
  * (36) D. Seuß, _Bridging the gap between explainable AI and uncertainty quantification to enhance trustability_ , _arXiv preprint arXiv:2105.11828_ (2021) .
  * (37) C. Grojean, A. Paul, Z. Qian and I. Strümke, _Lessons on interpretable machine learning from particle physics_ , _Nature Reviews Physics_ (2022) 1.
  * (38) J. Duarte, S. Han, P. Harris, S. Jindariani, E. Kreinar, B. Kreis et al., _Fast inference of deep neural networks in FPGAs for particle physics_ , _Journal of Instrumentation_ 13 (2018) P07027.
  * (39) Y. Iiyama, G. Cerminara, A. Gupta, J. Kieseler, V. Loncar, M. Pierini et al., _Distance-weighted graph neural networks on fpgas for real-time particle reconstruction in high energy physics_ , _Frontiers in big Data_ (2021) 44\.
  * (40) A. Heintz, V. Razavimaleki, J. Duarte, G. DeZoort, I. Ojalvo, S. Thais et al., _Accelerated charged particle tracking with graph neural networks on FPGAs_ , _arXiv preprint arXiv:2012.01563_ (2020) .
  * (41) A. Roy and M.S. Neubauer, “Interpretability of an interaction network for identifying $h\to b\bar{b}$ jets, https://agenda.infn.it/event/28874/contributions/169216/.”
  * (42) M.S. Neubauer and A. Roy, _Explainable AI for high energy physics_ , _arXiv preprint arXiv:2206.06632. Contribution to Snowmass 2021_ (2022) .
  * (43) E.A. Moreno, T.Q. Nguyen, J.-R. Vlimant, O. Cerri, H.B. Newman, A. Periwal et al., _Interaction networks for the identification of boosted $h\to b\bar{b}$ decays_, _Physical Review D_ 102 (2020) 012010.
  * (44) “Top tagging dataset, available at: https://desycloud.desy.de/index.php/s/llbX3zpLhazgPJ6.”
  * (45) T. Sjöstrand, S. Ask, J.R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Computer Physics Communications_ 191 (2015) 159.
  * (46) J. De Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaitre, A. Mertens et al., _DELPHES 3: a modular framework for fast simulation of a generic collider experiment_ , _Journal of High Energy Physics_ 2014 (2014) 1.
  * (47) M. Cacciari, G.P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _Journal of High Energy Physics_ 2008 (2008) 063.
  * (48) M. Cacciari, G.P. Salam and G. Soyez, _FastJet user manual_ , _The European Physical Journal C_ 72 (2012) 1.
  * (49) J. Thaler and K. Van Tilburg, _Identifying boosted objects with N-subjettiness_ , _Journal of High Energy Physics_ 2011 (2011) 1.
  * (50) S.D. Ellis and D.E. Soper, _Successive combination jet algorithm for hadron collisions_ , _Physical Review D_ 48 (1993) 3160.
  * (51) G.C. Blazeya, J.R. Dittmannb, S.D. Ellisc, V.D. Elvirab, K. Framed, S. Grinsteine et al., _Run II jet physics_ , _QCD and Weak Boson Physics in Run II_ (2000) 47.
  * (52) M. Zaheer, S. Kottur, S. Ravanbhakhsh, B. Póczos, R. Salakhutdinov and A.J. Smola, _Deep sets_ , in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , pp. 3394–3404, 2017.
  * (53) R. Wang and K. Tang, _Feature selection for maximizing the area under the ROC curve_ , in _2009 IEEE International Conference on Data Mining Workshops_ , pp. 400–405, IEEE, 2009.
  * (54) J. van der Waa, E. Nieuwburg, A. Cremers and M. Neerincx, _Evaluating XAI: A comparison of rule-based and example-based explanations_ , _Artificial Intelligence_ 291 (2021) 103404.
  * (55) S. Jesus, C. Belém, V. Balayan, J. Bento, P. Saleiro, P. Bizarro et al., _How can i choose an explainer? an application-grounded evaluation of post-hoc explanations_ , in _Proceedings of the 2021 ACM Conference on Fairness, Accountability, and Transparency_ , pp. 805–815, 2021.
  * (56) J. Tang, S. Alelyani and H. Liu, _Feature selection for classification: A review_ , _Data classification: Algorithms and applications_ (2014) 37.
  * (57) M.T. Ribeiro, S. Singh and C. Guestrin, _"why should I trust you?" explaining the predictions of any classifier_ , in _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_ , pp. 1135–1144, 2016.
  * (58) X.-w. Chen and M. Wasikowski, _Fast: a ROC-based feature selection metric for small samples and imbalanced data classification problems_ , in _Proceedings of the 14th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pp. 124–132, 2008.
  * (59) A.J. Serrano, E. Soria, J.D. Martin, R. Magdalena and J. Gomez, _Feature selection using ROC curves on classification problems_ , in _The 2010 international joint conference on neural networks (IJCNN)_ , pp. 1–6, IEEE, 2010\.
  * (60) S.M. Lundberg and S.-I. Lee, _A unified approach to interpreting model predictions_ , in _Advances in Neural Information Processing Systems_ , I. Guyon, U.V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan et al., eds., vol. 30, Curran Associates, Inc., 2017, https://proceedings.neurips.cc/paper/2017/file/8a20a8621978632d76c43dfd28b67767-Paper.pdf.
  * (61) M.T. Ribeiro, S. Singh and C. Guestrin, _Model-agnostic interpretability of machine learning_ , _arXiv preprint arXiv:1606.05386_ (2016) .
  * (62) A. Binder, S. Bach, G. Montavon, K.-R. Müller and W. Samek, _Layer-wise relevance propagation for deep neural network architectures_ , in _Information science and applications (ICISA) 2016_ , pp. 913–922, Springer (2016).
  * (63) G. Montavon, A. Binder, S. Lapuschkin, W. Samek and K.-R. Müller, _Layer-wise relevance propagation: an overview_ , _Explainable AI: interpreting, explaining and visualizing deep learning_ (2019) 193.
  * (64) S. Bach, A. Binder, G. Montavon, F. Klauschen, K.-R. Müller and W. Samek, _On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation_ , _PloS one_ 10 (2015) e0130140.
  * (65) T. Schnake, O. Eberle, J. Lederer, S. Nakajima, K.T. Schutt, K.-R. Mueller et al., _Higher-order explanations of graph neural networks via relevant walks_ , _IEEE Transactions on Pattern Analysis & Machine Intelligence_ (2021) 1.
  * (66) L. Toloşi and T. Lengauer, _Classification with correlated features: unreliability of feature ranking and solutions_ , _Bioinformatics_ 27 (2011) 1986.
  * (67) G. Montavon, S. Lapuschkin, A. Binder, W. Samek and K.-R. Müller, _Explaining nonlinear classification decisions with deep taylor decomposition_ , _Pattern recognition_ 65 (2017) 211.
  * (68) B.O. Ayinde, T. Inanc and J.M. Zurada, _Regularizing deep neural networks by enhancing diversity in feature extraction_ , _IEEE transactions on neural networks and learning systems_ 30 (2019) 2650.
  * (69) M. Cogswell, F. Ahmed, R. Girshick, L. Zitnick and D. Batra, _Reducing overfitting in deep networks by decorrelating representations_ , _arXiv preprint arXiv:1511.06068_ (2015) .
  * (70) H. Kaur, H. Nori, S. Jenkins, R. Caruana, H. Wallach and J. Wortman Vaughan, _Interpreting interpretability: understanding data scientists’ use of interpretability tools for machine learning_ , in _Proceedings of the 2020 CHI conference on human factors in computing systems_ , pp. 1–14, 2020.
  * (71) N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever and R. Salakhutdinov, _Dropout: a simple way to prevent neural networks from overfitting_ , _The journal of machine learning research_ 15 (2014) 1929.
  * (72) D.P. Kingma and M. Welling, _Auto-encoding variational bayes_ , _arXiv preprint arXiv:1312.6114_ (2013) .
  * (73) C.P. Burgess, I. Higgins, A. Pal, L. Matthey, N. Watters, G. Desjardins et al., _Understanding disentangling in $\beta$-VAE_, _arXiv preprint arXiv:1804.03599_ (2018) .
  * (74) G. Hadjeres, F. Nielsen and F. Pachet, _GLSR-VAE: Geodesic latent space regularization for variational autoencoder architectures_ , in _2017 IEEE Symposium Series on Computational Intelligence (SSCI)_ , pp. 1–7, IEEE, 2017\.
  * (75) C. Bajaj, A. Roy and H. Zhang, _Invariance-based multi-clustering of latent space embeddings for equivariant learning_ , _arXiv preprint arXiv:2107.11717_ (2021) .
  * (76) B. Bortolato, A. Smolkovič, B.M. Dillon and J.F. Kamenik, _Bump hunting in latent space_ , _Physical Review D_ 105 (2022) 115009\.
  * (77) Z. Liu, P. Luo, X. Wang and X. Tang, _Deep learning face attributes in the wild_ , in _Proceedings of International Conference on Computer Vision (ICCV)_ , December, 2015.
  * (78) B. Dillon, T. Plehn, C. Sauer and P. Sorrenson, _Better latent spaces for better autoencoders_ , _SciPost Physics_ 11 (2021) 061.
  * (79) I.T. Jolliffe and J. Cadima, _Principal component analysis: a review and recent developments_ , _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 374 (2016) 20150202\.


[a, b]Christopher Eckner

# Detecting and characterizing pulsar halos with the Cherenkov Telescope Array
Observatory

###### Abstract

The recently identified source class of pulsar halos may be populated and
bright enough at TeV energies to constitute a large fraction of the sources
that will be observed with the Cherenkov Telescope Array (CTA), especially in
the context of the planned Galactic Plane Survey (GPS). In this study, we
examine the prospects offered by CTA for the detection and characterization of
such objects. CTA will cover energies from 20 GeV to 300 TeV, bridging the
ranges already probed with the Fermi Large Area Telescope and High Altitude
Water Cherenkov Observatory, and will also have a better angular resolution
than the latter instruments, thus providing a complementary look at the
phenomenon. From simple models for individual pulsar halos and their
population in the Milky Way, we examine under which conditions such sources
can be detected and studied from the GPS observations. In the framework of a
full spatial-spectral likelihood analysis, using the most recent estimates for
the instrument response function and prototypes for the science tools, we
derive the spectral and morphological sensitivity of the CTA GPS to the
specific intensity distribution of pulsar halos. From these, we quantify the
physical parameters for which pulsar halos can be detected, identified, and
characterized, and what fraction of the Galactic population could be
accessible. We also discuss the effect of interstellar emission and data
analysis systematics on these prospects.

##  1 Introduction

Pulsars – high-energy astrophysical objects born in violent supernova
explosions – have claimed a prominent position in contemporary astroparticle
physics in the past decade. Not only are their emissions observable over a
broad range of the electromagnetic wave band, but they may also be one of the
main contributors to the observed local flux of cosmic-ray leptons [1]. The
Fermi-LAT has detected over 270 pulsars due to their gamma-ray emission at GeV
energies, far exceeding detections in the pre-Fermi era [2, 3]. Ground-based
Imaging Atmospheric Cherenkov Telescopes (IACTs), such as H.E.S.S., advanced
our understanding of the TeV emission of pulsars during the later stages of
their evolution, i.e. when they have formed so-called pulsar wind nebulae
(PWNe) [4, 5]. The gamma-ray emission patterns associated with PWNe are
influenced by the pulsar’s natal kick and the dynamics of the parent supernova
remnant. They become complex as pulsar, PWNe and supernova remnant evolve [6,
7]. Observations with the High-Altitude Water Cherenkov observatory (HAWC)
revealed extended gamma-ray emission, initially referred to as TeV halos [8]
but more recently replaced by the term _pulsar halos_ , around two middle-aged
pulsars – PSR J0633+1746 (Geminga) and PSR B0656+14 (Monogem) –, providing new
insights into particle acceleration and escape mechanisms [9]. The pulsar
halos’ emission appears to be fueled by the central pulsar since, in those two
cases, the central pulsar has left the region dynamically dominated by the
associated PWN. Interestingly, it is possible to describe the observed pulsar
halo emission with a phenomenological diffusion model featuring a suppressed
diffusion coefficient with respect to the values inferred in the remainder of
the Milky Way. The cause of the inhibited diffusion is not well understood and
subject of ongoing research [10, 11, 12]. Despite these observations, the
commonness of pulsar halos in the Galaxy is unknown. Upcoming observations
with HAWC, LHAASO, H.E.S.S., and the Cherenkov Telescope Array (CTA) are
expected to contribute significantly to this field [13, 14, 15]. This work –
whose unabridged content is published at [16] – aims to use the upcoming CTA
Galactic Plane Survey (GPS) [17] to constrain current phenomenological pulsar
halo models. We will determine the survey’s sensitivity to the $\gamma$-ray
spectrum and morphology of individual halos, and estimate the fraction of the
population detectable by CTA.

##  2 Sky model

This study focuses on the emission components considered in our simulations
and analyses of CTA observations. These components include the
phenomenological halo model, a synthetic halo population for the entire
Galaxy, and various models for the astrophysical background components. The
halo model has been utilized to predict the detectability of individual halos,
while the synthetic halo population is used to determine the potential
detectability fraction within the context of CTA’s GPS observations [18, 19].

Individual pulsar halo model. The individual halo model, introduced by [18],
simulates the behavior of electron-positron pairs accelerated and released
into the surrounding medium by a pulsar shortly after its birth. These
particles diffuse isotropically in a medium with a two-zone structure
representing different diffusion properties. Particles lose energy and radiate
via synchrotron and inverse-Compton scattering processes in assumed typical
magnetic and radiation fields of the Interstellar Medium (ISM) [19]. The
reference model is a pulsar with a current age of 200 kyr and a spin-down
power of $10^{35}$ erg/s, with particle injection assumed to start 60 kyr
after birth [20, 18]. The injected particles have a broken power-law spectrum
up to a cutoff at 1 PeV. Diffusion suppression is by a factor 500 within a
diffusion zone region of $r_{\mathrm{diff}}=50$ pc radius.

Pulsar halo population model. The model setups mentioned above are used to
ascertain the conditions for the potential observation of a typical halo in
the CTA GPS [21]. Using a synthetic population of young pulsars, we estimate
the number of objects potentially accessible to the survey. However, the
detectability criteria are only strictly valid for one set of halo model
parameters, therefore, the prospects at the population level should be taken
as approximations.

Large-scale diffuse backgrounds. In our analyses, we included a model for
large-scale interstellar emission (IE) from the Galactic population of cosmic
rays interacting with the ISM. While this component has been well mapped at
GeV energies with the Fermi-LAT, it is less established in the TeV range [22].
To model this, we refer to a recent study [23] based on available GeV to PeV
gamma-ray data and local, charged cosmic ray measurements. In terms of
consistency, we also consider a model adopted in a reference prospect study of
the CTA GPS [24].

##  3 Data Simulations and Analysis

###  3.1 Survey Simulation

The upcoming CTA GPS promises extensive coverage of the Galactic Plane at TeV
energies, comprising short and long-term programs totalling 1620 hours of
observation over a decade [15]. Our simulation follows the strategy outlined
by [17], ensuring optimal instrument performance by minimizing Moon light
contamination and zenith angle. The survey’s varying exposure is divided into
five distinct segments: Inner Galaxy, Cygnus/Perseus, Anticenter, and two
chunks for Vela/Carina [15].

The Inner Galaxy, the region with the youngest pulsars and thus pulsar halos,
receives the most extensive exposure. Our reference coordinates for this study
are $(l,b)=(-10^{\circ},0^{\circ})$ within this region, despite it being a
conservative choice given the intensity of IE.

We use the latest instrument response functions (IRFs) labeled prod5 for the
different configurations at the North and South sites. These files also
quantify the expected irreducible instrumental background (CR) of CTA due to
charged cosmic rays misclassified as gamma-ray events. The data files are
publicly available at [25]. We consider the initial "alpha" configuration to
reflect the real number of telescopes given the budget. We use the optimized
event reconstruction quality and background cuts derived from Monte Carlo
simulations of 50-hour observations. In Fig. 1 we display the exposure
achieved with the GPS in the inner Galaxy and overlay this map with the pulsar
halo population considered in this work.

Figure 1: A synthetic pulsar halo population overlaid on an exposure map for
the central regions surveyed in the GPS.

###  3.2 Analytical Framework

To analyze the potential for detecting and studying pulsar halos with CTA, we
apply the statistical inference framework of [26]. We employ a three-
dimensional likelihood analysis for binned data and Poisson statistics. Our
region of interest (ROI) is centered on $(l,b)=(-10^{\circ},0^{\circ})$ with a
size of $6^{\circ}\times 6^{\circ}$ and pixel size of $0.02^{\circ}\times
0.02^{\circ}$. The energy range is from 0.1 to 100 TeV, binned logarithmically
with 15 bins per decade.

Our signal model in the ROI is a linear combination of templates representing
the expected instrumental background $\\{B_{l}\\}_{l\in L}$ and astrophysical
signal components $\\{S_{k}\\}_{k\in K}$. The likelihood function is given by:

$\mathcal{L}\\!\left(\left.\bm{\mu}\right|\bm{n}\right)=\prod_{i,j}\frac{\mu_{ij}^{n_{ij}}}{\left(n_{ij}\right)!}e^{-\mu_{ij}}\mathrm{,}$
(1)

where $\bm{\mu}$ represents the model used to fit $\bm{n}$, the CTA
experimental data (mock data).

Our model is generally expressed as:

$\bm{\mu}=\sum_{k\in K}\bm{S}_{k}(\theta^{S}_{k})+\sum_{l\in
L}\bm{B}_{l}(\theta^{B}_{l})\mathrm{,}$ (2)

where $\\{\theta^{S}_{k}\\}_{k\in K}$ and $\\{\theta^{B}_{l}\\}_{l\in L}$
represent multiplicative parameters adjusting the spectral and angular
dependence of the signal and instrumental background templates. For
simplicity, we often consider a model with a global renormalization parameter
of predefined templates.

The capability to distinguish between alternative hypotheses regarding the
measured signal is assessed using the log-likelihood ratio test statistic
(TS). The TS is given by:

$\textrm{TS}=2\left(\ln\left[\frac{\mathcal{L}\left(\bm{\mu}({\hat{\theta}^{S_{\rm
test}}k},{\hat{\theta}^{B}_{l}})|\bm{n}\right)}{\mathcal{L}\left(\bm{\mu}({\hat{\theta}^{S_{\rm
null}}_{k}},{\hat{\theta}^{B}_{l}})|\bm{n}\right)}\right]\right)\mathrm{,}$
(3)

where the hatted quantities refer to the best-fit values for all model
parameters.

Estimation of Detectability. The detectability of astrophysical components can
be quantified using the test statistic (TS) given by Eq. 3 when assuming that
in the null hypothesis $\hat{\theta}_{k}^{S_{\rm null}}\equiv 0$. Under
relatively mild conditions satisfied here, the TS distribution follows a non-
central chi-square distribution with $K$ degrees of freedom. In most cases, we
examine a single signal template, which simplifies the problem to a half-chi-
square distribution with a single degree of freedom. Signal detection at the
5-sigma level corresponds to a TS value of approximately 25.

Spectral Sensitivity Analysis. We evaluated both spectral sensitivity and
model-independent angular sensitivity. As concerns the spectral sensitivity,
we use the entire energy range and derive the sensitivity energy bin by energy
bin. This sensitivity measures the flux level that results in detection with a
TS of 25 and at least ten detected gamma-ray events in each energy bin. The
process involves preparing mock datasets featuring the CR background as well
as IE and the pulsar halo source at a given normalization. We fit the full
model and background-only model (CR + IE) until we find the value of the
pulsar halo normalization that yields a TS of 25 over the null hypothesis.

Angular Decomposition of a Pulsar Halo Signal. We decomposed individual pulsar
halos using a step-by-step fitting process of growing concentric annuli
obtained from the original pulsar halo model. We evaluated Eq. 3 iteratively,
either expanding the annulus width or adding a new annulus to the model based
on the resulting TS value. This procedure continued until no significant
signal was observed. Finally, the parameters and errors for all significant
annuli were considered as the recovered angular decomposition of the input
signal.

Handling of Systematic Errors. In order to deal with systematic uncertainties,
we used an effective Poisson likelihood function [27, 26]. This function
incorporates nuisance parameters per pixel and energy bin of the binned model,
which mimic Gaussian noise affecting each pixel independently. However, this
method only accounts for uncertainties that linearly affect the number of
expected gamma rays, such as CTA’s effective area. Systematic errors were
integrated into our sensitivity derivations by adding these nuisance
parameters to the set of background parameters $\\{\theta^{B}_{l}\\}_{l\in
L}$.

##  4 CTA’s sensitivity to pulsar halos

###  4.1 Exploring the sensitivity to the phenomenological pulsar halo model

Spectral Sensitivity. The spectral sensitivity, visualized in the left panel
of Fig. 2, declines with an enlarged halo extension, a consequence of a closer
positioning. For a more intuitive understanding, the sensitivities to the
considered pulsar halo models align approximately with a Gaussian source of
the corresponding angular size. The sensitivity to the physical pulsar halo
model shows less energy variation due to the energy-dependent shrinking of the
source size. This effect is primarily due to the increase in energy losses
with increasing particle energy, offsetting the deteriorating instrument
sensitivity beyond roughly 10 TeV. At larger distances, exceeding
approximately 4-5 kpc, the sensitivity trend starts to flatten, as the halo’s
typical angular size becomes similar to the instrument’s angular resolution.
Meanwhile, the flux from a halo decreases inversely with the square of the
distance. The flux levels of our reference individual halo model are portrayed
in the left panel of Fig. 2. Spectral studies ranging from hundreds of GeV to
tens of TeV seem feasible for nearby halos at distances of 1-3 kpc, involving
injection powers a few times above that deduced for Geminga, or roughly
$10^{35}$ erg/s.

Impact of the IE Model: The right panel of Fig. 2 explores the effect of three
assumptions: the scenario without IE, the $\gamma$-optimized Min model [23],
and the IE model used in the GPS publication [17]. The IE impact is more
significant for an extended source, though limited to a $20\%$ impact.

Systematic Uncertainty: The analysis here relies heavily on accurate models
for emission components. However, many instrumental effects, particularly
those affecting smaller-scale field-of-view observations, are difficult to
model but could have a significant influence. Following [26], we consider an
impact scale of 0.1∘ and the impact magnitude of $1\%$ and $3\%$. The impact
of these effectively modeled systematic uncertainties is largest at the lowest
GeV energies considered in our study. This regime is already systematics-
dominated due to the large amount of available gamma-ray events. The effect of
these uncertainties is less than $10\%$ around 1 TeV, the energy range with
the highest sensitivity to this kind of gamma-ray source class.

Angular Sensitivity: In the left panel of Fig. 3 we present the results
regarding the angular decomposition of the baseline pulsar halo model with a
diffusion zone size of 50 pc at a distance of 1 kpc from Earth. The
experimental CTA mock data contains the halo source as well as CR background
and IE (Base Max). These results highlight the potential of CTA measurements
to study the angular profile of pulsar halos. In particular, the 0.1 - 1 TeV
energy range allows us to probe the halo’s morphology up to 30 pc from the
central pulsar and, in addition, it is possible to distinguish between
different parameters of the phenomenological pulsar halo model like the
diffusion zone size $r_{\mathrm{diff}}$.

###  4.2 Prospects for a Galactic pulsar halo population

We estimate the fraction of the pulsar halo population that should be
detectable with CTA, assuming Geminga-like diffusion properties and supposing
that all middle-aged pulsars develop a halo. The results of this population
study are shown in the right panel of Fig. 3. We assess different properties
of the halo population that CTA is capable of analyzing:

  1. 1.

Injection power such that the simulated halo signal is detected with a TS of
25 over the full energy range, using the true halo model in the fit process
(red line).

  2. 2.

Injection power such that a fit of the simulated halo signal with the true
halo model is significantly preferred over a simple energy-independent 2D
Gaussian intensity distribution (dashed purple line).

  3. 3.

Injection power such that a fit of the simulated halo signal with the true
halo model is significantly preferred over the true model truncated at 30 pc
from the pulsar (dashed orange line).

  4. 4.

Injection power such that a fit of the simulated halo signal with the true
halo model is significantly preferred over an alternative halo model having a
50% higher suppressed diffusion coefficient (dashed yellow line).

In this setup, CTA is expected to detect around 300 pulsar halos of which more
than 70 representatives will have a significant energy-dependent morphology.
In addition, up to 30 of these objects can be spatially decomposed to study
the diffusion conditions around the halo’s central pulsar. These prospects are
reduced by a factor of about four when assuming that all pulsar halos resemble
the case of Monogem. However, it must be noted that it is not clear at all
that the entirety of Galactic pulsars will develop halos. In such less
optimistic scenarios [18], the fraction of accessible pulsar halos may only be
5-10% of the stated values.

##  5 Conclusions

This study offers an evaluation of the potential for detecting and examining
the relatively novel astrophysical source class of pulsar halos through
planned CTA GPS observations. Using a basic phenomenological two-zone
diffusion model for individual pulsar halos and their population in the Milky
Way, we conduct a simulation to realistically study these entities through a
comprehensive spatial-spectral likelihood analysis of simulated survey
observations.

Assuming all middle-aged pulsars that exited their original PWN develop a
halo, and based on a halo model set up in line with HAWC observations of the
PSRJ0633+1746 halo, we estimate that around three hundred objects could
produce detectable emission components in the survey. However, only about a
third of these could be identified via their energy-dependent morphology, and
approximately a tenth would permit strong constraints on crucial physical
parameters, such as the magnitude or extent of suppressed diffusion around the
pulsar. These estimates are approximately reduced by four with a model setup
consistent with HAWC observations of PSRB0656+14.

Pulsar halos sustained by particle injection power in the range $10^{35-36}$
erg/s and situated a few kpc away should enable precise spectral studies from
a few hundreds of GeV to tens of TeV through CTA GPS observations. The $0.1-1$
TeV band accessible to CTA presents significant potential for constraining
halo transport properties. CTA is anticipated to complement HAWC and LHAASO by
extending the energy coverage below 1 TeV, where the emitting particles are
less influenced by energy losses and can thus explore a larger volume around
the pulsar.

Figure 2: (_Left_ :) Differential spectral energy sensitivity to a baseline
halo model, positioned at different distances from the observer. Results are
shown for analysis focusing on a $6^{\circ}$ region centred at
$(\ell,b)=(-10^{\circ},0^{\circ})$, and are overlaid with the halo model
intensities for halos at 1 and 13 kpc distance in orange and red dashed lines,
respectively. (_Right_ :) Differential sensitivity to our baseline halo model
at 1 kpc distance considering different diffuse emission models with respect
to our benchmark background model (cosmic rays + Base Max IEM). The dashed and
dash-dotted black lines show the benchmark model sensitivity with added
systematic uncertainties of 1% and 3%, respectively.

Figure 3: (_Left_ :) Model-dependent angular sensitivity to a pulsar halo at 1
kpc distance and $r_{\mathrm{diff}}=50$ pc studied under the conditions of the
GPS, overlaid with the baseline halo model predictions for a diffusion zone
size of 30/50/80 pc. The sensitivity analysis is decomposed into annuli with a
width of $0.2^{\circ}$, and split into three energy bands. (_Right_ :)
Sensitivity of the GPS in terms of particle injection power as a function of
distance. For each sensitivity curve, the numbers in parentheses indicate the
number of mock halos lying above the curve. Grey data points indicate the
positions and injection powers of known pulsars taken from the ATNF database
https://heasarc.gsfc.nasa.gov/W3Browse/all/atnfpulsar.html.

Acknowledgements. This research has made use of the CTA instrument response
functions provided by the CTA Consortium and Observatory, see http://www.cta-
observatory.org/science/cta-performance for more details. We gratefully
acknowledge financial support from the agencies and organizations listed here:
https://www.cta-observatory.org/consortium_acknowledgments/

## References

  * [1] C. Evoli, E. Amato, P. Blasi and R. Aloisio,  _Phys. Rev. D_ 103 (Apr., 2021) 083010, [2010.11955].
  * [2] D. J. Thompson,  _Reports on Progress in Physics_ 71 (Nov., 2008) 116901, [0811.0738].
  * [3] Fermi-LAT collaboration, July, 2022.
  * [4] H. Abdalla, A. Abramowski, F. Aharonian, F. Ait Benkhali, E. O. Angüner, M. Arakawa et al.,  _A &A_ 612 (Apr., 2018) A1, [1804.02432].
  * [5] H. Abdalla, A. Abramowski, F. Aharonian, F. Ait Benkhali, A. G. Akhperjanian, T. Andersson et al.,  _A &A_ 612 (Apr., 2018) A2, [1702.08280].
  * [6] H. Abdalla, F. Aharonian, F. Ait Benkhali, E. O. Angüner, M. Arakawa, C. Arcaro et al.,  _A &A_ 621 (Jan., 2019) A116, [1810.12676].
  * [7] G. Principe, A. M. W. Mitchell, S. Caroff, J. A. Hinton, R. D. Parsons and S. Funk,  _A &A_ 640 (Aug., 2020) A76, [2006.11177].
  * [8] T. Linden, K. Auchettl, J. Bramante, I. Cholis, K. Fang, D. Hooper et al.,  _Phys. Rev. D_ 96 (Nov., 2017) 103016, [1703.09704].
  * [9] A. U. Abeysekara, A. Albert, R. Alfaro, C. Alvarez, J. D. Álvarez, R. Arceo et al.,  _Science_ 358 (Nov., 2017) 911–914, [1711.06223].
  * [10] R. López-Coto and G. Giacinti,  _MNRAS_ 479 (Oct., 2018) 4526–4534, [1712.04373].
  * [11] K. Fang, X.-J. Bi and P.-F. Yin,  _MNRAS_ 488 (Sept., 2019) 4074–4080, [1903.06421].
  * [12] P. Mukhopadhyay and T. Linden,  _arXiv e-prints_ (Nov., 2021) arXiv:2111.01143, [2111.01143].
  * [13] R. López-Coto, E. d. O. Wilhelmi, F. Aharonian, E. Amato and J. Hinton,  _Nature Astronomy, 6:199, (2022)_ (2, 2022) , [2202.06899].
  * [14] H. Abdalla, F. Aharonian, F. Ait Benkhali, E. O. Angüner, C. Arcaro, C. Armand et al.,  _ApJ_ 917 (Aug., 2021) 6, [2107.01425].
  * [15] Cherenkov Telescope Array Consortium, B. S. Acharya, I. Agudo, I. Al Samarai, R. Alfaro, J. Alfaro et al., _Science with the Cherenkov Telescope Array_.  2019, 10.1142/10986#.
  * [16] C. Eckner, V. Vodeb, P. Martin, G. Zaharijas and F. Calore,  _Mon. Not. Roy. Astron. Soc._ 521 (2023) 3793–3809, [2212.11265].
  * [17] Q. Remy, L. Tibaldo, F. Acero, M. Fiori, J. Knödlseder, B. Olmi et al., in _37th International Cosmic Ray Conference. 12-23 July 2021\. Berlin_ , p. 886, Mar., 2022.  2109.03729.
  * [18] P. Martin, A. Marcowith and L. Tibaldo,  _Astron. Astrophys._ 665 (2022) A132, [2206.11803].
  * [19] X. Tang and T. Piran,  _MNRAS_ 484 (Apr., 2019) 3491–3501, [1808.02445].
  * [20] Y. Zhang, R.-Y. Liu, S. Z. Chen and X.-Y. Wang,  _arXiv e-prints_ (Oct., 2020) arXiv:2010.15731, [2010.15731].
  * [21] P. Martin, L. Tibaldo, A. Marcowith and S. Abdollahi,  _Astron. Astrophys._ 666 (2022) A7, [2207.11178].
  * [22] Fermi-LAT collaboration, M. Ackermann et al.,  _Astrophys. J._ 750 (2012) 3, [1202.4039].
  * [23] P. D. l. T. Luque, D. Gaggero, D. Grasso, O. Fornieri, K. Egberts, C. Steppa et al.,  _Astron. Astrophys._ 672 (2023) A58, [2203.15759].
  * [24] A. Dundovic, C. Evoli, D. Gaggero and D. Grasso,  _Astron. Astrophys._ 653 (2021) A18, [2105.13165].
  * [25] Cherenkov Telescope Array Observatory and Cherenkov Telescope Array Consortium, Sept., 2021.  10.5281/zenodo.5499840.
  * [26] CTA collaboration, A. Acharyya et al.,  _JCAP_ 01 (2021) 057, [2007.16129].
  * [27] H. Silverwood, C. Weniger, P. Scott and G. Bertone,  _JCAP_ 03 (2015) 055, [1408.4131].