# Classification of Fermi-LAT sources with deep learning using energy and time
spectra

Thorben Finke,1 Michael Krämer,1 Silvia Manconi,1
1Institute for Theoretical Particle Physics and Cosmology, RWTH Aachen
University, D-52056 Aachen, Germany
Preprint number: TTK-20-45 Corresponding author, E-mail<EMAIL_ADDRESS>aachen.de

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Machine learning techniques are powerful tools for the classification of
unidentified gamma-ray sources. We present a new approach based on dense and
recurrent deep neural networks to classify unidentified or unassociated gamma-
ray sources in the last release of the Fermi-LAT catalog (4FGL-DR2). Our
method uses the actual measurements of the photon energy spectrum and time
series as input for the classification, instead of specific, hand-crafted
features. We focus on the separation between extragalactic sources, i.e.
Active Galactic Nuclei, and Galactic pulsars, and on the further
classification of pulsars into young and millisecond pulsars. Our list of
high-confidence candidate sources labelled by the neural networks provides
targets for further multiwavelength observations to identify their nature. The
deep neural network architectures we develop can be easily extended to include
specific features as well as multiwavelength data on the source photon energy
and time spectra coming from different instruments.

###### keywords:

gamma-rays: general – galaxies: statistics – pulsar:general – catalogues –
software:data analysis

††pubyear: 2020††pagerange: Classification of Fermi-LAT sources with deep
learning using energy and time spectra–Classification of Fermi-LAT sources
with deep learning using energy and time spectra

##  1 Introduction

The Fermi Gamma-ray Space Telescope, now operating for more than ten years,
has revolutionised our knowledge of non-thermal gamma-rays produced in our
Galaxy and beyond. Its primary instrument, the Fermi Large Area Telescope
(LAT), is a pair-conversion telescope which surveys the entire sky, detecting
gamma-rays with energies from few tens of MeV up to the TeV range (Atwood et
al., 2009; Ackermann et al., 2012a). The data obtained by Fermi-LAT has led to
the discovery of more than five thousand Galactic and extragalactic pointlike
(and few extended) gamma-ray sources, enlarging the known gamma-ray source
populations to new objects such as radio-quiet pulsars, millisecond pulsars,
and globular clusters (Abdollahi et al., 2020; Ballet et al., 2020). All-sky
gamma-ray data permitted also to characterise Galactic diffuse emissions
(Ackermann et al., 2012b; Acero et al., 2016) and extragalactic diffuse
backgrounds with unprecedented accuracy (Ackermann et al., 2015; Fornasa &
Sánchez-Conde, 2015).

Observed sources are characterised according to their position and spectral
characteristics, and have been collected in a succession of public source
catalogs. The fourth Fermi-LAT gamma-ray catalog, Data Release 2 (4FGL-DR2)
(Ballet et al., 2020) is the latest release, based on the first ten years of
data taking with gamma-ray energies ranging from 50 MeV to 1 TeV. The 4FGL-DR2
catalog is obtained with the same analysis methods as the 4FGL (Abdollahi et
al., 2020), but considering two additional years of data and thus adding 723
new sources, resulting in a total of 5788 detected sources.

The most numerous and brightest source population is currently represented by
Active Galactic Nuclei (AGN), where jets originating from a supermassive black
hole at the center of a galaxy emit electromagnetic radiation in a broad
energy range, from radio frequencies to gamma-rays with TeV energies (Urry &
Padovani, 1995; Padovani et al., 2017). More information on the various
classes of AGN detected by Fermi-LAT is contained in dedicated catalogs (The
Fermi-LAT Collaboration, 2019).

Aside from AGN, the second main class of sources contained in present catalogs
consists of gamma-ray pulsars (PSR), divided further into young pulsars (YNG)
and millisecond pulsars (MSP) (Abdo et al., 2013). The gamma-ray emission of a
pulsar can be attributed to a rapidly rotating, highly magnetised neutron
star, surrounded by a plasma magnetosphere, although the exact emission
mechanism and the location of the emission region are still uncertain (Harding
& Muslimov, 1998; Romani, 2014). MSP, or so-called recycled pulsars, are
believed to be old pulsars for which the period has decreased via accretion
from a companion object (Abdo et al., 2009).

The AGN and PSR classes of gamma-ray emitters are characterised by different
timing and spectral properties. While AGN typically have a variable gamma-ray
flaring pattern on time scales of months, PSR are instead non-variable on long
time scales (Abdo et al., 2013). Moreover, the energy spectra of AGN are
generically well described by a power law shape, which is softer than that of
a PSR. The energy spectra of PSR are characterised by a spectral curvature,
where a break is observed at low and high energies in the Fermi-LAT energy
band, requiring the addition of an exponential cutoff at a few GeV.

In the 4FGL-DR2 catalog the gamma-ray sources are identified as, or associated
with, an AGN, a PSR, or some more rare population class. Sources are
identified according to a correlated timing signature at different
wavelengths, or only associated requiring positional coincidence with a
counterpart source at other wavelengths. A total of 1794 sources (about one
third) in the 4FGL-DR2 catalog remains unassociated or unidentified (for
brevity, we will label both unassociated and unidentified sources as UNC,
unclassified, in what follows). The fraction of UNC sources was similar in
previous catalogs (Acero et al., 2015). The distribution of their spectral
properties and their Galactic latitude reveal the composite nature of the UNC
population (Abdollahi et al., 2020).

Each source in a gamma-ray catalog is characterised by measurements as well as
human-crafted features (a total of 74 in the 4FGL-DR2) designed to represent
the main source properties. Apart from the features connected to the source
position in the sky, many features describe the energy and timing properties
of photon fluxes. The 4FGL-DR2 catalog reports the photon fluxes in seven
energy bands, as well as the integrated flux in the full energy range of the
analysis. The photon flux is then fitted with three different spectral shapes
over the full energy range, and the best fit parameters, together with their
uncertainties, are provided. As for the time variability, in 4FGL-DR2 one-year
light curves (integrated over the full energy range) are reported, along with
the variability index, which weights the relative variability between
different time bins with respect to the mean value of all bins.

Given the rapid growth in size and complexity of astronomical datasets,
machine learning algorithms have become increasingly powerful tools, in
particular for object classification tasks, see e.g. Ball & Brunner (2010);
Baron (2019). The number of gamma-ray sources detected by Fermi-LAT has
increased from 1451 in the first catalog obtained with 11 months of data
taking (Abdo et al., 2010a) to 5788 in the most recent 4FGL-DR2 catalog based
on ten years of data. Because of the large amount of data, and the subjective
nature of classification based on human-crafted features, various authors have
developed automated schemes based on machine learning techniques in order to
classify candidate gamma-ray sources. These approaches are generally based on
information collected in gamma-ray catalogs only, and provide a list of
particularly interesting candidate target sources for follow-up observations.
Machine learning has been applied to characterise the UNC sources in current
and past Fermi-LAT catalogs (Mirabal et al., 2012, 2016; Saz Parkinson et al.,
2016; Hui et al., 2020; Luo et al., 2020), or to predict the subclass of AGN
of uncertain type (Doert & Errando, 2014; Chiaro et al., 2016; Salvetti et
al., 2017; Kovačević et al., 2020). Furthermore, machine learning based
analyses are instrumental to population studies, pointing to new source
classes (Saz Parkinson et al., 2016) or to a more exotic origin of gamma-ray
emission, such as dark-matter sub-halos (Mirabal et al., 2012).

The studies thus far typically use standard classification algorithms, such as
classification trees, logistic regression, or random forest classifiers, and
are based on the features provided by the Fermi-LAT catalogs, or simple
additional features, such as hardness ratios, constructed from the energy- and
time dependent gamma-ray fluxes. Such an approach requires a careful choice of
the features to consider, introducing ambiguities and potential biases e.g.
due to existing knowledge of the source sample. To automatise and improve the
source classification performance, machine learning techniques for feature
selection were recently applied to the classification of 4FGL sources (Luo et
al., 2020). The analysis by Luo et al. (2020) demonstrates that feature
selection by machine learning can improve the classification performance
compared to analyses where the choice of features is guided by prior
knowledge. The result of this approach depends on the particular release of
the catalog used, as the number and choice of features may change (e.g. the
introduction of a new spectral form). Moreover, the features might not
represent all the information encoded in the spectral data, and the physical
interpretation of the importance of some features for the classifier, for
example errors in the flux measurement, may be difficult.

In this paper we propose a different approach to classify UNC gamma-ray
sources. We employ a fully connected deep neural network and a recurrent
neural network trained directly on the energy and time dependent fluxes
presented in the Fermi-LAT catalog, instead of using human-crafted features as
input. While a similar approach has been used to predict the subclass of AGN
of uncertain type in Ref.(Kovačević et al., 2020), deep learning based
directly on energy and time dependent fluxes has to our knowledge thus far not
been used for the classification of UNC sources. We demonstrate that deep
neural networks based on energy and time spectra provide a powerful classifier
to separate AGN from PSR, as well as YNG from MSP, with a performance that is
comparable with previous feature-based analyses presented in the literature.

In addition, our method is very flexible and can easily be generalised to
include multiwavelength data for the energy and time spectra coming from
different observatories.

We also reconsider the classification of gamma-ray sources based on features,
and propose a feature selection procedure based on autoencoders. The feature
selection with an autoencoder does not rely on the source labels and thus may
help to reduce the problem of a sample selection bias, which may occur when
the distributions of the training and the test data sets do not match.

The paper is organised as follows. The data processing, the main neural
network architectures and the training and testing strategies are presented in
Sec. 2. In Sec. 3 we discuss the classification performance on labelled
catalog sources. A cross-matching procedure using an older catalog is
presented in Sec. 4. In Sec. 5 we introduce an alternative approach for
feature selection using autoencoders. Our predictions for the unlabelled
sources in the current Fermi-LAT catalog are provided in Sec. 6. We conclude
in Sec. 7.

##  2 Classification with neural networks

The central idea of our classification procedure is to directly use the gamma-
ray energy spectra and time series, rather than a selection of derived
features. The reference dataset is the 4FGL-DR2
catalog111https://fermi.gsfc.nasa.gov/ssc/data/access/lat/10yr_catalog/
(gll_psc_v23.fit), which has been obtained with the same analysis methods as
the 4FGL Fermi-LAT catalog (Abdollahi et al., 2020), but considering two
additional years of data and thus updating spectral parameters, spectral
energy distributions and source associations.

We use the gamma-ray flux values in the seven energy bands as reported in the
4FGL-DR2 catalog (entry Flux_Band in the catalog tables), together with the
corresponding significance in each energy band, referred simply as
’significance’ in what follows (Sqrt_TS_Band).222The source significance is
quantified by the test statistics TS = $2\log(\mathcal{L}/\mathcal{L}_{0})$,
where $\mathcal{L}$ and $\mathcal{L}_{0}$ are the log-likelihood functions
with and without the source of interest included in the model, respectively.
These flux values are obtained from the spectral fit over the full energy
range, but adjusting the normalisation in each spectral band.

As for the time spectra, we select the ten one-year light curves and the
corresponding significance (integrated over the full energy range)
(Flux_History and Sqrt_TS_History). Note that two-month intervals were
provided in the first release of the 4FGL catalog, but removed in the 4FGL-DR2
release, since one-year light curves were demonstrated to capture most of the
variability information (Abdollahi et al., 2020).

A cross-match of the classification performance based on the 3FGL catalog
(Acero et al., 2015) will be presented in Sec. 4. Some differences in the
structure of the energy and time spectra are present in this previous catalog:
the energy bands are reduced to five, and the time spectra are represented by
48 monthly bins, instead of yearly bins.

As will be evident from the subsequent discussion, our neural network
architectures are very flexible and additional features, apart from the energy
and time spectra, can easily be included in the classification algorithms. We
will explore the impact of additional features on the classification
performance in Sec. 3, in particular for the task of differentiating between
the YNG and MSP classes of pulsars.

As a first step in data preparation, we select sources from the 4FGL-DR2
catalog falling in the two broad classes of AGN and PSR. The PSR class
includes all sources that are identified or associated as pulsars
(CLASS1=PSR,psr). As for the AGN class, we include all sources that are
identified or associated with any type of AGN, including blazars of unknown
type, i.e. with CLASS1=FSRQ, fsrq, BLL, bll, BCU, bcu, CSS, css, RDG, rdg,
NLSY1, nlsy1, agn, ssrq, sey. We will refer to identified/associated AGN and
PSR as labelled catalog sources. The UNC sources are instead those without any
identification or association, reported as a blank space in the 4FGL-DR2
catalog. The number of the corresponding sources in the 3FGL and the 4FGL-DR2
catalogs is reported in Table 1. We note an increase of almost a factor of two
in the number of sources between the 3FGL and the 4FGL-DR2 catalog. In Fig. 1
we display the number of sources in the 4FGL-DR2 catalog for each class as a
function of the significance (in $\sigma$). The number of UNC sources is
subdominant at high significances, but increases monotonically towards lower
significances, and is similar to the number of AGN detected at $5\sigma$.

Table 1: Composition of the 4FGL-DR2 and 3FGL catalogs. Source class | 4FGL-DR2 | 3FGL
---|---|---
AGN | 3508 | 1745
PSR | 259 | 167
UNC | 1679 | 1010
Tot | 5788 | 3034
Figure 1: Number of sources in the Fermi-LAT 4FGL-DR2 catalog as a function of
the source significance (in $\sigma$ units) in the 100 MeV to 1 TeV energy
band, for the three main source classes: AGN, PSR and UNC (see text for more
details).

Since we are also interested in the classification of pulsars as YNG or MSP,
we use the public list of LAT-detected gamma-ray
pulsars333https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-
Detected+Gamma-Ray+Pulsars to build a cross-match between this list and the
4FGL-DR2 catalog. We obtain a list of 230 known gamma-ray pulsars, split
between 123 YNG and 107 MSP. Note that energy and time spectra are given for
each source in the catalog, and thus there is no need to filter out sources
with missing information, as necessary when using features. The numbers in
Table 1 thus represent the full statistics of the source sample used for the
classification.

###  2.1 Deep learning architectures

For the classification of the gamma-ray sources we use two standard deep
neural network architectures: a fully connected network, also called dense
neural network (DNN), and a recurrent neural network (RNN) (for an
introduction on deep neural networks and other supervised learning techniques
in physical sciences see e.g. Carleo et al. (2019)). RNNs are designed for
analysing sequential data or time-series data, and can process variable length
sequences of inputs. They have been used successfully for classification tasks
for example in astronomy (Ishida, 2019; Möller & de Boissière, 2020; Naul et
al., 2018; Hinners et al., 2018; Becker et al., 2020) and in high-energy
physics (Guest et al., 2016; Egan et al., 2017; Louppe et al., 2019; Cheng,
2018; Fraser & Schwartz, 2018; Englert et al., 2020).

Both the DNN and RNN consist of two branches, one for the energy spectrum and
one for the time series, which are eventually merged into a fully connected
network for classification. When adding additional features, we append a third
branch to the architecture, which is merged for classification with the other
two branches in the same way as before. We have used a ten-fold cross-
validation procedure to optimise the number of hidden layers and nodes for
each architecture. For this procedure the training set is divided into ten
subsets. The networks with different hyperparameters are trained on nine of
these subsets and evaluated on the tenth. This is repeated ten times, always
taking another subset for the evaluation. The cross-validation procedure
allows us to evaluate on the full training set, while keeping a test set that
has not been used for hyperparameter optimisation. In the case of the DNN, a
good performance can be achieved with one hidden layer for the time series and
one hidden layer for the energy spectrum, each consisting of 16 nodes. For the
RNN, we use three recurrent layers per branch with eight units each. For both
the DNN and the RNN, the classification is done with a fully connected network
with two layers, with four and two nodes for the first and second layer,
respectively.

In the fully connected layers we use $\text{ReLU}(x)=\max(0,x)$ as activation
function, except for the final layer where the softmax activation function
converts our output into a probability distribution. The RNN is based on the
SimpleRNN layers from Tensorflow/Keras (Abadi et al., 2015; Chollet et al.,
2015) which use a hyperbolic tangent activation function. A visualisation of
the DNN architecture is shown in Fig. 2 for illustration.

Figure 2: Sketch of the dense neural network architecture for the
classification of Fermi-LAT sources.

###  2.2 Data preprocessing, training, and testing

The input to the neural networks consists of the energy spectra and the time
series, together with the corresponding significances. For the RNN this
amounts to two series with two entries each, i.e. the flux and the
significance for each bin. In the case of the DNN we reshape the input into
two one-dimensional vectors consisting of $2\times 7$ entries for the energy
spectrum and $2\times 10$ entries for the time series. We shift and re-scale
each entry such that the mean over the training set is zero and the standard
deviation is one.

We use the Adam optimiser (Kingma & Ba, 2014) with an adapted initial learning
rate for training and the categorical cross entropy loss function. We find
that the most efficient training is achieved for the fully connected network
with a batch size of 64 and an initial learning rate of 0.01. For the RNN a
batch size of 32 and an initial learning rate of 0.001 performs best. We train
the network for 50 and 150 epochs for the task of AGN vs. PSR and MSP vs. YNG
classification, respectively, see also Sec. 3.2. The testing consists of two
parts. First, we perform a ten-fold cross-validation to find the best
architecture and training hyperparameters, as described above. After the
architectures and hyperparameters are fixed, we train our networks on 70% of
the labelled sources and evaluate on the remaining 30%. The test split is
performed such that the relative amount of the two classes is the same in the
training and test sets. The performance measures and uncertainties stated in
this paper correspond to the mean values and the standard deviations of ten
training runs. For each run a new network is initialised and the training and
test split is performed randomly. The uncertainties thus combine uncertainties
from different network initialisations and from the selected training data.

##  3 Performance on labelled sources in the 4FGL-DR2 catalog

Due to the large imbalance of the two classes in our dataset (see the AGN and
PSR numbers in Table 1), the measures used to quantify the performance of the
classifier need to be selected with some care. Indeed, the global accuracy of
a classifier (defined as the ratio of correctly classified sources over the
total number of sources) is sensitive to this imbalance. For example, for our
specific dataset the accuracy would approach already $\approx 93$% if the
classifier simply labelled all the sources as AGNs.

In what follows, we call PSR the positive class and AGN the negative one. With
this definition, correctly classified PSR are true positives (TP), while false
positives (FP) are AGN incorrectly classified as PSR. Correctly classified AGN
are true negatives (TN), while PSR incorrectly classified as AGN are false
negatives (FN). We consider the following four performance measures:

  1. 1.

Accuracy = $\frac{\displaystyle{\rm TP}+{\rm TN}}{\displaystyle{\rm TP}+{\rm
FP}+{\rm TN}+{\rm FN}}$

  2. 2.

F1 score = $\frac{\displaystyle{\rm TP}}{\displaystyle{\rm TP}+0.5({\rm
FP}+{\rm FN})}$

  3. 3.

Equal accuracy

  4. 4.

Area under the Receiver Operating Characteristic (ROC) curve.

We choose PSR as the positive class since we want to further differentiate
between the pulsar subclasses MSP and YNG. For the MSP vs. YNG classification
we arbitrarily refer to MSP as the positive class.

The accuracy is calculated at the decision boundary that gives the highest
accuracy on the training set. The decision boundary refers to the threshold in
the positive class score necessary for a source to be labelled as belonging to
the positive class. The F1 score is the harmonic mean of precision (the ratio
of correct PSR classifications to all PSR classifications, i.e. TP/(TP+FP))
and recall (the ratio of correct PSR classifications to all actual PSR, i.e.
TP/(TP+FN)) at the same decision boundary. The “equal accuracy" measure is
given by the accuracy of the classifier when operating at the decision
boundary that gives the most similar accuracy in both classes on the training
set. By scanning over all possible decision boundaries of the classifier, one
can construct the ROC curve. Here, the ROC curve corresponds to the PSR
efficiency (true positive rate) as a function of the AGN efficiency (false
positive rate). We recall that the area under the ROC curve (AUC) is bound
between zero and one. A perfect classifier would have an AUC value of one,
while a random classifier would correspond to an AUC value of 0.5. AUC values
below 0.5 indicate that the classifier is switching the classes. An AUC value
of zero corresponds to a classifier that inverts all labels.

###  3.1 AGN vs. PSR classification

We first discuss the AGN vs. PSR classification. Table 2 summarises the
results on the labelled data of the 4FGL-DR2 catalog. The results for the four
performance measures are reported for both the DNN and RNN benchmark
architectures as described in Sec.2.1. We first of all note that the DNN and
RNN perform similar on the task of classifying AGN vs. PSR. Our benchmark
architectures use the energy and time spectra, and the corresponding
significances, for each source. We also test some variation of these
benchmarks, where one of the two input spectra is substituted by standard
features in the catalog, such as the Variability_Index which quantifies the
variability of the gamma-ray flux with one single estimator (see the
definition in Abdollahi et al. (2020)).

Table 2: Performance of AGN vs. PSR classification using the 4FGL-DR2 catalog.  Architecture | Accuracy [%] | F1 score [%] | Equal acc [%] | AUC [%]
---|---|---|---|---
DNN | $97.59\pm 0.35$ | $81.73\pm 3.43$ | $95.67\pm 0.83$ | $97.99\pm 0.75$
RNN | $97.31\pm 0.36$ | $79.27\pm 2.82$ | $93.86\pm 1.53$ | $97.40\pm 0.89$

Using the Flux_Band in combination with the Variability_Index feature instead
of the Flux_History data, we obtain for the DNN an accuracy of $97.90\pm
0.44$% and a F1 score of $84.09\pm 3.27$%. Using only the Flux_Band data and
no additional information about the time dependence, we see also a similar
performance (accuracy of $97.70\pm 0.32$% and a F1 score of $82.30\pm 2.48$%).
When using only the Flux_History, we see a decrease in performance to an
accuracy of $95.46\pm 0.33$% and a F1 score of $64.04\pm 3.92$%. We conclude
that the essential information for the AGN vs. PSR classification is contained
in the Flux_Band, i.e. the energy spectra integrated over time. The
Flux_History data does not provide significant additional information. The
same is observed for the RNN architecture.

We show the ROC curves for the DNN and RNN in the top panel of Fig. 3. The DNN
performs slightly better for all PSR efficiencies. Overall, the curves are
close to the step function corresponding to an ideal classifier, and the
fluctuations, denoted by the coloured bands, are rather small.

Figure 3: ROC curves for AGN vs. PSR classification (top) and MSP vs. YNG
classification (bottom). The shaded area is given by the standard deviation
for the ten runs at each position.

Removing the significances for both series gives the same results for the mean
values of our performance metrics, while their deviations slightly decrease:
accuracy $97.60\pm 0.20$ ($97.67\pm 0.24$), F1 score $81.68\pm 1.62$
($81.82\pm 2.36$) for the DNN (RNN). The slight decrease in uncertainty
originates from the fewer inputs when not using the significances. Since the
significances do not contribute to the overall class separation, they can be
seen as additional noise leading to some overfitting and thus a stronger
variation on the test sets. We nevertheless keep the significance as input in
our standard set-up, as it provides valuable information for the further
classification of the PSR class into MSP and YNG, see Sec. 3.2.

Luo et al. (2020) state an accuracy of $99.19\pm 0.16$% for the AGN vs. PSR
classification using automated feature selection and a random forest
classifier on the 4FGL sources. While Luo et al. (2020) do not provide any
other performance measures, we achieve comparable results on the accuracy,
without the need to extract and select features from the catalog.

One way to address the imbalance of the data set is by weighting the loss of
the AGN and PSR classes such that the corresponding contributions to the total
loss are equal. We have also explored techniques for augmenting data of the
subdominant PSR class. First, we simply reused the actual PSR data as often as
necessary to balance the two classes. Second, we reused the PSR data but added
some Gaussian noise to the augmented data points. Third, we adopted the SMOTE
algorithm (Chawla et al., 2002), which can be formulated as $S_{\rm
aug}=S_{1}+r(S_{2}-S_{1})$: we pick two PSR sources randomly and add to the
first source $S_{1}$ the difference to the second source $S_{2}$, multiplied
by a random number $r$, drawn from a uniform distribution in $(0,1)$. The
difference is taken for each input that we use individually. We thus obtain an
augmented source $S_{\rm aug}$ with spectra interpolating between those of the
two original sources. Overall, we find that these different data augmentation
techniques do not improve the classification. While the performance at a
generic threshold of 0.5 changes, i.e. the accuracy in the PSR and AGN classes
increases and decreases, respectively, these effects can be compensated for by
a change in the decision threshold. Since we adapt this decision threshold
according to the performance in the training set, we achieve no benefit from
using the above-mentioned data augmentation techniques.

###  3.2 MSP vs. YNG classification

The classification of MSP vs. YNG pulsars is significantly more difficult,
given the very small data set (only 230 sources, with 127 YNG and 103 MSP, of
which we use 161 for training). Keeping the batchsize as in Sec. 3.1, the
network parameters are updated only few times during one epoch. We thus train
the networks more epochs for the task of MSP vs. YNG classification (150
instead of 50). Additionally to the Flux_Band and Flux_History information,
the latitude of the pulsar in Galactic coordinates (GLAT) is a distinguishing
feature that can be directly measured, and it has been used in previous
attempts to distinguish these two source classes (Luo et al., 2020; Saz
Parkinson et al., 2016). Given our flexible neural network architectures, we
can easily include such features in the classification.

In table 3 we report our results for the MSP vs. YNG classification task. We
find that adding the latitude increases the different performance measures by
about 3% for the DNN and 4% for the RNN. We also observe that removing the
Flux_History has no strong effect on the classification power. In contrast to
the AGN vs. PSR classification, we see a decrease in performance when leaving
out the significance. Both, DNN and RNN, drop below 80% accuracy and 0.9 AUC.
To be compared with our results, the feature selection algorithm from Luo et
al. (2020) achieves $89.61\pm 2.34$% for the MSP vs. YNG classification task
using a boosted logistic regression as classifier.

Table 3: Performance of MSP vs. YNG classification using the 4FGL-DR2 catalog. Architecture | Inputs | Accuracy [%] | F1 score [%] | Equal acc [%] | AUC [%]
---|---|---|---|---|---
DNN | History, Band | $83.91\pm 3.69$ | $84.47\pm 3.50$ | $83.91\pm 3.97$ | $90.25\pm 2.65$
Band, GLAT | $87.10\pm 3.86$ | $87.82\pm 3.63$ | $86.67\pm 4.84$ | $92.84\pm 3.41$
History, Band, GLAT | $87.25\pm 3.82$ | $88.00\pm 3.66$ | $86.67\pm 4.14$ | $93.34\pm 2.09$
RNN | History, Band | $82.03\pm 5.15$ | $82.69\pm 4.50$ | $83.33\pm 5.11$ | $90.03\pm 4.11$
Band, GLAT | $85.51\pm 3.37$ | $85.72\pm 3.90$ | $85.51\pm 4.00$ | $92.61\pm 2.83$
History, Band, GLAT | $86.52\pm 2.83$ | $87.46\pm 2.75$ | $84.78\pm 2.36$ | $93.04\pm 2.32$

In the bottom panel of Fig. 3 we show the ROC curves for the classification of
MSP vs. YNG using our standard architectures (without Galactic latitude). The
performances of the DNN and RNN are very similar over the full range of false
positive rates. The comparison to the AGN vs. PSR classification (top panel of
Fig. 3) shows that the sub-classification task is harder. The corresponding
ROC curves are lower and the uncertainties are larger. One reason for this
behaviour is the much smaller amount of labelled sources available for
training as well as for evaluation.

##  4 Cross-match between the 3FGL and 4FGL-DR2 catalogs

Given the increased gamma-ray statistics, and novel associations based on
updated catalogs at other wavelengths, some of the sources labelled as UNC in
the previous 3FGL catalog have now a defined association in the 4FGL-DR2
catalog. In total, 294 unclassified sources in the 3FGL catalog are now
labelled as AGN (258) or PSR (36) in the 4FGL-DR2 catalog. These sources
constitute a cross-match test set to evaluate the capability of our
classifiers, by training on the 3FGL catalog and testing the classification
for the previously UNC objects against the classification in the 4FGL-DR2
catalog. In addition, the cross-match test data set allows for an estimate of
a potential sample selection bias (also called covariate shift or Malmquist
bias in astronomy), see Pasquet et al. (2019); Luo et al. (2020). Data sets in
astrophysics are typically based on brightness-limited surveys. Thus the
training data sets in supervised classification tasks may be biased towards
brighter sources, or sources easier to detect in general, while the
classification algorithm is then used to predict the properties of fainter
objects observed in more recent surveys. Within the context of our work, the
sample of bright AGN and PSR labelled in a given catalog might not necessarily
reflect the characteristics of dim sources, which are the dominant part of the
UNC sources we want to classify (see e.g. Fig. 1). This effect could bias our
classifier towards the characteristics of the brighter sources. With the
cross-match data set of the 3FGL and the 4FGL-DR2 catalogs we can test the
performance of our classifier to predict the actual labels of the UNC sources.

Additionally, the cross-matching between the 3FGL and 4FGL-DR2 catalogs allows
us to evaluate the behaviour of our classifier with an increasing amount of
training data. We train ten networks for each architecture on 60, 70, 80, 90
and 100 percent of the labelled sources in the 3FGL catalog, and use the
cross-matched data set of the 4FGL-DR2 catalog as test set.

Note that the data format of the 3FGL catalog is somewhat different from the
4FGL-DR2 catalog. First, the significance for the Flux_History of the gamma-
ray flux is not given. In addition, the Flux_History is not given in yearly
bins, but in monthly bins. Finally, given the lower statistics, the Flux_Band
is only provided for 5 energy bins. So the input dimensions are $(5,2)$ for
the energy spectrum including the significance, and $(48,1)$ for the timing
information (four years with 12 monthly bins each).

In table 4 we report the results of the 3FGL and 4FGL-DR2 catalog cross-
matching (with PSR as positive class) for the Accuracy and AUC measures
defined in Sec. 3. The error corresponds again to the standard deviation of
the ten trainings. The accuracies of our neural network classifiers are
comparable to the accuracies stated in Luo et al. (2020) for various machine
learning algorithms trained on selected features. By using 100 percent of the
sources in the 3FGL catalog to train our classifiers, we achieve an accuracy
of $94.73\pm 0.46$ with the RNN, and very similar values for the DNN. This
accuracy is lower with respect to the one we find when using the 4FGL-DR2
catalog for training and testing (accuracies of $97.59\pm 0.35$% and $97.31\pm
0.36$% for the DNN and RNN, respectively, see Table 2), or if we compute the
accuracy within the 3FGL catalog without using the cross-match sources. In
fact, for the training on 70% of 3FGL sources we get an accuracy of $96.99\pm
0.50$% when we test on the remaining 30% of labelled sources. This is likely
due to the sample selection bias, as the training set is biased towards the
brightest sources, while the cross-match set is composed of dim sources, which
might not follow exactly the same distribution in feature space. A possible
way to address this issue when classifying based on features is presented in
Sec. 5.

Table 4: Results for the cross-match classification of the 3FGL and 4FGL-DR2 catalogs. Architecture | Measure | 70% Training set | 100% Training set
---|---|---|---
DNN | Acc [%] | $94.63\pm 0.42$ | $94.42\pm 0.44$
AUC [%] | $92.69\pm 2.31$ | $92.96\pm 1.48$
RNN | Acc [%] | $94.05\pm 0.67$ | $94.73\pm 0.46$
AUC [%] | $94.84\pm 1.76$ | $95.13\pm 0.77$

We do not observe a significant increase in performance when increasing the
training set. That is somewhat surprising, as we expect a stronger dependence
of the performance of our neural networks on the amount of training data than
in the case of the standard machine learning algorithms used in previous
analyses. It is likely that larger increases of the amount of training data
are needed in order to observe a significant trend in the performance within
the uncertainties.

##  5 Feature selection based on autoencoders

One of the main problems when dealing with classification using machine
learning on small datasets is the tendency for overfitting of the training
data, resulting in bad generalisation to unseen data.444We indeed observe
overfitting also for our standard architectures, especially when looking at
the accuracy in the PSR class or the F1 score. Our standard DNN gives an F1
score of 81.73% on the test data and 92.28% on the training data. Similar
results are obtained for the RNN. Feature selection reduces the number of
inputs and can prevent the classification algorithm to learn training-set
specific information. In Luo et al. (2020) the authors apply a recursive
feature elimination procedure using a Random Forest (RF) algorithm as a
classifier. In this setup a RF classifier is trained on labelled data, and
each feature is assigned an importance score. The feature with the lowest
score is eliminated and another RF classifier is trained on the new subset.
This procedure is repeated until a stopping condition is fulfilled, which can
either be a number of remaining features or, as in Luo et al. (2020), a
performance ratio between two consecutive classifiers.

While this algorithm is selecting the best features for classification, it
might be affected by the sample selection bias (covariate shift). In fact, the
algorithm relies on labelled data, and the selected sample thus corresponds to
those sources that are observed with sufficient significance to be classified.
An alternative way to reduce the number of features used for classification is
feature reduction using an autoencoder (AE) (Pasquet et al., 2019), as we
outline below.

An AE is a neural network architecture that first encodes an input to a
reduced, latent dimension and then extracts the original input from this
encoding. The latent representation should then contain the most important
information about the input, and can thus be used for classification. Since no
labels are necessary to train the AE, this is an example of unsupervised
learning, and the AE can be trained on either the training set, the testing
set or both. By training the AE on the test set one selects features that are
most important to describe the test sources. Using these features for
training, the classifier can reduce the sample selection bias. Note that the
features are not selected according to the best classification performance,
but to encode the most relevant information about the sources.

Using the source features of the catalog as input to the AE requires some
additional preprocessing. First of all, there are some features that are not
measured for all sources. We adopt the approach used in Luo et al. (2020),
where features with missing values in more than five percent of the sources
are removed. We further remove those sources where more than five of the
remaining features are empty. The remaining missing values are replaced by the
mean value of all classified sources. This procedure leaves us with 40
features.

We use 70% of the remaining sources for training and 30% for testing. Each
feature is normalised such that the mean across the training set is zero and
the standard deviation is one. After training the autoencoder on the test set,
a reduced set of features is extracted from the latent dimension and used in
the final classifier. Our autoencoder consists of three fully connected hidden
layers, the second of which provides the encoding. The first and third layer
consist of 8 nodes and apply PReLU activation. We find that training on the
encoding allows a simple fully connected classifier with one hidden layer of
40 nodes to train longer without overfitting.

To estimate the performance and the uncertainty we train both networks ten
times, resampling the training and testing set before each training. Using the
40 features directly in a shallow classifier with one hidden layer of 40
nodes, we find an accuracy of $97.91\pm 0.48$% with an F1-score of $82.52\pm
3.92$% (for PSR as positive class). Using the same classifier on a five-
dimensional encoding of the features from the AE, we find $97.96\pm 0.26$%
accuracy and a F1-score of $84.42\pm 2.04$%.

While the overall accuracy is similar, we see that using the AE for feature
selection increases the performance stability. This is what we expect from the
reduction of features. Even if the performance remains below the performance
obtained with a recursive feature elimination procedure based on a Random
Forest classifier (Luo et al., 2020), we emphasise again that our AE is not
trained to select features separating the two classes and might thus be less
affected by the sample selection bias. Exploring the potential of autoencoders
for automated feature selection in more detail is left for future work.

##  6 Predictions for UNC sources in the 4FGL-DR2 catalog

We finally perform a two-step classification of the UNC sources in the 4FGL-
DR2 catalog. We first classify sources as either AGN or PSR, and then classify
the corresponding PSR sources as either MSP or YNG. We specify that a source
is classified as a PSR with confidence if the mean prediction of the ten
classifiers (obtained from ten training runs, see Sec. 3) for the PSR class is
above 0.9. The same holds for the AGN classification. In table 5 we show the
number of classified sources and an estimated accuracy for these
classifications. The accuracies are estimated as mean and standard deviation
of the accuracies obtained on the test sets for the ten classifiers. Results
in table 5 are reported for six different architectures (described in Sec.
3.2): DNN (RNN), DNN (RNN) adding the GLAT for the sub-classification of YNG
versus MSP, and DNN (RNN) Band GLAT, where we removed the Flux_History data
and use only the Flux_Band data in combination with the GLAT. As discussed in
Sec. 3.2, adding the latitude increases the performance in the MSP vs. YNG
classification. In table 5 we illustrate this effect for the classification of
UNC sources. Since we add the latitude only for the sub-classification of the
PSR sources, the classification of the UNC sources in terms of AGN and PSR is
not affected.

Table 5: Number of classified sources per class with an accuracy estimate from the performance on the labelled test sets.  Architecture | PSR (acc) | AGN (acc) | YNG (acc) | MSP (acc)
---|---|---|---|---
DNN | 78 ($92.75\pm 3.25$) | 1050 ($98.91\pm 0.32$) | 22 ($93.70\pm 5.20$) | 21 ($86.14\pm 4.79$)
DNN GLAT |  |  | 20 ($92.19\pm 5.26$) | 22 ($89.49\pm 4.49$)
DNN Band GLAT |  |  | 21 ($91.45\pm 5.66$) | 24 ($90.17\pm 4.93$)
RNN | 57 ($94.14\pm 2.64$) | 1076 ($98.86\pm 0.27$) | 5 ($97.94\pm 2.73$) | 0 (–)
RNN GLAT |  |  | 2 ($97.07\pm 2.79$) | 12 ($95.80\pm 5.81$)
RNN Band GLAT |  |  | 2 ($97.68\pm 4.61$) | 8 ($93.55\pm 6.71$)

The sky distribution of AGN, YNG and MSP candidates obtained with the DNN
classifier, together with the UNC sources in the 4FGL-DR2 catalog, is
illustrated in Fig. 4. The YNG and MSP sources in Fig. 4 are obtained by
adding the Galactic latitude as input to the classification. The map shows the
position of 4FGL-DR2 UNC sources, and candidate AGN and PSR in Galactic
coordinates in a Mollview projection. A complementary view of their latitude
distribution is presented in Fig. 5. In the left (right) panel we show the AGN
(PSR) Galactic latitude distribution for labelled (associated and identified)
sources in the catalog together with the distribution of source candidates
found by our two neural network architectures. The distribution of UNC sources
in the 4FGL-DR2 catalog is also reported in the left panel for comparison. As
expected from the observed distribution of labelled sources, PSR candidates
are clustered along the Galactic plane. AGN candidates are distributed more
isotropically, although slightly peaked at low latitudes where the majority of
UNC sources is located. This is compatible with an isotropic distribution of
AGN in the inner Galaxy. In fact, by summing the labelled and the candidates
sources, we find a similar value of about 350 to 400 AGN for each bin of
sin(GLAT)$\,<|0.25|$. Moreover, YNG pulsar candidates are concentrated at low
latitudes, while MSP candidates can be found also at high latitudes away from
the Galactic plane.

Figure 4: Sky distribution in Galactic coordinates of 4FGL-DR2 unclassified
sources (black dots). AGN candidates identified by our DNN classifier are
marked with red circles, while blue squares and diamonds indicate the YNG
pulsar and MSP candidates identified using our DNN together with the Galactic
latitude, respectively.

Figure 5: Galactic latitude distribution of AGN (left panel) and PSR (right
panel) of 4FGL-DR2 sources. In each panel, the distribution of AGN/PSR
labelled sources in the 4FGL-DR2 is reported together with the distribution of
source candidates as derived by our two neural network architectures. The
distribution of UNC sources is also shown in the left panel for comparison.

The DNN classifier predicts that among the 1679 UNC sources in the 4FGL-DR2
catalog, 1050 are candidates for AGN ($62$%), 78 for PSR ($4.6$%), while 551
remain unclassified. These 551 sources are mostly concentrated along the
Galactic plane (see Fig. 4), where the bright Galactic emission likely
prevents a precise measurement of their energy and time spectra. Indeed, more
than half of the 551 sources that remain unclassified have a low significance
(Sign_Avg=5-10) and some Flags555Sources have non-zero analysis flags when the
systematic uncertainties of features depend strongly on the uncertainties when
modelling the diffuse gamma-ray background, in the case of low source-to-
background ratio, source confusion etc., see Abdollahi et al. (2020).. Similar
results are found by the RNN, with 1076 AGN candidates ($64$%), 57 PSR
candidates ($3.4$%), and 546 sources remaining unclassified. We find that the
RNN classifies a smaller number of PSR candidate sources compared to the DNN,
however with higher estimated accuracy. Without adding the latitude
information, the RNN is not able to find any MSP, because the separation of
the two classes by the classifier is not strong enough to pass the threshold.
The percentage of AGN and PSR candidates is in line with similar estimates in
previous publications (Doert & Errando, 2014; Saz Parkinson et al., 2016), and
in fair agreement with the distribution among labelled sources in the 4FGL-DR2
catalog ($61$% AGN and $4.5$% PSR). The predicted labels for UNC sources
obtained with the two architectures are overall consistent: for 91.4% of the
sources classified as AGN or PSR by the DNN, the RNN obtains the same
classification. Note that the RNN finds only 57 PSR, and the score is thus
limited to 98.1%. While the two classifiers differ in the certainty of their
predictions, they still give consistent results. If we modify the threshold of
the RNN to a lower value of 0.8 (0.7) we obtain an agreement of 97.7% (99.2%).
In fact, we verified that all sources identified as PSR or AGN by the DNN are
also classified to belong to the same class by the RNN when reducing the
threshold to 0.5.

Although our DNN and RNN architectures use only the energy and time spectra
for training, testing, and in order to provide a prediction for the UNC
sources, the distribution of relevant features can be analysed afterwards. In
Fig. 6 we show a variability-curvature plot for the labelled sources in the
4FGL-DR2 catalog, where the red contours refer to AGN and the blue dashed
contours to PSR, by extracting the features Variability_index and LP_SignCurv
(the curvature significance for the log parabola spectral fit). In this plane,
AGN and PSR are known to cluster in different regions, as PSR tend to be non-
variable on long time scales and have spectra exhibiting significant curvature
(see e.g. Abdo et al. (2010b)). In Fig. 6 we also show the distribution of
candidate AGN (red circles) and PSR (blue squares) sources identified by the
DNN. Even though the DNN is not trained to separate sources based on these
features (as done instead when using feature selection, see e.g. Saz Parkinson
et al. (2016); Luo et al. (2020)), the distribution of candidate sources
nicely follows the expected clustering in the variability-curvature plane,
confirming that the neural networks extract the relevant features from the
energy and time spectra. Similar results are obtained for the RNN.

The separation of the two classes from these features looks very
straightforward. However, the sources included in Fig. 6 are classified with
high certainty, and a clear separation is thus expected. Moreover, we chose
the features specifically to illustrate such a separation. For other features
the separations are not as obvious.

Table 6: List of MSP candidates selected from the UNC sources in the 4FGL-DR2 catalog by our DNN and RNN GLAT architectures. The last two columns contain a checkmark if the source is found also in the MSP candidate list of Saz Parkinson et al. (2016) (SP16) or Luo et al. (2020) (L20).  Source name | GLON [deg] | GLAT [deg] | Architecture | SP16 | L20
---|---|---|---|---|---
4FGL J0312.1-0921 | 191.496 | -53.37 | DNN/RNN | $\checkmark$ | -
4FGL J0953.6-1509 | 251.903 | 29.6028 | DNN/RNN | $\checkmark$ | -
4FGL J1120.0-2204 | 276.492 | 36.0559 | DNN/RNN | $\checkmark$ | -
4FGL J1221.4-0634 | 289.711 | 55.5178 | DNN/RNN | - | -
4FGL J1225.9+2951 | 185.427 | 83.7729 | DNN/RNN | - | -
4FGL J1400.0-2415 | 322.353 | 36.0036 | DNN/RNN | - | -
4FGL J1400.6-1432 | 326.982 | 45.0692 | DNN/RNN | - | $\checkmark$
4FGL J1627.7+3219 | 52.9838 | 43.2368 | DNN/RNN | $\checkmark$ | -
4FGL J2043.9-4802 | 351.629 | -38.2837 | DNN/RNN | - | -
4FGL J2112.5-3043 | 14.9039 | -42.4435 | DNN/RNN | $\checkmark$ | $\checkmark$
4FGL J2133.1-6432 | 328.748 | -41.2888 | DNN/RNN | $\checkmark$ | -
4FGL J2212.4+0708 | 68.7857 | -38.4799 | RNN | $\checkmark$ | $\checkmark$

Let us now discuss the classification of MSP vs. YNG pulsars in more detail.
Among the 21 sources classified as MSP by our standard DNN, 16 are contained
within the 22 and 24 MSP sources found when adding the latitude information
and only using the Flux_Band data, respectively. In total, these three
different classifications give 29 possible MSP candidates. As for the RNN, the
eight sources classified as MSP using the Flux_Band data in combination with
the GLAT are part of the 12 sources classified using RNN with energy, time
spectra and GLAT; 11 of these are included in those found by the DNN.
Considering the YNG class, the five candidate sources found by the RNN are all
contained within the YNG candidate sources obtained by the DNN.

In table 6 we provide a list of 11 promising MSP candidates identified by both
our DNN GLAT and RNN GLAT architectures, plus an additional source found with
high certainty by the RNN only. The last two columns illustrate the cross-
match of our list with previous analyses searching for MSP candidates (Saz
Parkinson et al., 2016; Luo et al., 2020).

An overlap of seven out of 13 MSP candidates is found with the list provided
in Saz Parkinson et al. (2016). The sources 4FGL J1221.4-0634 (3FGL
J1221.5-0632), 4FGL J1225.9+2951 (3FGL J1225.9+2953), 4FGL J1400.0-2415 (3FGL
J1400.2-2413), 4FGL J2043.9-4802 (3FGL J2043.8-4801) and 4FGL 1400.6-1432
(3FGL 1400.5-1437) are not present in the list of high-confidence candidates
in Saz Parkinson et al. (2016), most probably since these sources had a low
significance in the 3FGL; because of the low significance, they have not been
classified as probable AGN. Our overlap with the list provided in Luo et al.
(2020) is smaller, and only three sources are found in both lists. Since the
two lists are drawn according to different classification thresholds to define
high-confidence candidates, it is likely that a larger overlap could be found
taking more similar cuts. The first two and the last three sources in table 6
are also found in the list of Galactic candidates provided by Mirabal et al.
(2016). Future multiwavelength observation campaigns targeted towards these
candidate sources are required to finally clarify their association with MSP,
YNG or other objects. The full list of the DNN and RNN predictions for the
classes of UNC sources in the 4FGL-DR2 catalog is available at
https://github.com/manconi/agn-psr-nn-classification.

Figure 6: Sqrt(LP SignCurv) vs. Log(Variability Index) for labelled AGN (red
contours) and PSR (blue dashed contours) sources in the 4FGL-DR2 catalog. The
distribution of candidate AGN (red circles) and PSR (blue squared) obtained
from the DNN is superimposed.

In order to check the robustness of our MSP vs. YNG classifiers we performed a
further test. Instead of using the sources classified as PSR as a starting
point for our YNG vs. MSP classification, we start with the UNC sources
classified as AGN. The DNN tends to label these sources as mostly MSP (16 YNG,
414 MSP and 621 below threshold). This underlines the need for a proper
definition of an uncertainty for each classified source. This could for
example be achieved by transforming our classifiers into Bayesian neural
networks, which would provide an uncertainty for each source and might result
in high uncertainties for classes that are not used for training, see e.g.
Möller & de Boissière (2020). We leave this extension to future work.

##  7 Conclusions

In this paper we have presented a new approach to classify
unassociated/unidentified gamma-ray sources with deep neural networks. The
neural networks have been trained directly on the energy and time-dependent
photon fluxes as provided in the 4FGL-DR2 Fermi-LAT catalog, without relying
on specific, human-crafted features.

We have focused on the separation between Active Galactic Nuclei (AGN) and
Galactic pulsars (PSR), and the sub-classification of pulsars into young (YNG)
and millisecond pulsars (MSP). We have explored two different deep neural
networks, a fully connected dense neural network (DNN) and a recurrent neural
network (RNN). These neural network architectures have been shown to provide a
powerful classifier. They give compatible results, with a performance that is
comparable with previous analyses based on human-crafted features. For our
benchmark DNN we reach an accuracy of $97.59\pm 0.35$% for the classification
of AGN versus PSR, and of $87.25\pm 3.82$% for the sub-classification of
pulsars into YNG and MSP, when adding the Galactic latitude.

The accuracy of our classifiers has also been tested with a cross-match
procedure, where we have compared our predictions for the unlabelled sources
from the previous Fermi-LAT 3FGL source catalog with the classification
provided by the current 4FGL-DR2 catalog. We find an accuracy of $94.42\pm
0.44$% when testing against the cross-match data set. The lower accuracy with
respect to the accuracy found using the 4FGL-DR2 sources for training and
testing may point to a sample selection bias, as the training set of the 4FGL-
DR2 catalog may be biased towards the brightest sources.

We have proposed a method to address the sample selection bias when performing
source classification based on features. Using an autoencoder one can select
those features for classification that encode the most relevant information
about the sources, irrespective of the source labels and of the feature’s
classification performance. With such a technique we find an accuracy of
$97.96\pm 0.26$% for the classification of AGN versus PSR, comparable to what
we found using the energy and time dependent fluxes.

Finally, we have presented our predictions for the classification of
unlabelled sources in the 4FGL-DR2 catalog as either AGN or PSR, and of the
PSR sources as candidate YNG or MSP pulsars. Considering our benchmark DNN, we
identify 1050 candidate AGN and 78 candidate PSR sources. The distribution of
candidate sources follows the expected Galactic latitude distribution, and
most importantly, the clustering in the variability-curvature plane,
confirming that neural networks extract the relevant features from the energy
and time spectra. We have also provided a list of 12 MSP candidates selected
from unlabelled sources and discussed the overlap with results from previous
classifications based on other machine learning techniques.

The sources have been labelled as candidate AGN, PSR, YNG or MSP based on
their gamma-ray properties only. They should thus not be considered identified
sources, but rather as sources likely to be associated/identified in the
future. We encourage multiwavelength follow-up observations targeting the
candidate AGN and PSR, searching for example for radio and gamma-ray
pulsations, or spectral variability on larger time scales.

We emphasise that our method is very flexible and can easily be generalised to
include multiwavelength data on the energy and time spectra coming from
different observatories, as well as additional selected features. For example,
spectra can be complemented with data at higher energies coming from
forthcoming observations of Fermi-LAT sources with CTA (Acharya et al., 2018).

Up to now gamma-ray source classification techniques do not provide reliable
error estimates for the assignment of sources. Bayesian inference methods, for
example, may provide more information about the accuracy and the robustness of
the neural networks. We leave the exploration of this aspect to future work.

## Acknowledgements

We would like to thank Anja Butter, Mattia Di Mauro, Fiorenza Donato,
Alexander Mück and Tilman Plehn for insightful discussions. The work of TF is
supported by the Research Training Group GRK 2497 of the German Research
Foundation DFG. Computations were performed with computing resources granted
by RWTH Aachen University.

## Data Availability

The networks used in this paper are implemented using Tensorflow 2.1.0 (Abadi
et al., 2015) and the built in version of Keras (Chollet et al., 2015). The
reference dataset is the 4FGL-DR2 Fermi-LAT catalog (gll_psc_v23.fit)
available at: https://fermi.gsfc.nasa.gov/ssc/data/access/lat/10yr_catalog/.
The public list of LAT-detected gamma-ray pulsars used to build a cross-match
between this list and the 4FGL-DR2 catalog is available at
https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-
Detected+Gamma-Ray+Pulsars.

The complete list of the neural network predictions for the classes of UNC
sources in the 4FGL-DR2 catalog obtained in this paper is available at
https://github.com/manconi/agn-psr-nn-classification.

## References

  * Abadi et al. (2015) Abadi M., et al., 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, https://www.tensorflow.org/
  * Abdo et al. (2009) Abdo A. A., et al., 2009, Science, 325, 848
  * Abdo et al. (2010a) Abdo A. A., et al., 2010a, The Astrophysical Journal Supplement Series, 188, 405
  * Abdo et al. (2010b) Abdo A. A., et al., 2010b, ApJS, 188, 405
  * Abdo et al. (2013) Abdo A., et al., 2013, Astrophys. J. Suppl., 208, 17
  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, Astrophys. J. Suppl., 247, 33
  * Acero et al. (2015) Acero F., et al., 2015, Astrophys. J. Suppl., 218, 23
  * Acero et al. (2016) Acero F., Ackermann M., Ajello M., Albert A., Baldini L., et al., 2016, ApJS, 223, 26
  * Acharya et al. (2018) Acharya B., et al., 2018, Science with the Cherenkov Telescope Array. WSP (arXiv:1709.07997), doi:10.1142/10986
  * Ackermann et al. (2012a) Ackermann M., Ajello M., Albert A., Allafort A., Atwood W. B., et al., 2012a, ApJS, 203, 4
  * Ackermann et al. (2012b) Ackermann M., Ajello M., Atwood W. B., Baldini L., Ballet J., et al., 2012b, ApJ, 750, 3
  * Ackermann et al. (2015) Ackermann M., Ajello M., Albert A., Atwood W. B., Baldini L., et al., 2015, ApJ, 799, 86
  * Atwood et al. (2009) Atwood W. B., Abdo A. A., Ackermann M., Althouse W., Anderson B., et al., 2009, ApJ, 697, 1071
  * Ball & Brunner (2010) Ball N. M., Brunner R. J., 2010, Int. J. Mod. Phys. D, 19, 1049
  * Ballet et al. (2020) Ballet J., Burnett T., Digel S., Lott B., 2020
  * Baron (2019) Baron D., 2019, Machine Learning in Astronomy: a practical overview (arXiv:1904.07248)
  * Becker et al. (2020) Becker I., Pichara K., Catelan M., Protopapas P., Aguirre C., Nikzat F., 2020, MNRAS, 493, 2981
  * Carleo et al. (2019) Carleo G., Cirac I., Cranmer K., Daudet L., Schuld M., Tishby N., Vogt-Maranto L., Zdeborová L., 2019, Reviews of Modern Physics, 91
  * Chawla et al. (2002) Chawla N. V., Bowyer K. W., Hall L. O., Kegelmeyer W. P., 2002, Journal of artificial intelligence research, 16, 321
  * Cheng (2018) Cheng T., 2018, Computing and Software for Big Science, 2
  * Chiaro et al. (2016) Chiaro G., Salvetti D., La Mura G., Giroletti M., Thompson D., Bastieri D., 2016, Mon. Not. Roy. Astron. Soc., 462, 3180
  * Chollet et al. (2015) Chollet F., et al., 2015, Keras, https://keras.io
  * Doert & Errando (2014) Doert M., Errando M., 2014, The Astrophysical Journal, 782, 41
  * Egan et al. (2017) Egan S., Fedorko W., Lister A., Pearkes J., Gay C., 2017, Long Short-Term Memory (LSTM) networks with jet constituents for boosted top tagging at the LHC (arXiv:1711.09059)
  * Englert et al. (2020) Englert C., Fairbairn M., Spannowsky M., Stylianou P., Varma S., 2020, Sensing Higgs cascade decays through memory (arXiv:2008.08611)
  * Fornasa & Sánchez-Conde (2015) Fornasa M., Sánchez-Conde M. A., 2015, Phys. Rep., 598, 1
  * Fraser & Schwartz (2018) Fraser K., Schwartz M. D., 2018, Journal of High Energy Physics, 2018
  * Guest et al. (2016) Guest D., Collado J., Baldi P., Hsu S.-C., Urban G., Whiteson D., 2016, Physical Review D, 94
  * Harding & Muslimov (1998) Harding A. K., Muslimov A. G., 1998, ApJ, 508, 328
  * Hinners et al. (2018) Hinners T. A., Tat K., Thorp R., 2018, The Astronomical Journal, 156, 7
  * Hui et al. (2020) Hui C., et al., 2020, Mon. Not. Roy. Astron. Soc., 495, 1093
  * Ishida (2019) Ishida E. E. O., 2019, Nature Astronomy, 3, 680
  * Kingma & Ba (2014) Kingma D. P., Ba J., 2014
  * Kovačević et al. (2020) Kovačević M., Chiaro G., Cutini S., Tosti G., 2020, Mon. Not. Roy. Astron. Soc., 493, 1926
  * Louppe et al. (2019) Louppe G., Cho K., Becot C., Cranmer K., 2019, Journal of High Energy Physics, 2019
  * Luo et al. (2020) Luo S., Leung A. P., Hui C. Y., Li K. L., 2020, Monthly Notices of the Royal Astronomical Society, 492, 5377–5390
  * Mirabal et al. (2012) Mirabal N., Frias-Martinez V., Hassan T., Frias-Martinez E., 2012, Mon. Not. Roy. Astron. Soc., 424, L64
  * Mirabal et al. (2016) Mirabal N., Charles E., Ferrara E. C., Gonthier P. L., Harding A. K., Sánchez-Conde M. A., Thompson D. J., 2016, The Astrophysical Journal, 825, 69
  * Möller & de Boissière (2020) Möller A., de Boissière T., 2020, MNRAS, 491, 4277
  * Naul et al. (2018) Naul B., Bloom J. S., Pérez F., van der Walt S., 2018, Nature Astronomy, 2, 151
  * Padovani et al. (2017) Padovani P., et al., 2017, Astron. Astrophys. Rev., 25, 2
  * Pasquet et al. (2019) Pasquet J., Pasquet J., Chaumont M., Fouchez D., 2019, A&A, 627, A21
  * Romani (2014) Romani R. W., 2014, Science, 344, 159
  * Salvetti et al. (2017) Salvetti D., Chiaro G., La Mura G., Thompson D. J., 2017, Monthly Notices of the Royal Astronomical Society, 470, 1291–1297
  * Saz Parkinson et al. (2016) Saz Parkinson P. M., Xu H., Yu P. L. H., Salvetti D., Marelli M., Falcone A. D., 2016, The Astrophysical Journal, 820, 8
  * The Fermi-LAT Collaboration (2019) The Fermi-LAT Collaboration 2019
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, PASP, 107, 803