# BlaST - A Machine-Learning Estimator for the Synchrotron Peak of Blazars

T. Glauch  T. Kerscher  P. Giommi

###### Abstract

Active Galaxies with a jet pointing towards us, so-called blazars, play an
important role in the field of high-energy astrophysics. One of the most
important features in the classification scheme of blazars is the peak
frequency of the synchrotron emission ($\nu_{\rm peak}^{S}$) in the spectral
energy distribution (SED). In contrast to standard blazar catalogs that
usually calculate the $\nu_{\rm peak}^{S}$manually, we have developed a
machine-learning algorithm - BlaST - that not only simplifies the estimation,
but also provides a reliable uncertainty evaluation. Furthermore, it naturally
accounts for additional SED components from the host galaxy and the disk
emission, which may be a major source of confusion. Using our tool, we re-
estimate the synchrotron peaks in the Fermi 4LAC-DR2 catalog. We find that
BlaST improves the $\nu_{\rm peak}^{S}$ estimation especially in those cases
where the contribution of components not related to the jet is important.

###### keywords:

Machine learning , galaxies: jets , galaxies: active , BL Lacertae objects:
general , methods: data analysis , astronomical databases: miscellaneous

††journal: Astronomy & Computing

[inst1]organization=Technische Universität München, Physik-Department,
addressline=James-Frank-Str. 1, city=Garching bei München, postcode=D-85748,
country=Germany

[inst2]organization=Institute for Advanced Study, Technische Universität
München, addressline=Lichtenbergstrasse 2a, city=Garching bei München,
postcode=D-85748, country=Germany

[inst3]organization=Center for Astro, Particle and Planetary Physics (CAP3),
addressline=New York University Abu Dhabi, city=Abu Dhabi, postcode=PO Box
129188, country=United Arab Emirates

[inst4]organization=Associated to Agenzia Spaziale Italiana, ASI,
addressline=via del Politecnico s.n.c., city=Roma, postcode=I-00133,
country=Italy

[inst5]organization=Heidelberg University, Institute of Environmental Physics,
addressline=Im Neuenheimer Feld 229, city=Heidelberg, postcode=D-69120,
country=Germany

The synchrotron peak is an important measure for the classification of blazars

BlaST uses machine learning to determine the peak in the spectral energy
distribution

Consistent & fast estimation that accounts for emission from the host galaxy

Re-estimation of all synchrotron peaks in the Fermi 4LAC catalog

##  1 Introduction

Blazars are among the most powerful and variable persistent sources in the
Universe. These remarkable objects are numerically a tiny minority compared to
stars, galaxies, and even quasi stellar objects. Still, they are the subject
of very active contemporary research, because of their dominant role in
current high-energy astrophysics, and in the newly emerging field of multi-
messenger astronomy. Blazars are a special type of active galactic nuclei –
galaxies hosting a very luminous core – (see Padovani et al., 2017, for a
recent review) with a jet of relativistic plasma pointing towards the Earth.
Due to this property blazars are strong emitters across the entire
electromagnetic spectrum, from radio waves up to TeV $\gamma$-rays. Blazars
are sub-classified on the basis of their optical spectrum as Flat-Spectrum
Radio Quasars (FSRQs) and BL Lacertae objects (or BL Lacs), with FSRQs showing
broad emission lines, and BL Lacs being featureless or displaying very weak
emission lines (Falomo et al., 2014) The sketch shown in Figure 1 illustrates
how the textbook example spectral energy distribution (SED) of blazars is
characterized by two broad bumps, which are thought to be the signature of the
non-thermal radiation emitted by the relativistic jet. The first bump, which
peaks between the far infrared and the X-ray bands, depending on the blazar
type, is commonly attributed to synchrotron radiation emitted by highly-
accelerated electrons and potentially protons in the magnetic field of the
jet. The second, more energetic, bump peaks in different parts of the
$\gamma$-ray band and has been attributed to different mechanisms, such as
inverse Compton scattering, decay of mesons or proton synchrotron radiation.
Additionally, the SED of a Blazar can also include other features such as
thermal radiation from the host galaxy (infrared bump or stellar emission),
and radiation from the accretion disk around the central black hole (blue
bump) and from the broad-line region (e.g. Giommi et al., 2012b). Fig. 2 shows
the real SEDs of two bright and well-monitored blazars (MRK 501 and 3C 279)
characterised by largely different energy distributions. For more examples of
recent well populated SEDs of blazars of all types; see Giommi et al.
(2021)111see also https://openuniverse.asi.it/blazars/swift/.

The frequency where the synchrotron component peaks in the observer frame
222In order to account for the redshift a factor of (1+z) should be used, if
the redshift $z$ is known. ($\nu_{\rm peak}^{S}$ ) is believed to carry
important information about the acceleration power of a source. Hence, blazars
are historically separated into low, intermediate and high-energy peaked
sources or LBLs, IBLs and HBLs, see (Padovani and Giommi, 1995), depending on
the rest-frame position of their $\nu_{\rm peak}^{S}$ according to the
following definitions,

  * 1.

LBLs: $\nu_{\rm peak}^{S}$ $<10^{14}$ Hz [$<$ 0.41 eV],

  * 2.

IBL: $10^{14}$ Hz$~{}<$ $\nu_{\rm peak}^{S}$ $<10^{15}$ Hz [0.41 – 4.1 eV],

  * 3.

HBLs: $\nu_{\rm peak}^{S}$ $>10^{15}$ Hz [$>$ 4.1 eV].

Abdo et al. (2010) modified the naming to LSP, ISP and HSP, maintaining the
same frequency (or energy) limits, to describe the SED of FSRQs as well. We
feel that the SED classification should be independent from the one based on
the presence of broad lines in the optical spectrum. In the following we adopt
the original nomenclature for all blazars.

Since then the $\nu_{\rm peak}^{S}$ parameter has been acquiring increasing
importance and is now widely used in the literature. Recent studies suggest
that the position of the peak is more than just a feature in the SEDs, but
that objects with different $\nu_{\rm peak}^{S}$ values belong to separate
source populations with different intrinsic properties (Giommi and Padovani,
2021).

Figure 1: A schematic view of the SED of different types of blazars. The
vertical axis shows the energy flux against the emission frequency (or energy)
on the horizontal axis. The peak of the synchrotron component ($\nu_{\rm
peak}^{S}$ ) spans a wide range of frequencies, from the far infrared in LBL
objects (red curves) to the X-ray band in HBL sources (blue curves).

Figure 2: Examples of well populated (time-integrated) SEDs of blazars with
high $\nu_{\rm peak}^{S}$ (MRK 501, top panel) and low $\nu_{\rm peak}^{S}$
(3C 279, lower panel), corresponding to the blue and red lines in the
schematic representations of Fig.1. Note the large flux variability in both
cases.

A number of recent works rely on $\nu_{\rm peak}^{S}$ values that have been
estimated in a variety of ways using different tools and data sets, to
identify interesting sub-samples of blazars that are used for statistical
investigations. Examples include the $\gamma$-ray 4FGL-DR2 sample (Abdollahi
et al., 2020; Ballet et al., 2020), the radio selected MOJAVE sample (Homan et
al., 2021), the multi-frequency selected 3HSP catalogue (Chang et al., 2019b),
the sample of blazars frequently observed by the Neil Gehrels Swift
observatory (Gehrels et al., 2004; Giommi et al., 2021), and the sample of IBL
blazars of Giommi and Padovani (2021).

Recently, IceCube announced the first evidence for neutrino emission from a
$\gamma$-ray bright IBL-type blazar, TXS 0506+056 (Aartsen et al., 2018b, a;
Padovani et al., 2018), indicating also the production of cosmic-rays, i.e.,
protons and heavier nuclei, with energies at least up to $\sim$10 PeV. Other
studies indicate that especially those blazars with $\nu_{\rm peak}^{S}$above
$10^{14}$ Hz are connected to neutrino emission (Giommi et al., 2020).

A reliable estimation of $\nu_{\rm peak}^{S}$ is important for the
understanding of the emission mechanisms in blazars as it is an observational
parameter that strongly constrains theoretical models (e.g. Ghisellini et al.,
1985; Giommi et al., 2012b; Stathopoulos et al., 2021). Precise values of
$\nu_{\rm peak}^{S}$ of the blazars included in unbiased samples, may also
help settling long-standing issues like the existence of the so called blazar
sequence, originally introduced by Fossati et al. (1998) and debated over the
years in a number of papers (e.g. Ghisellini and Tavecchio, 2008; Giommi et
al., 2012a; Cerruti et al., 2017; Keenan et al., 2021; Prandini and
Ghisellini, 2022).

Several tools exist to build SEDs suitable for the estimation of the $\nu_{\rm
peak}^{S}$ . The most frequently used are the SSDC SED
builder333https://tools.ssdc.asi.it/SED/, the VOU-Blazars
tool444https://hub.docker.com/r/chbrandt/voublazars/,555https://github.com/ecylchang/VOU_Blazars
(Chang et al., 2020), the SED data
service666https://ned.ipac.caltech.edu/forms/photo.html, and the Vizier
photometry tool777http://vizier.u-strasbg.fr/vizier/sed/, each providing
different multi-frequency data sets, different conversion methods from
integrated flux measurements to $\nu$f$({\nu})$ fluxes, and different ways of
de-reddening fluxes at frequencies where absorption in our Galaxy is not
negligible.

Despite the relative abundance of SED data, the precise determination of
$\nu_{\rm peak}^{S}$ is not a straightforward task as it is subject to a high
level of heterogeneity. This has several reasons, including the different
methods and tools used to build the SEDs, different multi-frequency data sets,
large flux and spectral variability, different approaches in the de-reddening
of the observed flux at optical, UV, and soft X-ray frequencies, etc. Another
complication in the estimation of the $\nu_{\rm peak}^{S}$ is the frequent
presence of spectral components that are not related to the jet like the blue
bump, the dusty torus, and the stellar emission of the Host galaxy. These
features are often difficult to identify by non-experts and are not taken into
account in a uniform way by different teams, potentially leading to large
errors in the estimations of the $\nu_{\rm peak}^{S}$ .

To make future works more accurate and comparable and to remove the
inhomogeneities connected to human components in the interpretation of the
data, it is useful to homogenise the way the $\nu_{\rm peak}^{S}$ is
calculated. In this paper we propose a new method to estimate the $\nu_{\rm
peak}^{S}$ based on a machine learning tool that is trained on 3,793 blazar
SEDs with the $\nu_{\rm peak}^{S}$ranging between $10^{12}$ Hz and $10^{18}$
Hz. In contrast to the human estimation, the machine learning algorithms is
able to include all the information of the SED and to better recognize
features that are not related to the synchrotron emission, namely the emission
from the accretion (blue bump), the dust in the galaxy (infrared bump) and
emission from the broad-line region. Furthermore it allows the determination
of a confidence interval for the estimated $\nu_{\rm peak}^{S}$.

This paper is structured as follows: First, we describe the training dataset
of blazars, then we explain the different machine learning algorithms that
have been tested as well as their performance, and finally we compare and
validate the best algorithms on the 4LAC-DR2 catalog (Lott et al., 2020).

##  2 Data sample of blazars

Currently, the largest collection of blazars is the master
list888https://openuniverse.asi.it/OU4Blazars/MasterListV2/ compiled as part
of the Open Universe for Blazars initiative. This table includes over 6,000
blazars assembled combining the BZCAT 5th edition (Massaro et al., 2015,
2016), the sample of blazars listed in the Fermi LAT catalogues (Ackermann et
al., 2015a, b; Ballet et al., 2020), and the 3HSP sample of high $\nu_{\rm
peak}^{S}$ blazars (Chang et al., 2019c, a). The master list includes blazars
discovered in radio, X-ray, and $\gamma$-ray surveys, as well as objects
discovered for their multi-frequency SED properties, and therefore it covers
blazars of all known types.

The dataset that has been used to train our $\nu_{\rm peak}^{S}$ estimator -
called BlaST 999Acronym for Blazar Synchrotron Tool or Blazar SED Tool \- is a
subsample of the blazars master list, selected on the basis of the
availability of sufficient multi-frequency data and to ensure that all blazars
types (LBL, IBLs and HBLs) and data combinations (jet emission plus other non-
jet related components) are adequately represented. The SED of each blazar in
the sample was assembled using the latest version (V1.94) of the VOU-Blazars
tool (Chang et al., 2020), which retrieves multi-frequency data from 71
catalogues and spectral databases from different on-line services using
Virtual Observatory101010https://ivoa.net/ protocols. Besides including all
the classical multi-frequency catalogs, this version of VOU-Blazars also
includes recent radio surveys, such as the VLASSQL catalog (Gordon et al.,
2020) and the most up to date catalogs of the Swift (Evans et al., 2020;
Giommi et al., 2019), XMM (Webb et al., 2020), and Fermi-LAT 4FGL-DR2 (Lott et
al., 2020) detectors. Once the data is downloaded, VOU-Blazars automatically
converts it to homogeneous SED units before de-reddening the optical to soft-
X-ray measurements to remove the effects of absorption in our Galaxy. The de-
reddening is especially important for sources close to the Galactic plane and
automatically done by standard tools as VOU-Blazars or the SSDC SED
Tool111111https://tools.ssdc.asi.it/SED/. In case other tools are used to
retrieve the SED, a user should ensure that the data points are correctly de-
reddened before using BlaST.

The SEDs assembled in the previously described procedure contain data from
different time periods. This is unavoidable as measurements are rarely
performed simultaneously. Hence, the $\nu_{\rm peak}^{S}$predicted by the tool
is a time-averaged one as typically used in the literature.

The value of $\nu_{\rm peak}^{S}$ of each blazar in the sample was then
estimated by the authors via a manual process, also by fitting the data to
polynomial functions, taking great care in avoiding the SED components not
related to the emission from the jet, such as the host galaxy and the blue
bump. This value is used as the ground truth to train our machine learning
algorithms as described in the following sections.

##  3 Developing a machine-learning estimator

###  3.1 Data preparation

For each SED in the training sample, we download the available multi-frequency
flux measurements using the VOU-Blazars tool to a human-readable file. In rare
cases VOU-Blazar returns non-physical sampling points such as zero fluxes or
best-fit frequencies that are outside their own error margins. We remove those
points. Additionally, all values (fluxes and frequencies) are converted into
decadic log space as it proved to be beneficial for the training procedure
later. 700 of the total 3,793 dataset entries were put aside as test set for
later comparison.

Figure 3: Bins used for the data augmentation process. Edges have been chosen
for uniform number of data points and thus do not have a unique width. From
the initial 34 bins, only 26 show reduction in the bias after being augmented,
see Figure 4. The remaining 8 biased bins are hatched. Figure 4: Comparison of
the bias in the raw dataset (left) versus the augmented dataset (right). The
bin IDs on the horizontal axis correspond to the (non-uniform) bins in Figure
3 in increasing order. The darker the color, the more data points are
available in the respective bin. The augmented dataset shows a visible
reduction in overall bias, but still suffers from it in certain frequency
ranges, especially between $\sim 10^{10}$ Hz and $\sim 10^{12.5}$ Hz.

In order to evaluate the distribution of available frequencies in the training
sample, the data is sorted in 34 frequency bins according to Figure 3. To
define the frequency bins, we have started with a fine binning and merged
sparse bins until we achieve a uniform number of data point across the bins.
We have validated that a further reduction of bins would reduce the final
performance of the estimator. Clearly, the dataset generated in this way has a
significant bias, as some emission frequencies are dominantly present for
specific $\nu_{\rm peak}^{S}$ , see Figure 4.

The bias was taken care of by augmenting the data set through oversampling:
The samples were split into 32 uniform bins based on their $\nu_{\rm
peak}^{S}$ (vertical axis in Figure 4). Each of these $\nu_{\rm peak}^{S}$
bins were then filled with randomly-selected copies of its members to achieve
90 samples per bin. To make these copies distinct, for each copy 15 (not
necessarily empty) frequency bins were selected, through random sampling with
replacement, and then deleted. The probability for a frequency bin to be
chosen is proportional to its overall abundance to the power of four, i.e.,
more prominent frequencies are much more likely to be deleted. Copies with
less than 5 remaining non empty bins were assumed to no longer hold enough
information and were thus discarded. Through this augmentation process the
amount of SEDs grew to a total of 10,233, while on average 12 bins were
deleted. Figure 4 shows that this procedure significantly reduces the biases
of the dataset.

By construction, the number of available measurements varies strongly between
the SEDs. To create a fixed size representation, as needed by the machine-
learning algorithms, three different approaches were originally tested: Simple
binning, auto-encoder (AE) as depicted by Goodfellow et al. (2016), and
compressed sensing (CS) as introduced by Candes et al. (2006).

AE are neural networks consisting of the encoder transforming the input in a
lower dimensional representation followed by the decoder which reconstructs
the original input from it. After minimizing the reconstruction error, the
decoder part can be used to compress the input while retaining all the
essential information.

CS exploits the sparse nature most signals inhibit, e.g. music consist of
specific sine waves and is thus sparse in the Fourier basis, by finding the
coefficients in such a basis that minimizes the absolute error to the original
signal while penalizing the sum of absolute coefficients. This algorithm is
widely known as Lasso (Tibshirani, 1996) and produces a sparse set of
coefficients. These coefficients can finally be used to reconstruct the
complete signal.

All three methods require a binned input. The simple binning approach uses the
non-uniform bin width in Figure 3 in order to achieve uniform number of data
points, while the autoencorder and compressed sensing use a fixed bin width of
0.1 in decadic log space. Bins that are biased after augmentation are ignored
in the case of the simple binning and the AE. Additionally, AE ignores also
empty bins. Compressed sensing, on the contrary, uses all bins. The reason is
that CS completes the sparse SED data and thus eliminates the correlation
between missing bins and $\nu_{\rm peak}^{S}$. For the simple binning case, we
have overall 34 bins. However, 8 of these bins still show a bias after
augmentation and are thus dropped, leading to 26 bins. In each case the bin
value is the mean of all the measurements within it.

Since CS exploits the sparsity of a signal it requires a custom vector space
basis of SEDs. We construct this basis through a process known as dictionary
learning (Mairal et al., 2009) using 500,000 synthetic SEDs generated with
naima (Zabalza, 2015).

For the AE, both the en- and decoder consist of a two layer, bidirectional
gated recurrent unit neural network, while the decoder additionally has a
final feed forward layer. The encoder iterates only over non empty, unbiased
bins in ascending order with frequency and flux as input, compressing into an
encoded vector of size 64 as its final output. Besides the encoded vector, the
decoder predicts one flux value for each frequency with non-empty input. The
squared difference to the original value is used as a loss metric. The size of
the encoded vector was chosen in a hyperparameter search such that the loss is
minimal.

In summary, each of the methods follows a different approach to cope with the
sparsity and variable length of the input: Binning uses non uniform bin widths
and deliberately accepts empty bins, while trying to reduce their occurrence.
CS can handle smaller bins and still avoid empty ones by filling in missing
data. AE also uses smaller bins, but ignores the empty ones altogether by
providing an encoded vector.

###  3.2 Training & Architecture

In total, three different and widely known machine learning algorithms have
been examined: Random Forest (RF), Gradient Boosting (GB) and Neural Networks
(NN). For RF and GB the implementation provided by the python package _scikit-
learn_ (Pedregosa et al., 2011), while for the NN _pytorch_ (Paszke et al.,
2019) was used. For the neural networks the input data - as described in the
previous section - has been scaled to unit variance with zero mean. We only
considered algorithms which support prediction intervals (PI). For GB this can
be achieved by utilizing the quantile loss function. Likewise a method for
random forest was introduced by Meinshausen (2006) and one for NN by
Lakshminarayanan et al. (2017).

Figure 5: Distribution of estimation error with a specific PI noted as 1$\sigma$ alongside the expected Gaussian distribution as solid line. Random Forest | Gradient Boosting | Neural Network
---|---|---
Num. Trees | Num. Trees | Layer1
| Learning rate | Layer2
| Max. tree depth | Layer3
| Min. samples leaf | Layer4
| Min. samples split | Learning Rate
|  | Weight Decay
Table 1: List of hyperparameters that are optimized for the respective
algorithms. For each fully-connected layer in the neural network case the
number of neurons and the dropout rate are optimized. For more information we
refer to the respective documentation in (Pedregosa et al., 2011) (RF, BG) and
(Paszke et al., 2019) (NN).

Finally, a distinct model has been trained on each combination of machine-
learning algorithm (RF, GB, NN) and feature set (Binning, AE, CS) totalling in
9 $\nu_{\rm peak}^{S}$estimators. For each estimator the corresponding
hyperparameters have been determined by an automatic search (Liaw et al.,
2018). A list of the hyperparameters in each model is found in Table 1. In
particular, RF and GB were tuned using a 5-fold cross validation, while NN
used a validation set consisting of 1,000 samples put aside beforehand. Using
the previously created test set their performances are evaluated and listed in
Table 2(c).

| Binning | CS | AE
---|---|---|---
RF | $0.300^{+0.313}_{-0.170}$ | $0.376^{+0.324}_{-0.189}$ | $0.400^{+0.310}_{-0.200}$
GB | $0.469^{+0.427}_{-0.257}$ | $0.356^{+0.330}_{-0.195}$ | $0.403^{+0.333}_{-0.220}$
NN | $0.265^{+0.242}_{-0.133}$ | $0.314^{+0.257}_{-0.173}$ | $0.316^{+0.255}_{-0.167}$
(a) Median Absolute Error

| Binning | CS | AE
---|---|---|---
RF | $2.583^{+1.325}_{-0.683}$ | $2.301^{+1.084}_{-0.501}$ | $3.100^{+0.973}_{-0.900}$
GB | $4.675^{+0.445}_{-1.176}$ | $4.264^{+0.544}_{-1.158}$ | $3.804^{+1.058}_{-0.835}$
NN | $1.687^{+0.307}_{-0.197}$ | $2.074^{+0.307}_{-0.197}$ | $2.085^{+0.417}_{-0.261}$
(b) Median PI width

| Binning | CS | AE
---|---|---|---
RF | $2.992^{+0.982}_{-1.040}$ | $2.666^{+0.734}_{-0.816}$ | $3.195^{+0.905}_{-0.993}$
GB | $4.311^{+0.823}_{-0.809}$ | $3.943^{+0.893}_{-0.816}$ | $3.932^{+1.034}_{-0.933}$
NN | $1.910^{+0.132}_{-0.415}$ | $2.277^{+0.270}_{-0.635}$ | $2.302^{+0.246}_{-0.472}$
(c) Mean Interval Score (Lower better)

Table 2: Summary of the model performances on the test set

in terms of the the median absolute prediction error (top), the median PI
(95%) width (middle) and the mean interval score (bottom) including their 25%
and 75% quantile, respectively. The interval score as described by Gneiting
and Raftery (2007) is a strictly proper scoring rule for each prediction
rewarding narrow PI while penalizing if the true values is not covered by it
and shows a unique minimum for a perfect point like prediction.

The best performance is reached for the NN trained on binned data. Thus, this
method is used for the final version of the tool. It will be discussed in
detail in the following. It is noteworthy, however, that the random forest on
binned data was not only fast to train, but also showed very good performance
”out-of-the-box”, i.e., without any hyperparameter tuning.

For the final model the training procedure is slightly altered: The test set
is replaced by bagging without replacement, i.e., similar to cross-validation
the dataset is divided into 5 splits from which 5 new datasets are created
each time using 4 splits and leaving 1 out. Subsequently, 5 neural-network
models with different initializations were trained for each of these five new
datasets. In total, we therefore train an ensemble of 25 networks. This
ultimately allows to create genuine estimates for any SED by either using the
complete ensemble for truly unseen data, or instead using the sub-ensemble
consisting of only the five NN which haven’t seen that particular SED during
the training.

Note that, after generation, each of the 5 new datasets has been augmented
independently from another using the methods describe above to avoid a set
consisting of only augmented samples. This increased their initial size on
average from 750 to 2,600 samples and thus the bagged sets’ to 10,400 samples,
from which 1,500 each were put aside as validation set.

As mentioned earlier, the binned input is scaled to unit variance with zero
mean. However, the zero flux of empty bins after this transformation are far
away from zero and thus may harm the networks performance by saturating the
non-linearities. To aid the network a mask indicating empty bins with zero and
filled ones with one is attached, thus doubling the number of inputs to 52.

All NN of the ensemble share the same architecture of 4 hidden layers followed
by a single linear layer with the sizes 52-152-80-72-45-2. Each hidden layer
thereby consists of a linear layer, followed by batch normalization, dropout
and ReLU as non-linear activation function. The drop-out rates are 11.62%,
14.95%, 2.56% and 3.21%, respectively. Following Lakshminarayanan et al.
(2017) we assume the negative log-likelihood of a Gaussian distribution as
loss function. Hence our model estimates the mean and standard deviation of
this distribution. To enforce the latter to be positive, a softplus is applied
in the final layer. In order to combine the predictions ($\hat{\mu}_{m}$,
$\sigma_{m}$) of the $M$ neural nets, we use the methods presented in
(Lakshminarayanan et al., 2017). Hence, for the mean

$\hat{\mu}(x)=\frac{1}{M}\sum_{m=1}^{M}\hat{\mu}_{m}(x)$ (1)

and the variance

$\sigma^{2}(x)=\frac{1}{M}\left(\sum_{m=1}^{M}\sigma_{m}^{2}(x)+\hat{\mu}_{m}^{2}(x)\right)-\hat{\mu}^{2}(x).$
(2)

Figure 5 shows that the Gaussian is indeed a good description for the variance
of the prediction.

The hyperparameters, including the depth, were reused from the process of the
best-algorithm selection above.

Each model of the ensemble has been trained on their respective bagged set
independently from another using Adam with a learning-rate of $1.840\times
10^{-3}$ and a weight decay of $1.635\times 10^{-5}$. If the performance on
the validation set plateaued for 10 iterations, the learning rate was reduced
to a tenth. The models were trained for up to 800 epochs, but only the one
with the best validation loss was kept - in most cases after a few hundred
epochs.

###  3.3 Validation & Performance

Figure 6: Histogram of the prediction on the data set for validation. The
solid black line shows the median prediction, while the dashed ones show the
10% and 90% quantiles. The expectation for a perfect prediction is given by
the grey line for comparison.

The splitting method of the previous section allows to use the whole training
data to determine BlaST’s performance, as it ensures that there is always an
ensemble of 5 NN which haven’t been trained on a specific SED. Hence, each SED
has an independent $\nu_{\rm peak}^{S}$ estimate.

Figure 6 shows the precision of the algorithm by comparing the ground truth
against the prediction. In the relevant range of the parameter space the
median prediction is well compatible with the expectation. Hence, we consider
the estimator unbiased. Small deviations are visible towards the upper and
lower end of the catalog $\nu_{\rm peak}^{S}$. Those can be attributed to
lacking data in the far-infrared and X-ray, the algorithm itself and imprecise
$\nu_{\rm peak}^{S}$ determinations in the training dataset. The latter are
discussed in more detail in section 4. In addition, deviations at the
boundaries can be caused through threshold effects. The network is penalized
if it predicts values below the lower $\nu_{\rm peak}^{S}$ limit of the
training dataset imposing essentially a one-sided Gaussian loss function at
the boundary. This leads to an overall over-estimations of the $\nu_{\rm
peak}^{S}$ and vice versa an underestimation for the upper $\nu_{\rm
peak}^{S}$limit. Overall, the median absolute error between prediction and
ground truth is $0.260$.

Figure 7: Prediction Intervals (PI) in blue alongside the actual ground truth
marked as orange dots. To make the plot more readable the PI’s edges and the
ground truth are shifted by their center, i.e. a perfect score would be a
straight line at zero. On the horizontal axis the samples are ordered by their
PI width in ascending order. The fraction of ground truths above and below
their respective PIs are 2.4% and 2.69%, respectively.

In addition to the point estimation, BlaST also provides a 95% prediction
interval (PI), which equals 1.96 times the estimate standard deviation. We
have validated the coverage of the PI as shown in Figure 7. With 5.1% of the
predictions outside the blue range the 95% coverage is met. Overall, the PI
width does not vary strongly for the majority of predictions. A minor fraction
is significantly larger than the rest. In summary, the PI range has a median
total width of $1.7$ in $\log_{10}$($\nu_{\rm peak}^{S}$/Hz) space.

Both metrics are also shown as a function of $\nu_{\rm peak}^{S}$in Figure 8.
It shows that BlaST performs well on all ranges of $\nu_{\rm peak}^{S}$ and
especially for values above $10^{12.5}\,$\mathrm{H}\mathrm{z}$$ and below
$10^{17.5}\,$\mathrm{H}\mathrm{z}$$.

Figure 8: Median absolute prediction error and median prediction interval
width as a function of the $\nu_{\rm peak}^{S}$. Both curves are produced by a
moving average and smoothed using a cubic spline. Additionally, the
distribution of $\nu_{\rm peak}^{S}$in the data set is shown as black-dotted
line. We note that the median absolute error and median PI width do not depend
much on this quantity.

For application on unknown data it is crucial that BlaST recognizes non-blazar
data by assigning them with wider PIs. This property of a (machine-learning)
estimator is sometimes referred to as knowing-what-it-knows. Figure 9 shows
the distribution of PI width for three different scenarios: True blazar data
from the 4LAC-DR2 catalogue, blazar-like data simulated by adjusting Gaussian
noise to frequency and flux ranges typical for blazars, and finally the pure
Gaussian noise without any scaling. For noise that is not in the range of
typical blazar data, the algorithm predicts PI values close to 4. In case the
noise resembles the typical frequencies and fluxes of blazars, the algorithm
tends to provide smaller PI width with a long tail. Nevertheless even in this
case the PI width tends to be much larger than for real blazars. We hence
conclude that the PI width is a decent figure of merit to separate noise from
true blazar data. In the following section, and especially Figure 12 we
discuss how the PI behaves for different, $\gamma$-ray detected, source
classes.

Figure 9: Distribution of the PI widths for true blazar data, similar data
simulated by mapping normal distributed random numbers to frequency ranges and
flux levels typical for blazars, and true noise using these random numbers
without mapping.

##  4 Application to the 4LAC catalog

With an overall of 3696 objects the Fermi 4LAC-DR2 catalog (Lott et al., 2020)
121212While writing this paper the 4LAC-DR3 catalog (Abdollahi et al., 2022)
was released. We keep the comparisons to the well established DR2 here, but
provide the $\nu_{\rm peak}^{S}$for the DR3 on our github page is the largest
collection of active galaxies detected at GeV $\gamma$-ray energies. In total,
98% of the objects in the catalog are blazars. For most of those objects also
an estimation of the $\nu_{\rm peak}^{S}$ is provided. It is determined
manually by the Fermi Collaboration on a source-by-source basis, using the
spectral energy distribution available from the interactive web-tool of the
Italian Space Agency (ASI) Space Science Data Center (SSDC). Data points are
removed if they are clearly inconsistent with the overall SED or appear to be
dominated by thermal emission from the accretion disk or the host galaxy.
Subsequently, the remaining data are fitted using a third-order polynomial in
the log-log plane of the SED. Systematic uncertainties in this procedure are
mainly caused by human errors, i.e., a bad choice of the fit ranges, and the
unaccounted contamination from thermal emission (Lott et al., 2020).

In contrast to this manual approach, the BlaST tool enables an automatized
procedure. This also implies that the estimations are consistent throughout
the sample and insusceptible to human error. As the algorithm has been trained
on real and complete blazar SEDs, it is robust against thermal contributions
and faulty measurements that rarely appear in VOU-Blazars. Hence, there is no
need for manual data cleaning.

In order to run BlaST on the 4LAC-DR2 SEDs, we first download their spectral
energy distributions using the VOU blazar tool. Subsequently, the $\nu_{\rm
peak}^{S}$ is estimated on a source-to-source basis. This takes onlya a few
seconds for the entire catalog. In addition to the best-fit $\nu_{\rm
peak}^{S}$ , BlaST also provides an uncertainty estimation. This not only
accounts for missing or noisy data, but also for the non-simultaneous
observations. In the following, we study how our $\nu_{\rm peak}^{S}$
estimation compares to the ones presented in the 4LAC-DR2 paper. The full
catalog with our estimations can be accessed on
GitHub131313https://github.com/tkerscher/blast/blob/master/4LAC.csv.

Figure 10: Comparison between the estimated $\nu_{\rm peak}^{S}$of the Fermi
4LAC-DR2 catalog and the BlaST tool. The distributions are normalized along
each column, i.e., along the 4LAC-DR2 $\nu_{\rm peak}^{S}$axis.

In Figure 9 we have analyzed how the PI width behaves for artificial data and
noise. We have concluded that the PI width becomes significantly larger when
the SED does not fit to the typical blazar shape. In Figure 12 we now consider
the specific case of different source classes and $\gamma$-ray luminosities
for the sources in the 4LAC-DR2 catalog. In consistency with our previous
findings, the PI width increases for radio galaxies - as non-blazar objects -
compared to the different blazar classes. Moreover, the distributions of PI
width as a function of the $\gamma$-ray luminosity show a trend to more
certain predictions the larger the luminosity. This is expected, as the
growing $\gamma$-ray luminosity also indicates a stronger contribution of the
non-thermal part of the emission and hence less confusion with the galactic
emission. Moreover, luminous objects usually have a better data coverage
across the entire wavelength band. Note that also the pronounced peak for
FSRQs is related to the fact that FSRQs (in the definition of the 4LAC-DR2
catalog, which is disputed (Padovani et al., 2022)) tend to have a larger
$\gamma$-ray luminosity.

Figure 10 shows a comparison between the estimated $\nu_{\rm peak}^{S}$ by
BlaST and in the 4LAC-DR2 catalog. Overall, we observe a good agreement
between the two with 70% (87%) of the predictions deviating less than 0.5
(1.0) in $\log_{10}$($\nu_{\rm peak}^{S}$/Hz). Larger deviations are mainly
observed for a Fermi $\nu_{\rm peak}^{S}$ between
$10^{14}\,$\mathrm{H}\mathrm{z}$$ and $10^{15}\,$\mathrm{H}\mathrm{z}$$. This
is the region of major confusion between emission from the accretion disk and
(electron-)synchrotron processes. Figure 11 shows exemplary the top two
redshift-known SEDs with the strongest upward (downward) deviation between the
4LAC-DR2 and the BlaST $\nu_{\rm peak}^{S}$ prediction. In all those cases,
BlaST correctly accounts for emission from the host galaxy (from thermal gas
or accretion) and provides convincing $\nu_{\rm peak}^{S}$estimates and
errors. Through visual inspection of many SEDs we confirm that those are
typical problem that are cured with our novel approach. In the range of
$\nu_{\rm peak}^{S}$$\gtrsim 10^{16}$, Figure 10 also shows that BlaST tends
to estimate slightly higher peak values. Note that, as Figure 11 shows,
accounting for the emission from the galaxy and the accretion disk can shift
the estimated $\nu_{\rm peak}^{S}$ in both directions - to lower and higher
values.

Figure 11: Collection of spectral energy distributions, showing typical cases
of disagreement between BlaST and the 4LAC-DR2. The name of the respective
blazar is given in the title at the top left. Available measurements are shown
as black error bars. The $\nu_{\rm peak}^{S}$estimation of the 4LAC-DR2 is
shown in orange and the BlaST estimation with 95% prediction interval in blue.
The gray solid lines indicate the model expectations of the galactic
contributions (Mannucci et al., 2001; Vanden Berk et al., 2001). In all cases
there are clear indications that BlaST correctly recognizes galactic features
and provides the more accurate $\nu_{\rm peak}^{S}$.

Figure 12: Distribution of the PI width for different source classes (top; RDG
= Radio Galaxy, BCU=Blazar of uncertain type) and $\gamma$-ray luminosities
(bottom; in erg/s ).

In the following, we reproduce the $\nu_{\rm peak}^{S}$-dependent source count
distributions from the 4LAC-DR2 paper and discuss how they change under the
BlaST estimation. First, Figure 13 shows the distributions of synchrotron
peaks for all blazars which have an estimation in the 4LAC-DR2 catalog. In
addition, we also show a number density estimation based on the best-fit value
$\hat{\mu}$ and the standard deviation $\sigma$ predicted by BlaST , i.e.,

$N(x)=\Delta
x\cdot\sum\frac{1}{\sqrt{2\pi}\hat{\sigma}}\exp\left(-\frac{1}{2}\left(\frac{x-\hat{\mu}}{\sigma}\right)^{2}\right)$
(3)

where $\Delta x$ represents the bin width of the histogram. Most prominently,
the new estimation produces a stronger peak of sources for $\nu_{\rm
peak}^{S}$$\lesssim 10^{13}\,$\mathrm{H}\mathrm{z}$$. While further studies
are needed, this might be a hint in the direction that blazars can be
separated into two distinct classes - LBLs with $\nu_{\rm peak}^{S}$
$<10^{13.5}$ Hz and IHBLs with $\nu_{\rm peak}^{S}$ $\geq 10^{13.5}$ Hz
(Giommi and Padovani, 2021). The peak is a bit wider if we consider the kernel
density estimation. This is expected as the width of the $\nu_{\rm peak}^{S}$
prediction intervals in Fig. 8 behaves as a smearing parameter. In the future,
this smearing might be reduced with additional data points in the SED and
better instruments.

Finally, in Figure 14 we study the connection between the power-law spectral
index of the $\gamma$-ray emission and the estimated $\nu_{\rm peak}^{S}$. We
find that, overall, the their relation is comparable for the 4LAC-DR2 catalog
and BlaST, with slightly fewer outliers in the latter case.

Figure 13: Distribution of the synchrotron peak for the identified blazars in
the 4LAC-DR2 catalog. Only objects that have a defined $\nu_{\rm peak}^{S}$in
the 4LAC-DR2 are shown. The distribution based on the BlaST estimations are
shown in blue. In addition to a histogram also a kernel density estimation is
shown based on the best-fit value and uncertainty estimation of each sources.
For comparison the 4LAC-DR2 distribution is shown in orange. Figure 14: The
distribution of the $\gamma$-ray photon index as a function of the estimated
$\nu_{\rm peak}^{S}$in the 4LAC-DR2 catalog (left) and using the BlaST tool
(right). The solid lines indicate the median values.

##  5 Conclusion

The synchrotron emission peak ($\nu_{\rm peak}^{S}$) is an important measure
for the classification and characterization of blazar jets. In this paper we
have presented a new tool - BlaST - that utilizes machine learning methods to
estimate the $\nu_{\rm peak}^{S}$ directly from the binned spectral energy
distribution without any need for manual data preparation. By construction,
the tool accounts for emission from the dust in the Galaxy (infrared bump) and
the accretion disk (blue bump). In contrast to other methods, it also provides
a reliable uncertainty estimation. The training has been done on a sample
containing the blazars from the BZCAT-5th edition, 4LAC-DR2, and the 3HSP
catalog using the measurements available through the Open Universe VOU-Blazars
tool. The final predictions of the tool are nearly instantaneous.

In order to validate and check the BlaST estimation and its improvements, we
have compared them against the $\nu_{\rm peak}^{S}$ given in the 4LAC-DR2
catalog. While the majority of predictions are well compatible with the 4LAC-
DR2, we find that BlaST is more robust against confusions with the blue- and
infrared bump. A refined estimation of the $\nu_{\rm peak}^{S}$ of the entire
4LAC-DR2 catalog can be found on the project GitHub
page141414https://github.com/tkerscher/blast. In addition, we consider BlaST
particularly useful for the extension of existing blazar catalogs (as for
example in Giommi et al. (2020)), fast estimations of the $\nu_{\rm peak}^{S}$
in real-time astronomy pipelines, and the robust estimation of $\nu_{\rm
peak}^{S}$ for non-experts.

##  6 Acknowledgements

We thank Elisa Resconi and Martina Karl for proofreading the paper and
acknowledge the use of data and software facilities from the SSDC, managed by
the Italian Space Agency, and the United Nations “Open Universe” initiative.
We also thank the unknown reviewers for reading the paper thoroughly and
providing valuable comments. PG and TG acknowledge the support of the
Technische Universität München - Institute for Advanced Study, funded by the
German Excellence Initiative (and the European Union Seventh Framework
Programme under grant agreement n. 291763).

## References

  * Aartsen et al. (2018a) Aartsen, M.G., et al. (IceCube, Fermi-LAT, MAGIC, AGILE, ASAS-SN, HAWC, H.E.S.S., INTEGRAL, Kanata, Kiso, Kapteyn, Liverpool Telescope, Subaru, Swift NuSTAR, VERITAS, VLA/17B-403), 2018a.  Multimessenger observations of a flaring blazar coincident with high-energy neutrino IceCube-170922A.  Science 361, eaat1378.  doi:10.1126/science.aat1378, arXiv:1807.08816.
  * Aartsen et al. (2018b) Aartsen, M.G., et al. (IceCube), 2018b.  Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert.  Science 361, 147–151.  doi:10.1126/science.aat2890, arXiv:1807.08794.
  * Abdo et al. (2010) Abdo, A.A., Ackermann, M., Agudo, I., Ajello, M., Aller, H.D., Aller, M.F., Angelakis, E., , et al., 2010\.  The Spectral Energy Distribution of Fermi Bright Blazars.  ApJ 716, 30–70.  doi:10.1088/0004-637X/716/1/30, arXiv:0912.2040.
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., Ajello, M., Atwood, W.B., Axelsson, M., Baldini, L., et al., 2020\.  Fermi Large Area Telescope Fourth Source Catalog.  ApJS 247, 33.  doi:10.3847/1538-4365/ab6bcb, arXiv:1902.10045.
  * Abdollahi et al. (2022) Abdollahi, S., et al. (Fermi-LAT), 2022\.  Incremental Fermi Large Area Telescope Fourth Source Catalog arXiv:2201.11184.
  * Ackermann et al. (2015a) Ackermann, M., Ajello, M., Atwood, W.B., Baldini, L., Ballet, J., Barbiellini, G., Bastieri, D., et al., 2015a.  The Third Catalog of Active Galactic Nuclei Detected by the Fermi Large Area Telescope.  ApJ 810, 14.  doi:10.1088/0004-637X/810/1/14, arXiv:1501.06054.
  * Ackermann et al. (2015b) Ackermann, M., Ajello, M., Atwood, W.B., Baldini, L., Ballet, J., Barbiellini, G., Bastieri, D., et al., 2015b.  VizieR Online Data Catalog: Third catalog of LAT-detected AGNs (3LAC) (Ackermann+, 2015).  VizieR Online Data Catalog , J/ApJ/810/14.
  * Ballet et al. (2020) Ballet, J., Burnett, T.H., Digel, S.W., Lott, B., 2020\.  Fermi Large Area Telescope Fourth Source Catalog Data Release 2.  arXiv e-prints , arXiv:2005.11208arXiv:2005.11208.
  * Candes et al. (2006) Candes, E., Romberg, J., Tao, T., 2006.  Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information.  IEEE Transactions on Information Theory 52, 489–509.  doi:10.1109/TIT.2005.862083.
  * Cerruti et al. (2017) Cerruti, M., Benbow, W., Chen, X., Dumm, J.P., Fortson, L.F., Shahinyan, K., 2017\.  Luminous and high-frequency peaked blazars: the origin of the $\gamma$-ray emission from PKS 1424+240.  A&A 606, A68.  doi:10.1051/0004-6361/201730799, arXiv:1707.00804.
  * Chang et al. (2019a) Chang, Y.L., Arsioli, B., Giommi, P., Padovani, P., Brandt, C., 2019a.  VizieR Online Data Catalog: 3HSP blazars catalog (Chang+, 2019).  VizieR Online Data Catalog , J/A+A/632/A77.
  * Chang et al. (2019b) Chang, Y.L., Arsioli, B., Giommi, P., Padovani, P., Brandt, C.H., 2019b.  The 3HSP catalogue of extreme and high-synchrotron peaked blazars.  A&A 632, A77.  doi:10.1051/0004-6361/201834526, arXiv:1909.08279.
  * Chang et al. (2019c) Chang, Y.L., Arsioli, B., Giommi, P., Padovani, P., Brandt, C.H., 2019c.  The 3HSP catalogue of extreme and high-synchrotron peaked blazars.  A&A 632, A77.  doi:10.1051/0004-6361/201834526, arXiv:1909.08279.
  * Chang et al. (2020) Chang, Y.L., Brandt, C.H., Giommi, P., 2020.  The Open Universe VOU-Blazars tool.  Astronomy and Computing 30, 100350\.  doi:10.1016/j.ascom.2019.100350, arXiv:1909.11455.
  * Evans et al. (2020) Evans, P.A., Page, K.L., Osborne, J.P., Beardmore, A.P., Willingale, R., Burrows, D.N., Kennea, J.A., Perri, M., Capalbi, M., Tagliaferri, G., Cenko, S.B., 2020.  2SXPS: An Improved and Expanded Swift X-Ray Telescope Point-source Catalog.  ApJS 247, 54.  doi:10.3847/1538-4365/ab7db9, arXiv:1911.11710.
  * Falomo et al. (2014) Falomo, R., Pian, E., Treves, A., 2014.  An optical view of BL Lacertae objects.  A&A Rev. 22, 73.  doi:10.1007/s00159-014-0073-z, arXiv:1407.7615.
  * Fossati et al. (1998) Fossati, G., Maraschi, L., Celotti, A., Comastri, A., Ghisellini, G., 1998.  A unifying view of the spectral energy distributions of blazars.  MNRAS 299, 433–448.  doi:10.1046/j.1365-8711.1998.01828.x, arXiv:astro-ph/9804103.
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., Mason, K.O., Nousek, J.A., Wells, A.A., White, N.E., et al., 2004\.  The Swift Gamma-Ray Burst Mission.  ApJ 611, 1005–1020.  doi:10.1086/422091.
  * Ghisellini et al. (1985) Ghisellini, G., Maraschi, L., Treves, A., 1985.  Inhomogeneous synchrotron-self-compton models and the problem of relativistic beaming of BL Lac objects.  A&A 146, 204–212.
  * Ghisellini and Tavecchio (2008) Ghisellini, G., Tavecchio, F., 2008\.  The blazar sequence: a new perspective.  MNRAS 387, 1669–1680.  doi:10.1111/j.1365-2966.2008.13360.x, arXiv:0802.1918.
  * Giommi et al. (2019) Giommi, P., Brandt, C.H., Barres de Almeida, U., Pollock, A.M.T., Arneodo, F., Chang, Y.L., Civitarese, O., De Angelis, M., D’Elia, V., Del Rio Vera, J., Di Pippo, S., Middei, R., Penacchioni, A.V., Perri, M., Ruffini, R., Sahakyan, N., Turriziani, S., 2019.  Open Universe for Blazars: a new generation of astronomical products based on 14 years of Swift-XRT data.  A&A 631, A116.  doi:10.1051/0004-6361/201935646, arXiv:1904.06043.
  * Giommi et al. (2020) Giommi, P., Glauch, T., Padovani, P., Resconi, E., Turcati, A., Chang, Y.L., 2020\.  Dissecting the regions around IceCube high-energy neutrinos: growing evidence for the blazar connection.  MNRAS 497, 865–878.  doi:10.1093/mnras/staa2082, arXiv:2001.09355.
  * Giommi and Padovani (2021) Giommi, P., Padovani, P., 2021\.  Astrophysical Neutrinos and Blazars arXiv:2112.06232.
  * Giommi et al. (2012a) Giommi, P., Padovani, P., Polenta, G., Turriziani, S., D’Elia, V., Piranomonte, S., 2012a.  A simplified view of blazars: clearing the fog around long-standing selection effects.  MNRAS 420, 2899–2911.  doi:10.1111/j.1365-2966.2011.20044.x, arXiv:1110.4706.
  * Giommi et al. (2021) Giommi, P., Perri, M., Capalbi, M., D’Elia, V., Barres de Almeida, U., Brandt, C.H., Pollock, A.M.T., et al., 2021\.  X-ray spectra, light curves and SEDs of blazars frequently observed by Swift.  MNRAS 507, 5690–5702.  doi:10.1093/mnras/stab2425, arXiv:2108.07255.
  * Giommi et al. (2012b) Giommi, P., Polenta, G., Lähteenmäki, A., Thompson, D.J., Capalbi, M., Cutini, S., Gasparrini, D., et al., 2012b.  Simultaneous Planck, Swift, and Fermi observations of X-ray and $\gamma$-ray selected blazars.  A&A 541, A160.  doi:10.1051/0004-6361/201117825, arXiv:1108.1114.
  * Giommi et al. (2020) Giommi, P., et al., 2020.  Open Universe survey of Swift-XRT GRB fields: Flux-limited sample of HBL blazars.  Astron. Astrophys. 642, A141.  doi:10.1051/0004-6361/202037921, arXiv:2003.05153.
  * Gneiting and Raftery (2007) Gneiting, T., Raftery, A.E., 2007\.  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association 102, 359–378.  URL: https://doi.org/10.1198/016214506000001437, doi:10.1198/016214506000001437, arXiv:https://doi.org/10.1198/016214506000001437.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., Courville, A., 2016.  Deep Learning.  MIT Press.  http://www.deeplearningbook.org.
  * Gordon et al. (2020) Gordon, Y.A., Boyce, M.M., O’Dea, C.P., Rudnick, L., Andernach, H., Vantyghem, A.N., Baum, S.A., Bui, J.P., Dionyssiou, M., 2020.  A Catalog of Very Large Array Sky Survey Epoch 1 Quick Look Components, Sources, and Host Identifications.  Research Notes of the American Astronomical Society 4, 175.  doi:10.3847/2515-5172/abbe23.
  * Homan et al. (2021) Homan, D.C., Cohen, M.H., Hovatta, T., Kellermann, K.I., Kovalev, Y.Y., Lister, M.L., Popkov, A.V., Pushkarev, A.B., Ros, E., Savolainen, T., 2021\.  MOJAVE XIX: Brightness Temperatures and Intrinsic Properties of Blazar Jets.  arXiv e-prints , arXiv:2109.04977arXiv:2109.04977.
  * Keenan et al. (2021) Keenan, M., Meyer, E.T., Georganopoulos, M., Reddy, K., French, O.J., 2021.  The relativistic jet dichotomy and the end of the blazar sequence.  MNRAS 505, 4726–4745.  doi:10.1093/mnras/stab1182, arXiv:2007.12661.
  * Lakshminarayanan et al. (2017) Lakshminarayanan, B., Pritzel, A., Blundell, C., 2017.  Simple and scalable predictive uncertainty estimation using deep ensembles, in: Proceedings of the 31st International Conference on Neural Information Processing Systems, Curran Associates Inc., Red Hook, NY, USA. p. 6405–6416.
  * Liaw et al. (2018) Liaw, R., Liang, E., Nishihara, R., Moritz, P., Gonzalez, J.E., Stoica, I., 2018\.  Tune: A research platform for distributed model selection and training.  arXiv preprint arXiv:1807.05118 .
  * Lott et al. (2020) Lott, B., Gasparrini, D., Ciprini, S., 2020.  The Fourth Catalog of Active Galactic Nuclei Detected by the Fermi Large Area Telescope – Data Release 2.  arXiv e-prints , arXiv:2010.08406arXiv:2010.08406.
  * Mairal et al. (2009) Mairal, J., Bach, F., Ponce, J., Sapiro, G., 2009\.  Online dictionary learning for sparse coding.  Proceedings of the 26th Annual International Conference on Machine Learning - ICML 09 doi:10.1145/1553374.1553463.
  * Mannucci et al. (2001) Mannucci, F., Basile, F., Poggianti, B.M., Cimatti, A., Daddi, E., Pozzetti, L., Vanzi, L., 2001.  Near-infrared template spectra of normal galaxies: k-corrections, galaxy models and stellar populations.  Mon. Not. Roy. Astron. Soc. 326, 745\.  doi:10.1046/j.1365-8711.2001.04628.x, arXiv:astro-ph/0104427.
  * Massaro et al. (2015) Massaro, E., Maselli, A., Leto, C., Marchegiani, P., Perri, M., Giommi, P., Piranomonte, S., 2015.  The 5th edition of the Roma-BZCAT. A short presentation.  Ap&SS 357, 75.  doi:10.1007/s10509-015-2254-2, arXiv:1502.07755.
  * Massaro et al. (2016) Massaro, E., Maselli, A., Leto, C., Marchegiani, P., Perri, M., Giommi, P., Piranomonte, S., 2016.  VizieR Online Data Catalog: The Roma BZCAT - 5th edition (Massaro+, 2015).  VizieR Online Data Catalog , VII/274.
  * Meinshausen (2006) Meinshausen, N., 2006.  Quantile regression forests.  Journal of Machine Learning Research 7, 983–999.  URL: http://jmlr.org/papers/v7/meinshausen06a.html.
  * Padovani et al. (2017) Padovani, P., Alexander, D.M., Assef, R.J., De Marco, B., Giommi, P., Hickox, R.C., Richards, G.T., Smolčić, V., Hatziminaoglou, E., Mainieri, V., Salvato, M., 2017.  Active galactic nuclei: what’s in a name?  A&A Rev. 25, 2.  doi:10.1007/s00159-017-0102-9, arXiv:1707.07134.
  * Padovani and Giommi (1995) Padovani, P., Giommi, P., 1995\.  The connection between x-ray- and radio-selected BL Lacertae objects.  ApJ 444, 567–581.
  * Padovani et al. (2022) Padovani, P., Giommi, P., Falomo, R., Oikonomou, F., Petropoulou, M., Glauch, T., Resconi, E., Treves, A., Paiano, S., 2022.  The spectra of IceCube neutrino (SIN) candidate sources – II. Source characterization.  Mon. Not. Roy. Astron. Soc. 510, 2671–2688.  doi:10.1093/mnras/stab3630, arXiv:2112.05394.
  * Padovani et al. (2018) Padovani, P., Giommi, P., Resconi, E., Glauch, T., Arsioli, B., Sahakyan, N., Huber, M., 2018.  Dissecting the region around IceCube-170922A: the blazar TXS 0506+056 as the first cosmic neutrino source.  Mon. Not. Roy. Astron. Soc. 480, 192–203.  doi:10.1093/mnras/sty1852, arXiv:1807.04461.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chintala, S., 2019.  Pytorch: An imperative style, high-performance deep learning library, in: Advances in Neural Information Processing Systems 32. Curran Associates, Inc., pp. 8024–8035.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E., 2011\.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research 12, 2825–2830.
  * Prandini and Ghisellini (2022) Prandini, E., Ghisellini, G., 2022\.  The Blazar Sequence and Its Physical Understanding.  Galaxies 10, 35\.  doi:10.3390/galaxies10010035, arXiv:2202.07490.
  * Stathopoulos et al. (2021) Stathopoulos, S.I., Petropoulou, M., Giommi, P., Vasilopoulos, G., Padovani, P., Mastichiadis, A., 2021\.  Probing Neutrino Emission from X-ray Blazar Flares observed with Swift-XRT.  arXiv e-prints , arXiv:2107.14632arXiv:2107.14632.
  * Tibshirani (1996) Tibshirani, R., 1996.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological) 58, 267–288.  URL: http://www.jstor.org/stable/2346178.
  * Vanden Berk et al. (2001) Vanden Berk, D.E., et al. (SDSS), 2001\.  Composite quasar spectra from the Sloan Digital Sky Survey.  Astron. J. 122, 549–564.  doi:10.1086/321167, arXiv:astro-ph/0105231.
  * Webb et al. (2020) Webb, N.A., Coriat, M., Traulsen, I., Ballet, J., Motch, C., Carrera, F.J., Koliopanos, F., Authier, J., de la Calle, I., Ceballos, M.T., Colomo, E., Chuard, D., Freyberg, M., Garcia, T., Kolehmainen, M., Lamer, G., Lin, D., Maggi, P., Michel, L., Page, C.G., Page, M.J., Perea-Calderon, J.V., Pineau, F.X., Rodriguez, P., Rosen, S.R., Santos Lleo, M., Saxton, R.D., Schwope, A., Tomás, L., Watson, M.G., Zakardjian, A., 2020.  The XMM-Newton serendipitous survey. IX. The fourth XMM-Newton serendipitous source catalogue.  A&A 641, A136.  doi:10.1051/0004-6361/201937353, arXiv:2007.02899.
  * Zabalza (2015) Zabalza, V., 2015.  naima: a python package for inference of relativistic particle energy distributions from observed nonthermal spectra.  Proc. of International Cosmic Ray Conference 2015 , 922arXiv:1509.03319.