11institutetext: Department of Physics, National and Kapodistrian University
of Athens, University Campus Zografos, GR 15784, Athens, Greece
22institutetext: Institute of Accelerating Systems & Applications, University
Campus Zografos, Athens, Greece

# Application of neural networks to synchro-Compton blazar emission models

A. Tzavellas 11<EMAIL_ADDRESS>G. Vasilopoulos 1122
<EMAIL_ADDRESS>M. Petropoulou 1122 A. Mastichiadis 11 S. I. Stathopoulos
1122

(Received … ; accepted …)

###### Abstract

Context. Jets from supermassive black holes in the centers of active galaxies
are the most powerful persistent sources of electromagnetic radiation in the
Universe. To infer the physical conditions in the otherwise out-of-reach
regions of extragalactic jets we usually rely on fitting of their spectral
energy distribution (SED). The calculation of radiative models for the jet
non-thermal emission usually relies on numerical solvers of coupled partial
differential equations.

Aims. In this work machine learning is used to tackle the problem of high
computational complexity in order to significantly reduce the SED model
evaluation time, which is needed for SED fitting with Bayesian inference
methods.

Methods. We compute SEDs based on the synchrotron self-Compton model for
blazar emission using the radiation code ATHE$\nu$A, and use them to train
Neural Networks exploring whether these can replace the original computational
expensive code.

Results. We find that a Neural Network with Gated Recurrent Unit neurons can
effectively replace the ATHE$\nu$A leptonic code for this application, while
it can be efficiently coupled with MCMC and nested sampling algorithms for
fitting purposes. We demonstrate this through an application to simulated data
sets and with an application to observational data. We offer this tool in the
community through a public repository.

Conclusions. We present a proof-of-concept application of neural networks to
blazar science. This is the first step in a list of future applications
involving hadronic processes and even larger parameter spaces.

###### Key Words.:

Methods: numerical, statistical; Radiation mechanisms: non-thermal, Radiative
transfer

##  1 Introduction

Active galactic nuclei (AGN) are known to launch relativistic and collimated
plasma outflows, known as jets. While jetted AGN consist a small fraction of
the AGN population (Padovani et al., 2017), they have attracted a lot of
attention over the past decades because they are the most powerful and
persistent sources of non-thermal electromagnetic radiation (for recent
reviews, see Blandford et al., 2019; Hovatta & Lindfors, 2019). AGN with jets
that are directed toward us are known as blazars. The combination of the small
viewing angle with the relativistic motion of the emitting plasma results in
strong Doppler beaming, making blazars ideal sources for studying non-thermal
radiation processes in jets.

Central to the study of blazars is their spectral energy distribution (SED),
which can be thought of as a metric that maps the energy output of the jet
across different photon energies, presenting a panoramic view of the radiative
processes within the jet. Decoding the SED is pivotal for uncovering the
prevailing physical conditions and radiative processes at play in the
otherwise out-of-reach regions of extragalactic jets. A significant hurdle in
this process, however, is the fitting of the SED itself, which relies on a
juxtaposition of the observed energy distributions with model predictions.
Radiative models for jet emission are primarily numerical, as they require the
solution of a system of stiff coupled partial differential equations (PDEs)
describing the evolution with time of the distribution functions of radiating
particles and photons (e.g. Dimitrakoudis, S. et al., 2012; Gao et al., 2017;
Cerruti et al., 2015; Gasparyan et al., 2022). Such numerical codes tend to
have a high computational complexity, especially in the case of lepto-hadronic
models111Models that account for the radiation processes of relativistic
electrons and protons are known as (lepto-)hadronic., where each SED
computation could last from a few minutes up to an hour (depending on the
numerical schemes used by the PDE solver, but also on the approximations made
for each operator in the PDE). If such numerical model were to be used in a
Markov Chain Monte Carlo (MCMC) algorithm to determine the posterior
distributions of the model parameters, the overall execution time could become
prohibitively long. Addressing this computational bottleneck is central to
accelerating the pace of discovery in blazar science.

Although the idea behind a Neural Network (NN) has been around for more than
50 years (McCulloch & Pitts, 1943), it’s extended application has been
feasible with hardware advantages over the last decade. Inspired by recent
trends in computational biology Wang et al. (2019) and similar applications in
other disciplines (e.g. Wiecha & Muskens, 2020; Alarfaj et al., 2022), this
work investigates the possibility of replacing a numerical model of average
computational complexity, i.e. the synchrotron self-Compton222This is a
leptonic model where low-energy synchrotron photons produced by relativistic
electrons are up-scattered to high-energies via inverse Compton off the same
electron population.  (SSC, e.g. Maraschi et al., 1992; Bloom & Marscher,
1996; Mastichiadis & Kirk, 1997) with a NN (see Abiodun et al., 2018). In
pursuit of this objective, the initial phase involves the meticulous assembly
of a dataset derived from the numerical model. While the intricate processes
of dataset generation are set aside for brevity, it is imperative that this
dataset is comprehensive and uniformly spans the parameter space. Furthermore,
this compilation should be both ample in size to enable rigorous training and
analysis, and efficiently constructed to conserve time and resources. With a
robust dataset in hand, the exploration of various NN architectures takes
center stage.

The realm of machine learning (ML) is characterized by its blend of scientific
rigor and artistic intuition; there exists no monolithic model that stands as
the unequivocal solution. Each model is a unique amalgamation of structural
nuances and algorithmic intricacies, tailored to address specific aspects of
complex problems. Consequently, we will investigate a suite of NN
configurations to identify those offering optimal performance and accuracy in
the context of the defined problem space.

This paper is structured as follows. In Sec. 2 we outline the methodology used
to construct a model based on a NN, which includes the use of a numerical code
to construct a dataset and the various NN tuning and evaluation tests. In Sec.
3 we present a few scenarios where the NN trained model is coupled with a
Bayesian interface to fit simulated data and actual blazar observations. We
also provide notes on the efficiency of our approach compared to numerical
codes. Finally, we present the conclusions of this work and some future
aspects in Sec. 4.

Figure 1: Implementation flow chart describing our methodology.

##  2 Methods

Our approach, which is schematically shown in Fig. 1, is comprised of the
following steps: (i) we built a database of blazar SED models, (ii) we train
the NN using available datasets, and (iii) we evaluate our results based on
datasets not used in the training. In the following paragraphs we elaborate on
the above steps.

###  2.1 Generation of dataset

To generate blazar SEDs we use the radiative transfer code ATHE$\nu$A
(Mastichiadis & Kirk, 1995; Dimitrakoudis et al., 2012) that solves a system
of coupled PDEs describing the evolution of relativistic particle
distributions and photons contained in a non-thermal emitting source. For this
project we use the leptonic module of the code that describes a magnetized
source containing relativistic electron-positron pairs and photons. Electrons
are being injected into the source at a fixed rate (or injection luminosity),
and escape on an energy-independent timescale that is equal to the light-
crossing time of the source. Photons are being produced via electron
synchrotron radiation and scattered to higher energies through inverse Compton
scattering, while they escape on the light-crossing time from the
source333Additional processes, like photon-photon pair production and
synchrotron self-absorption, that may affect the lowest and highest energy
parts of the photon spectrum are not included in this proof-of-concept work.
Trained NNs with these procedures will be gradually added to the GitHub
repository of the project.. For a detailed description of the modeling of each
physical process, we refer the reader to Mastichiadis & Kirk (1995);
Dimitrakoudis et al. (2012)).

The SSC model parameters, which are used as an input to ATHE$\nu$A, are:

  * •

the characteristic size of the emitting region $R$ (i.e., its radius assuming
a spherical source in the jet rest frame),

  * •

the magnetic field strength $B$ (measured in the jet rest frame),

  * •

the electron injection compactness $\ell_{e}$, which is a dimensionless
measure of the power injected into relativistic electrons (defined as
$\ell_{e}=\sigma_{T}L_{e}/(4\pi Rm_{e}c^{3})$),

  * •

the power-law slope $p$ of the electron distribution at injection (i.e.,
$dN_{e}/dtd\gamma_{e}\propto\gamma_{e}^{-p}$), and

  * •

the Lorentz factor range of radiating electrons,
$[\gamma_{min},\gamma_{max}]$, also at injection.

For the generation of the dataset we use parameter values that are motivated
by SSC models of blazar emission (e.g. Celotti & Ghisellini, 2008). For
instance, we take the lower bound of the electron distribution to range
between $10^{0.1}$ and $10^{4}$. To avoid very narrow electron distributions,
or cases where $\gamma_{max}<\gamma_{min}$ that will lead to code crashing, we
set the upper cutoff of the distribution to lie between $10^{2}\gamma_{min}$
and $\min(10^{8},\gamma_{H})$; here $\gamma_{H}=eBR/(m_{e}c^{2})$ is an upper
limit imposed by the Hillas condition (Hillas, 1984), which ensures magnetic
confinement of the most energetic electrons in the source. Each variable is
randomly drawn from a uniform distribution in logarithmic space (expect for
the slope $p$) – see Table 1. Given the theoretically known dependence of the
synchrotron/SSC fluxes and the characteristic photon energies on model
parameters, like $R$ and $B$ (see e.g. Mastichiadis & Kirk, 1997; Finke et
al., 2008), the selected ranges are wide enough to produce a diverse dataset,
which is crucial for the NN training (see next subsections).

Another parameter that should be taken into account for the SED fitting is the
Doppler factor $\delta$ of the emitting region, which is used to transform the
photon spectra from the jet comoving frame to the observer’s frame. However,
this is a nuisance parameter for the dataset creation, since the Doppler
boosting can be applied at a later stage to transform the comoving photon
spectra (unprimed quantities) to the observer’s frame, i.e. $\nu_{\rm
obs}=\delta\nu$ and $L_{\nu,\rm obs}=\delta^{3}L_{\nu}$ (here, $\nu$ is the
photon frequency and $L_{\nu}$ is the specific luminosity).

To generate the dataset we sample the space of leptonic parameters and feed
each set to ATHE$\nu$A. The system evolves for $5\leavevmode\nobreak\ R/c$,
which is sufficiently long for the establishment of a steady state
(equilibrium solution of the PDEs). The parameters along with the photon PDE
solution are concatenated to form the entry of each dataset. All SEDs were
computed at the same 500 grid points of photon energies, logarithmically
spaced between $10^{-15}$ and $10^{10}$ (in units of $m_{\rm e}c^{2}$). The
SED fluxes returned by ATHE$\nu$A at each grid point ($\nu F_{\nu}$ in code
units) were re-normalized using Min-Max scaling between the values of 0 and 1,
before being used for the training. Normalization of inputs is standard
practice (Chris Bishop et al., 1995). In total, we created a dataset of 10,000
samples, for the range of parameters listed in Table 1, within $\sim$500 CPU
hours (i.e. about 3 min per model) on a desktop computer (AMD Ryzen 5950X).
The generated dataset is split into 80% for training, 10% for validation and
10% for testing.

Table 1: Input parameter ranges used for the generation of training and validation datasets with the ATHE$\nu$A code. $U()$ stands for a uniform distribution. Parameter [unit] | Range
---|---
$\log R$ [cm] | $U(14,17)$
$\log B$ [G] | $U(-2,2)$
$\log\gamma_{min}$ | $U(0.1,4)$
$\log\gamma_{max}$ | $U(\log 10^{2}\gamma_{min},\min\left(8,\log\gamma_{H}\right))$
$\log\ell_{e}$ | $U(-5,-1)$
$p$ | $U(1.5,3)$

###  2.2 Training of the NNs

Once a dataset has been obtained, it is essential to select the appropriate
type of NN neuron and its corresponding structure. The prediction of the SED
can be perceived as a regression problem, making the simple artificial neuron
an apparent choice. However, it is also possible to approach the bins of the
SED as sequential data. In this case, all three types of recurrent neurons
become potential options for constructing the NN. We employed several NN
topologies that varied with respect to the following parameters: number of
hidden layers, number of neurons per hidden layer, and three types of neurons:
an artificial NN with a deep stack of hidden layers, i.e. Artificial Neural
Network (ANN), a Gated Recurrent Unit (GRU, Cho et al., 2014) and a Long
Short-Term Memory (LSTM, Hochreiter & Schmidhuber, 1997).

Currently, there are a few ML frameworks available for several programming
languages that can be used for practical applications. Python programming
language is a very popular choice in the astrophysics community, and as a
result there are a lot of scientific software packages available for use. The
major ML packages available in Python are: Tensorflow created by Google Brain
(Abadi et al., 2015), PyTorch developed by Facebook’s AI Research lab (Paszke
et al., 2019), scikit-learn built by David Cournapeau (Pedregosa et al.,
2011), Theano developed by the Montreal Institute for Learning Algorithms (The
Theano Development Team et al., 2016). Tensorflow and Theano use Keras
(Chollet et al., 2015), an open-source software library that provides a Python
interface for NNs. It was developed with a focus on enabling fast
experimentation and prototyping, being both user-friendly and modular. In this
work we opted for Tensorflow.

We trained the 3 NNs using the 10k dataset; the hyper-parameters and
complexity of each NN are presented in Table 2. An additional necessary
parameter is the number of epochs, which refers to the number of training
cycles. The appropriate number of epochs must be carefully selected for each
network to address the issue of over-fitting. To determine the optimal number
of epochs for each network, an initial training phase is conducted with a
substantially large number of epochs (set at 3,000), since experience
demonstrated that over-fitting, if it occurs, typically starts after 200
epochs. Finally, one can choose a loss function other than the typical mean-
squared error (MSE). Since our networks handle normalized data it was decided
to use instead the mean squared logarithmic error (MSLE) as a loss function.
The training process of these networks is observed in real-time using
Tensorboard. This enables the monitoring of whether a network is over-fitting
or if it has achieved a satisfactory loss error, which is
$\mathcal{O}(10^{-4})$. Whichever of these conditions is met first dictates
the stopping point. The epoch at which training is stopped is noted, and then
the network undergoes retraining for this specific number of epochs.

Table 2: Hyper-parameters for the NNs.

Parameter Cell Type | ANN | GRU | LSTM
---|---|---|---
Hidden Layers | 6 | 6 | 6
Neurons per Layer | 184 | 184 | 184
Learning Rate ($\times 10^{-4}$) | $1.7$ | $1.7$ | $1.7$
Batch Size | 32 | 32 | 32
NN Complexity* | 300k | 1.4M | 1.9M

  * *

Total number of model parameters in thousands (k) or millions (M).

###  2.3 Evaluation of results

Table 3: Evaluation of results.

| Cell Type
---|---
|  ANN |  GRU |  LSTM
Metric | mean score*
R2 |  0.902 |  0.985 |  0.989
MSE |  1.84 |  0.23 |  0.20
MAE |  0.89 |  0.13 |  0.13
KS |  0.256 |  0.070 |  0.068
DTW |  8.86 |  3.12 |  4.37
| Ranking percentages (%)**
MSE |  0/1.3/98 |  42/56/1 |  57/42/0
MAE |  0/0 /100 |  79/ 21/0 |  21/79/0
KS |  0.5/2.9/97 |  54/44 /1.8 |  46/53/1.6
DTW |  0.00/1.2/99 |  87/12/0.7 |  13/86/0.4

  * *

Median metric value for each NN based on the test sample.

  * **

Percentage of times where each NN ranked 1st/2nd/3rd in their between
comparison for the test sample.

(a)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[hc,up][vc,hc]
(b)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(c)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]

Figure 2: Visual comparison of GRU (magenta points) and LSTM (maroon points)
NNs on top of the ATHE$\nu$A SED. Both NN generated models capture the overall
trends, however GRU offers a better description of the ATHE$\nu$A SED with
less scatter.

(a)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(b)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(c)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(d)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(e)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(f)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(g)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(h)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(i)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(j)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(h)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(k)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(l)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(m)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(n)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(o)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(p)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(q)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(r)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(s)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]
(t)left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]

Figure 3: Examples of SEDs from the test dataset (i.e., not used for the
training). Predictions of the trained GRU NN are overplotted (colored markers)
to the numerical models created with ATHE$\nu$A (solid lines).

After training our 3 NNs we proceed to evaluation of the models by comparing
the NN predictions to the SEDs from the test sample. We first use standard
metrics, namely the R2 score, the MSE, the mean-squared absolute error (MAE),
the Kolmogorov-Smirnov (KS) criterion, and Dynamic Time Warping (DTW), to
compare the performance of the trained NNs. In terms of R2 the GRU and LSTM
NNs have significantly better values than the ANN network. The ranking
results, which are listed in Table 3, show that the GRU is ranked first among
the three by three metrics, while the ANN is ranked third by all criteria. We
also visually inspect the test dataset, where we compare the prediction of the
NNs against the result of ATHE$\nu$A– see Fig. 3. From this qualitative test,
it became clear that ANN failed to describe the SED, suggesting that more
complexity has to be added to the network. Meanwhile, both GRU and LSTM
delivered acceptable solutions (based on the eyeball test). The main
difference was that GRU offered more smooth solutions for the SED, while the
LSTM predictions had noticeably small amplitude scatter across neighboring
values. In combination with the evaluation results listed in Table 3, we
select the GRU NN for the applications that follow.

In Fig. 3 we present random examples drawn from the test dataset that
highlight the wide variety of SEDs (in terms of flux and shape) reproduced by
the GRU NN, including spectra that are not typical of blazars. For instance,
there are cases without a clear two-hump morphology due to the smooth
superposition of the synchrotron and SSC components (see e.g. panels a and b),
or due to a very extended synchrotron spectrum (panel q), which reminds more
of spectra of pulsar wind nebulae (see e.g., Amato & Olmi, 2021, for Crab
nebula). Moreover, the sample contains spectra produced in the so-called
inverse Compton catastrophe limit (Readhead, 1994; Petropoulou et al., 2015b);
these SEDs consist of three distinct components attributed to synchrotron,
first-order and second-order SSC emission, with the latter carrying most of
the bolometric power of the source (see panels l and o). In most cases, the
trained NN describes accurately the numerical model (e.g. panels a, h and i).
There are also cases where the NN prediction misses the details of the model,
but still captures the general trend (e.g. panels t, o, and p). We will show
in the following section, that such differences are not crucial for the
interpretation of observed spectra in the context of fitting with Bayesian
methods.

Case 1left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc] Case
2left25pt+/2right-3pt-/2up-5pt-/2down3pt+/2[left,up][vc,hc]

Figure 4: Simulated blazar SEDs fitted with the GRU NN generated model. Simulated data are shown with markers, and with dashed magenta line is plotted the SED from the ATHE$\nu$A code. Dark- and light-shaded areas denote the range of GRU NN models with parameters from the 50% and 90% posterior distributions, respectively. Table 4: Values from GRU NN fitting to simulated SEDs with emcee. | Case 1 | Case 2
---|---|---
Parameter [unit] | True | Fit | True | Fit
$\log R$ [cm] | 16.60 | $16.04^{+0.61}_{-0.97}$ | 15.42 | $14.78^{+0.66}_{-0.53}$
$\log B$ [G] | -0.34 | $-0.49^{+0.36}_{-0.37}$ | -0.12 | $0.24^{+0.38}_{-0.38}$
$\log\gamma_{min}$ | 3.37 | $3.13^{+0.35}_{-0.24}$ | 2.31 | $2.20^{+0.31}_{-0.26}$
$\log\gamma_{max}$ | 5.55 | $5.51^{+0.25}_{-0.19}$ | 4.76 | $5.57^{+1.07}_{-0.98}$
$\log\ell_{e}$ | -4.57 | $-5.08^{+0.40}_{-0.49}$ | -1.82 | $-1.83^{+0.46}_{-0.29}$
$p$ | 2.20 | $2.76^{+0.17}_{-0.34}$ | 2.62 | $2.64^{+0.20}_{-0.22}$
$\delta$ | 1 | $1.3^{+0.2}_{-0.2}$ | 1 | $1.00^{+0.16}_{-0.14}$

##  3 Results

###  3.1 Simulated SEDs

A fundamental question is whether a trained NN can be used to recover the
physical parameters corresponding to an observed SED, provided one knows the
“correct” answer. The only way to test this is to use simulated SEDs based on
an SSC model with known parameters, which are then fitted by the NN; the
fitted values can be then compared to the true values of the SSC model to
address the question of whether the NN can constrain the parameters with
sufficient accuracy.

To create a simulated SED we select 25 points of an SSC model computed with
ATHE$\nu$A from a wide frequency range and add Gaussian noise with 0.1
standard deviation to the logarithmic flux values. We then attribute an error
to each flux point between 0.1 and 0.4 in logarithm to mimic statistical
uncertainties in the measurement. As blazars are observed in specific energy
bands, we select points that broadly correspond to radio, UV, X-ray and
$\gamma$-ray bands. By creating an incomplete SED with gaps in regions where
the solution may have large derivatives (e.g., spectral cutoffs), we increase
the difficulty for the NN fitting. We also note that the simulated flux points
may overshoot the baseline model of ATHE$\nu$A in certain frequency bands due
to the randomness of the simulation and the small number of selected points.
This mimics observational effects related to non-simultaneous observations and
intrinsic source variability.

We used the trained GRU NN to fit the data in Python with emcee444A Python
implementation of the Affine invariant Markov chain Monte Carlo (MCMC)
ensemble sampler. (Foreman-Mackey et al., 2013). For constructing the log-
likelihood function and fitting the data we followed an often-adopted
methodology in similar problems (e.g. Karaferias et al., 2023; Stathopoulos et
al., 2023). We also added a term $\ln{f}$ as a free parameter to account for
the systematic scatter and noise of the simulated data and the NN model. By
including $\ln{f}$ we get an excess variance compared to statistical
uncertainties, i.e. $\sigma_{\rm tot,\rm i}^{2}=\sigma_{\rm
i}^{2}+e^{2\ln{f}}$, where $\sigma_{\rm i}$ is the error of individual flux
points.

Here we present results for two characteristic cases. The simulated SEDs
(markers) together with the corresponding SED from the ATHE$\nu$A code (dashed
line) and the fitted models (shaded regions) are shown in Fig. 4. The corner
plots of the posterior distributions are presented in Figs. 6 and 7. The NN
model is able to successfully fit the simulated data and capture the overall
shape of the true model. In terms of parameters estimations all parameters
were consistent within 2$\sigma$ with the original values used for simulating
the data sets (see Table 4).

###  3.2 Observed SEDs

Figure 5: Observed SED of 3HSP J095507.9+355101 (colored markers) fitted with
the GRU NN generated model, with emcee (left) and UltraNest (right). Dark- and
light-shaded areas denote respectively the 50% and 90% confidence regions of
the fitted models.

We demonstrate the capabilities of the trained NN model to fit a real blazar
SED with emcee and UltraNest (Buchner, 2021), which employs nested sampling to
scan multi-modal posterior distributions. The latter allows us to better
explore possible degeneracies in multi-parameter problems.

For this application, we choose 3HSP J095507.9+355101 (Giommi et al., 2020;
Paliya et al., 2020), which is a BL Lac object at redshift $z=0.557$ (Paiano
et al., 2020; Paliya et al., 2020) that belongs to the rare class of extreme
blazars (Biteau et al., 2020). This blazar has been detected in GeV
$\gamma$-rays by the _Fermi_ Large Area Telescope (LAT) and is part of the
4FGL catalog (Abdollahi et al., 2020). 3HSP J095507.9+355101 is an ideal test-
bed because it combines measurements with small errors at lower energies with
several upper limits in the $\gamma$-rays. Another reason for our source
selection was that this target was also recently modeled with an emcee
adaptation of the leptonic module of LeHaMoC
555https://github.com/mariapetro/LeHaMoC, a Python-based lepto-hadronic code
that offers a factor of $\sim 30-100$ speed improvement compared to
ATHE$\nu$A, while producing consistent results (Stathopoulos et al., 2023).
Therefore, we can compare the efficiency and accuracy of SED fitting performed
with a fast leptonic code and a trained NN.

We used the same dataset as Stathopoulos et al. (2023). The quasi-simultaneous
observations in the UV, soft X-rays, and hard X-rays provide a detailed
picture of the low-energy part of the spectrum. On the contrary, the high-
energy part of the spectrum is less constrained observationally. For fitting
the data and constructing the log-likelihood function we followed Stathopoulos
et al. (2023), taking into account asymmetric errors and upper limits in the
Fermi data. For our application we defined uniform priors to cover all the
parameter space used to build the datasets, and a Doppler boost factor (in
logarithm) in the range [0,3]. We run emcee with 48 walkers and 10,000 steps
each. We started with a guess value that resembles the general spectral shape
and discard the first 1,000 steps of each chain as burn-in. We also ran
UltraNest with 400 live points, a maximum number improvement loops of 3 and a
target evidence uncertainty of 0.5.

The results of the fits are presented in Fig. 5. We find no noticeable
differences between the emcee and UltraNest results, demonstrating that the
SSC model has no multimodal posterior distributions (these are presented in
Figs. 8 and 9). When compared to the fitting results obtained with a leptonic
code, the posterior distributions of parameters obtained by the GRU NN and
LeHaMoC (Stathopoulos et al., 2023) are identical in terms of median values,
standard deviation and overall distribution shape as shown in the corner plot
of Fig. 8. Moreover, the physical processes that were neglected in the NN
training (i.e. synchrotron self-absorption and $\gamma\gamma$ pair production)
turn out to be negligible for this application, since both processes were
included in the fitting with LeHaMoC.

###  3.3 Efficiency of NNs

We briefly comment on the efficiency of the NN generated model in comparison
to the other numerical codes used in our tests.

The ATHE$\nu$A code, which was first presented in the mid-nineties
(Mastichiadis & Kirk, 1995), did not target fitting large samples of data. As
a result, the numerical scheme adopted in ATHE$\nu$A is not optimized for
speed but for accuracy (see also Cerruti et al., 2022, for comparison of
ATHE$\nu$A against other proprietary codes). The typical execution time of
ATHE$\nu$A for an SSC model is about 1-4 min. Moreover, the code may fail for
certain input parameters (when these lead to very large derivatives in the
PDEs), which cannot be determined a priori. For instance, about 4% of the
cases failed when creating the dataset for the NN training. LeHaMoC, on the
other hand, is more flexible for the generation of big datasets and its use by
fitting algorithms (Stathopoulos et al., 2023): it delivers a factor of $\sim
30-100$ improvement in speed compared to ATHE$\nu$A, and does not crash
regardless of the input parameters. When combined with the parallelization
option of emcee LeHaMoC can fit the SED of 3HSP J095507.9+355101 within 1 day
(using 16/32 CPU cores/threads). On the contrary, fitting of 3HSP
J095507.9+355101 using the GRU NN model was completed in approximately 20 min
on a single core, which is about 1,000 faster compared to LeHaMoC.

Therefore, implementing a NN like GRU offers a huge improvement in
computational time as it makes MCMC fitting with emcee in a matter of minutes.
When combined with nested sampling algorithms like UltraNest, the efficiency
of the sampling can drop significantly, as it is common for high-dimension
problems. For our application to 3HSP J095507.9+355101 about 20 million GRU NN
models were evaluated to search all the parameter space (see Table 1) within
21 hr. In Table 5 we list execution times for the various cases discussed in
this work. All computational tasks were performed with the same hardware (AMD
Ryzen 5950X). Quoted times were scaled to a single CPU core when parallel
computing was used. Times may vary depending on the hardware used, and
therefore the listed values in Table 5 are of comparative value.

Table 5: Typical CPU times for SSC computations.

| ATHE$\nu$A | LeHaMoC | GRU NN
---|---|---|---
Single SSC model | 1-4 min | 2-3 s | 3 ms
emcee* | – | 20 d | 20 min
UltraNest | – | – | 21 hr

  * *

For 48 walkers and 10,000 steps for each chain.

##  4 Conclusions

We have presented a proof-of-concept study of blazar SED modeling that relies
on the use of neural networks (NNs) and Bayesian inference. We have
demonstrated how computationally expensive numerical models can be substituted
with trained NNs. By testing 3 typical configurations of neurons we conclude
that a GRU NN offers optimal results for the astrophysical problem at hand. We
have tested the efficiency of our approach against simulated datasets and
observational data. Our results demonstrate the big leaps in computing
efficiency that can be achieved compared to state-of-the-art radiative
numerical codes. This in turn makes the process of fitting about 1,000 times
faster, enabling the use of not only typical MCMC methods but also nesting
sampling algorithms like in UltraNest.These findings will provide extra
motivation in the development and testing of different NNs or neurons and
experimenting with more complex configurations.

A natural next step would be the training of NNs against more complex
radiative models that include lepto-hadronic processes. A typical example is
the case of 3HSP J095507.9+355101 where lepto-hadronic models computed with
ATHE$\nu$A have only been tested by eye against the data (Petropoulou et al.,
2020) – see also Petropoulou et al. (2015a) for similar applications. Given
the large number of parameters (at least 11) and the execution time (tens of
minutes, even with faster codes like LeHaMoC) a statistical fitting with
Bayesian methods is prohibitive. In lepto-hadronic models the number of PDEs
that need to be solved increases from two to at least five. Moreover, the PDE
describing the evolution of each particle species in the source is non-
linearly coupled to the PDEs of other species. Due to the non-linear coupling
small change in the model parameters may lead to drastic changes in the
resulting SED. Complexity of lepto-hadronic models is further increased by the
fact that not all equilibrium solutions of the problem are constant in time.
There are regions of the high-dimension parameter space, as extensively
discussed by Mastichiadis et al. (2020), that produce oscillatory equilibrium
solutions (Petropoulou & Mastichiadis, 2012; Mastichiadis et al., 2005;
Petropoulou & Mastichiadis, 2018), known as limit cycles in non-linear
dynamics (Strogatz, 2000). The methodology presented here can also be extended
to other astrophysical problems, like spectral modeling of X-ray pulsars (see
examples of available numerical models Wolff et al., 2016; West et al., 2017;
Becker & Wolff, 2022) or GRBs (e.g. Rudolph et al., 2023).

Code Availability. – The results of our work will be made available in a
GitHub repository (https://github.com/tzavellas/blazar_ml) upon acceptance of
the paper, this includes; (a) the NN and accompanied code produced to train
them, (b) code for visualization of results in python and jupyter notebooks
with instructions, and (c) part of the ATHE$\nu$A datasets that can be used
for evaluation and plotting examples.

Note. – While we were finalizing this work, the paper by Bégué et al. (2023)
appeared. While the main concept and the radiation model are the same, the two
studies are independent and complementary. They focus on a different NN,
evaluation and modeling examples.

###### Acknowledgements.

M.P. acknowledge support from the Hellenic Foundation for Research and
Innovation (H.F.R.I.) under the “2nd call for H.F.R.I. Research Projects to
support Faculty members and Researchers” through the project UNTRAPHOB
(Project ID 3013). GV acknowledges support by H.F.R.I. through the project
ASTRAPE (Project ID 7802).

## References

  * Abadi et al. (2015) Abadi, M., Agarwal, A., Barham, P., et al. 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, software available from tensorflow.org
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33
  * Abiodun et al. (2018) Abiodun, O. I., Jantan, A., Omolara, A. E., et al. 2018, Heliyon, 4, e00938
  * Alarfaj et al. (2022) Alarfaj, F. K., Khan, N. A., Sulaiman, M., & Alomair, A. M. 2022, Symmetry, 14, 2482
  * Amato & Olmi (2021) Amato, E. & Olmi, B. 2021, Universe, 7
  * Becker & Wolff (2022) Becker, P. A. & Wolff, M. T. 2022, ApJ, 939, 67
  * Bégué et al. (2023) Bégué, D., Sahakyan, N., Dereli Bégué, H., et al. 2023, arXiv e-prints, arXiv:2311.02979
  * Biteau et al. (2020) Biteau, J., Prandini, E., Costamante, L., et al. 2020, Nature Astronomy, 4, 124
  * Blandford et al. (2019) Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467
  * Bloom & Marscher (1996) Bloom, S. D. & Marscher, A. P. 1996, ApJ, 461, 657
  * Buchner (2021) Buchner, J. 2021, The Journal of Open Source Software, 6, 3001
  * Celotti & Ghisellini (2008) Celotti, A. & Ghisellini, G. 2008, MNRAS, 385, 283
  * Cerruti et al. (2022) Cerruti, M., Kreter, M., Petropoulou, M., et al. 2022, in 37th International Cosmic Ray Conference, 979
  * Cerruti et al. (2015) Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910
  * Cho et al. (2014) Cho, K., Van Merriënboer, B., Bahdanau, D., & Bengio, Y. 2014, arXiv preprint arXiv:1409.1259
  * Chollet et al. (2015) Chollet, F. et al. 2015, Keras, https://keras.io
  * Chris Bishop et al. (1995) Chris Bishop, J., Bishop, C., Hinton, G., & Bishop, P. 1995, Neural networks for pattern recognition. Advanced texts in econometrics
  * Dimitrakoudis et al. (2012) Dimitrakoudis, S., Mastichiadis, A., Protheroe, R. J., & Reimer, A. 2012, A&A, 546, A120
  * Dimitrakoudis, S. et al. (2012) Dimitrakoudis, S., Mastichiadis, A., Protheroe, R. J., & Reimer, A. 2012, A&A, 546, A120
  * Finke et al. (2008) Finke, J. D., Dermer, C. D., & Böttcher, M. 2008, ApJ, 686, 181
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gao et al. (2017) Gao, S., Pohl, M., & Winter, W. 2017, ApJ, 843, 109
  * Gasparyan et al. (2022) Gasparyan, S., Bégué, D., & Sahakyan, N. 2022, MNRAS, 509, 2102
  * Giommi et al. (2020) Giommi, P., Padovani, P., Oikonomou, F., et al. 2020, A&A, 640, L4
  * Hillas (1984) Hillas, A. M. 1984, ARA&A, 22, 425
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. & Schmidhuber, J. 1997, Neural computation, 9, 1735
  * Hovatta & Lindfors (2019) Hovatta, T. & Lindfors, E. 2019, New A Rev., 87, 101541
  * Karaferias et al. (2023) Karaferias, A. S., Vasilopoulos, G., Petropoulou, M., et al. 2023, MNRAS, 520, 281
  * Maraschi et al. (1992) Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5
  * Mastichiadis et al. (2020) Mastichiadis, A., Florou, I., Kefala, E., Boula, S. S., & Petropoulou, M. 2020, MNRAS, 495, 2458
  * Mastichiadis & Kirk (1995) Mastichiadis, A. & Kirk, J. G. 1995, A&A, 295, 613
  * Mastichiadis & Kirk (1997) Mastichiadis, A. & Kirk, J. G. 1997, A&A, 320, 19
  * Mastichiadis et al. (2005) Mastichiadis, A., Protheroe, R. J., & Kirk, J. G. 2005, A&A, 433, 765
  * McCulloch & Pitts (1943) McCulloch, W. S. & Pitts, W. 1943, The bulletin of mathematical biophysics, 5, 115
  * Padovani et al. (2017) Padovani, P., Alexander, D. M., Assef, R. J., et al. 2017, A&A Rev., 25, 2
  * Paiano et al. (2020) Paiano, S., Falomo, R., Padovani, P., et al. 2020, Monthly Notices of the Royal Astronomical Society: Letters, 495, L108
  * Paliya et al. (2020) Paliya, V. S., Böttcher, M., Olmo-García, A., et al. 2020, ApJ, 902, 29
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., et al. 2019, in Advances in Neural Information Processing Systems 32 (Curran Associates, Inc.), 8024–8035
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Petropoulou et al. (2015a) Petropoulou, M., Dimitrakoudis, S., Padovani, P., Mastichiadis, A., & Resconi, E. 2015a, MNRAS, 448, 2412
  * Petropoulou & Mastichiadis (2012) Petropoulou, M. & Mastichiadis, A. 2012, MNRAS, 421, 2325
  * Petropoulou & Mastichiadis (2018) Petropoulou, M. & Mastichiadis, A. 2018, MNRAS, 477, 2917
  * Petropoulou et al. (2020) Petropoulou, M., Oikonomou, F., Mastichiadis, A., et al. 2020, ApJ, 899, 113
  * Petropoulou et al. (2015b) Petropoulou, M., Piran, T., & Mastichiadis, A. 2015b, MNRAS, 452, 3226
  * Readhead (1994) Readhead, A. C. S. 1994, ApJ, 426, 51
  * Rudolph et al. (2023) Rudolph, A., Petropoulou, M., Bošnjak, Ž., & Winter, W. 2023, ApJ, 950, 28
  * Stathopoulos et al. (2023) Stathopoulos, S. I., Petropoulou, M., Vasilopoulos, G., & Mastichiadis, A. 2023, arXiv e-prints, arXiv:2308.06174
  * Strogatz (2000) Strogatz, S. H. 2000, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering (Westview Press)
  * The Theano Development Team et al. (2016) The Theano Development Team, Al-Rfou, R., Alain, G., et al. 2016, arXiv e-prints, arXiv:1605.02688
  * Wang et al. (2019) Wang, S., Fan, K., Luo, N., et al. 2019, Nature Communications, 10, 4354
  * West et al. (2017) West, B. F., Wolfram, K. D., & Becker, P. A. 2017, ApJ, 835, 129
  * Wiecha & Muskens (2020) Wiecha, P. R. & Muskens, O. L. 2020, Nano Letters, 20, 329
  * Wolff et al. (2016) Wolff, M. T., Becker, P. A., Gottlieb, A. M., et al. 2016, ApJ, 831, 194

##  Appendix A Corner plots

Figure 6: Corner plot showing the posterior distributions of the parameters of
the GRU NN model presented in Fig. 4 (left). Contours denote Parameters of the
initial model are marked with magenta lines. Figure 7: Corner plot showing the
posterior distributions of the parameters of the GRU NN model presented in
Fig. 4 (right). Parameters of the initial model are marked with magenta lines.
Figure 8: Corner plot showing the posterior distributions of the GRU NN model
parameters for 3HSP J095507.9+355101as derived from emcee. Dashed lines in the
histograms indicate the median and 68 % range of values. With magenta color we
overplot the corner plot of the SSC model fitted by LeHaMoC as presented by
Stathopoulos et al. (2023). Figure 9: Corner plot showing the posterior
distributions of the GRU NN model parameters for 3HSP J095507.9+355101as
derived from UltraNest. Dashed lines in the histograms indicate the median and
68 % range of values.