# A thirteen-year-long broadband view of BL Lac

N. Sahakyan,1,2,3 P. Giommi2,4,5,6
1ICRANet-Armenia, Marshall Baghramian Avenue 24a, Yerevan 0019, Armenia
2ICRANet, P.zza della Repubblica 10, 65122 Pescara, Italy
3 ICRA, Dipartimento di Fisica, Sapienza Universita‘ di Roma, P.le Aldo Moro
5, 00185 Rome, Italy
4 Associated to Italian Space Agency, ASI, via del Politecnico snc, 00133
Roma, Italy
5Institute for Advanced Study, Technische Universität München,
Lichtenbergstrasse 2a, D-85748 Garching bei München, Germany
6 Center for Astro, Particle and Planetary Physics (CAP3), New York University
Abu Dhabi, PO Box 129188 Abu Dhabi, United Arab Emirates E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We report optical, ultraviolet, X-ray and $\gamma$-ray observations of BL Lac
objects prototype BL Lacertae carried out over a period of nearly 13 years,
between 2008 and 2021. The source is characterized by strongly variable
emission at all frequencies, often accompanied by spectral changes. In the
$\gamma$-ray band several prominent flares have been detected, the largest one
reaching the flux of $F_{\rm\gamma}(>196.7\>{\rm MeV})=(4.39\pm 1.01)\times
10^{-6}\>{\rm photon\>cm^{-2}\>s^{-1}}$, corresponding to the isotropic
luminosity of $\simeq 10^{47}\>{\rm erg\>s^{-1}}$. At the time of the first
major X-ray flare, on MJD 56268.65, a traditional harder-when-brighter trend
was observed. During the brightest flare, which occurred on MJD 59128.18, the
source instead exhibited a softer-when-brighter trend due to a shift of the
synchrotron peak to $10^{16}$ Hz, well into the HBL domain. The widely
changing multi-wavelength emission of BL Lacertae was investigated by fitting
one-zone and two-zone leptonic models to 520 high-quality and quasi-
simultaneous broad-band spectral energy distributions. The HBL behaviour
observed during the brightest X-ray flare is interpreted as due to the
emergence of synchrotron emission from freshly accelerated particles in a
second emission zone.

###### keywords:

quasars: individual: BL Lacertae– galaxies: jets – X-rays: galaxies- - gamma-
rays: galaxies

††pubyear: 2015††pagerange: A thirteen-year-long broadband view of BL Lac–A
thirteen-year-long broadband view of BL Lac

##  1 Introduction

Radio-loud Active Galactic Nuclei (AGNs) are characterized by two-sided narrow
relativistic jets that originate from the central supermassive black hole.
Blazars are the subclass of radio loud AGNs in which one of the jets happens
to make a small angle ($<10^{\circ}$) to the line of sight of the observer
(Urry & Padovani, 1995). These jets transport a large amount of power in the
form of particles, radiation and magnetic field and are strong sources of non-
thermal emission. Due to the small viewing angle and the relativistic motion
the emission in blazars is strongly Doppler boosted, a special situation that
makes these sources detectable up to large redshifts (e.g., Ackermann et al.,
2017; Sahakyan et al., 2020) and is responsible for the observed extreme
properties that characterizes them, like superluminal motion and rapid
variability across the electromagnetic spectrum. Historically blazars are
classified as BL Lacertae objects (BL Lacs), which exhibit an optical spectrum
that is completely featureless or at most shows very weak emission lines
(equivalent width ${\rm EW}\leq 5{\AA}$), and as flat spectrum radio quasars
(FSRQs) when the emission lines are stronger and quasar-like (Urry & Padovani,
1995). Blazars are generally assumed to be persistent sources, however
recently a case of a transient blazar, 4FGL J1544.3-0649, was observed. This
object was never detected in the X-ray and $\gamma$-ray bands until May 2017
when it raised above detactability and for a few months it became one of the
brightest X-ray blazars (Sahakyan & Giommi, 2021). If this is not an isolated
case, but rather a common phenomenon, it would significantly change our
current understanding of the extragalactic sky.
The broadband spectral energy distribution (SED) of blazars, in a ${\rm
log(\nu F\nu)}$ vs. ${\rm log(\nu)}$ representation, shows two prominent broad
components, one (low-energy component) peaking form far infrared frequencies
to X-ray energies and another (high energy component) peaking at MeV/GeV
energies. The peak of the low-energy component ($\nu_{\rm s}$) is used to
further classify blazars as high synchrotron peaked (HSP or HBL when $\nu_{\rm
s}>10^{15}$ Hz), intermediate synchrotron peaked (ISP or IBL when
$10^{14}<\nu_{\rm s}<10^{15}$ Hz), or low synchrotron peaked (LSP or LBL when
$\nu_{\rm s}<10^{14}$ Hz) objects (Padovani & Giommi, 1995; Abdo et al.,
2010). Sometimes the synchrotron peak can reach energies as high as $\sim$1
keV, ($\sim 2\times 10^{17}$ Hz) or beyond, showing what is considered to be
extreme behaviour, even for these highly peculiar sources (e.g. Giommi et al.,
1999; Costamante et al., 2001; Biteau et al., 2020). Such a high synchrotron
peak was first observed during a flare of Mkn 501 (Pian et al., 1998), and
subsequently in many other objects (e.g., Costamante et al., 2018; Sahakyan,
2020a). Independently of the location of the peak, the low-energy part of the
SED is generally interpreted as synchrotron emission from the relativistic
electrons in the jet. A proton synchrotron origin of the high energy end of
this component during X-ray flares has also been considered (Mastichiadis &
Petropoulou, 2021; Stathopoulos et al., 2021). The nature of the high energy
(HE; $>100$ keV) component is instead still under debate. Within one-zone
leptonic scenarios, the second component originates from inverse Compton
scattering of the synchrotron photons (SSC) by the electron population
producing the low-energy component (Ghisellini et al., 1985; Bloom & Marscher,
1996; Maraschi et al., 1992). Depending on the location of the emission
region, the photons external to the jet (e.g., photons from the disc, or those
reflected from the broad-line region or those from the infrared torus) can up-
scatter, producing the second component (external inverse Compton (EIC)
Błażejowski et al., 2000; Ghisellini & Tavecchio, 2009; Sikora et al., 1994).
On the other hand, the HE component can be also produced from the interaction
of relativistic protons either from their synchrotron emission (Mücke &
Protheroe, 2001) or from the secondary particles from pion decay (Mannheim,
1993; Mannheim & Biermann, 1989; Mücke & Protheroe, 2001; Mücke et al., 2003;
Böttcher et al., 2013). Recently, after associating TXS 0506+056 with the
IceCube-170922A neutrino event (IceCube Collaboration et al., 2018a, b;
Padovani et al., 2018) these scenarios have became more attractive. These
hadronic models also predict very high energy (VHE; $>100$ GeV) neutrinos
observable by the IceCube detector (Ansoldi et al., 2018; Keivani et al.,
2018; Murase et al., 2018; Padovani et al., 2018; Sahakyan, 2018; Righi et
al., 2019; Cerruti et al., 2019; Sahakyan, 2019; Gao et al., 2019).

Blazars, being powerful sources of strongly variable non-thermal emission, are
often targets of multiwavelength observations. The resulting data have been
accumulating over time enriching the archives with very valuable information
that can be used for detailed energy and time-domain investigations of the
origin of their emission. BL Lacertae (BL Lac) is one of these frequently
studied blazars; at $z=0.069$ it is a prototype of the BL LAC subclass of
blazars. BL Lac is usually classified as an LBL (Nilsson et al., 2018), but is
sometimes listed as an IBL (Ackermann et al., 2011). BL Lac is well known for
its prominent variability in a wide energy range, especially in the optical
(Larionov et al., 2010) and radio bands (Wehrle et al., 2016). BL Lac has been
a target of many multiwavelength campaigns ranging from the radio to the HE or
VHE $\gamma$-ray bands (Marscher et al., 2008; Raiteri et al., 2009, 2013;
MAGIC Collaboration et al., 2019; Weaver et al., 2020) which resulted in a
deep understanding of its properties in different bands. For example, in the
X-ray band, BeppoSAX observations in June 1999 showed that the 0.3-2 keV flux
of BL Lac doubled in $\sim 20$ min and the spectrum was concave with a very
hard component above 5-6 keV (Ravasio et al., 2002). In the $\gamma$-ray band,
the EGRET observations in 1995 showed an average $\gamma$-ray flux above 100
MeV of $(40\pm 12)\times 10^{-8}\>{\rm photon\>cm^{-2}\>s^{-1}}$ (Catanese et
al., 1997) which increased up to $(171\pm 42)\times 10^{-8}\>{\rm
photon\>cm^{-2}\>s^{-1}}$ during the flare in 1997 (Bloom et al., 1997).
Afterwards, Fermi Large Area Telescope (Fermi-LAT) observations showed that
during flaring periods the average $\gamma$-ray flux above 100 MeV can be as
high as $(2.3\pm 0.3)\times 10^{-6}\>{\rm photon\>cm^{-2}\>s^{-1}}$ (Cutini,
2012). VHE $\gamma$-rays above 1 TeV from BL Lac were initially reported by
the Crimean Observatory in 1998 (Neshpor et al., 2001) and later, in 2005, the
MAGIC telescope discovered a VHE $\gamma$-ray signal with an integral flux of
3% of the Crab Nebula flux above 200 GeV (Albert et al., 2007). The source is
flaring also in the VHE $\gamma$-ray band; for example, on June 28 2011, a
very rapid TeV $\gamma$-ray flare was detected by VERITAS when the integral
flux above 200 GeV reached roughly 125% of the Crab Nebula flux (Arlen et al.,
2013), or on June 15 2015 MAGIC detected a flare with a maximum flux of
$(1.5\pm 0.3)\times 10^{-10}{\rm photons\>cm^{-2}\>s^{-1}}$ and halving time
of $26\pm 8$ min (MAGIC Collaboration et al., 2019).
BL Lac shows a peculiar behavior both in terms of its classification and
interpretation of the observed broadband SED. First, the observation of
$H\alpha$ and $H\beta$ lines ($\sim 10^{41}\>{\rm erg\>s^{-1}}$) (Corbett et
al., 1996; Capetti et al., 2010) in different periods is quite unusual for
this type of blazars. This might indeed indicate a presence of a broad-line
region structure. On the other hand, the single-zone SSC models, usually
successful for explaining the TeV BL Lac spectrum, have a difficulty in
reproducing the variability of this source in different bands and taking into
account the emission in all the bands. When the spectrum extends to the VHE
$\gamma$-ray band or when a large Compton dominance is observed, the SED of BL
Lac can be modeled only by considering an EIC component added to SSC or by
using two-zone models (e.g., Bloom et al., 1997; Madejski et al., 1999;
Böttcher & Bloom, 2000; Abdo et al., 2011; MAGIC Collaboration et al., 2019).
This illustrates that different models/components are contributing in the
overall complex broadband spectrum of BL Lac.

Over the past decade, BL Lac was constantly monitored in the HE $\gamma$-ray
band by Fermi-LAT and AGILE (Bulgarelli et al., 2019) and frequently observed
in the optical/UV and X-ray bands by Neil Gehrels Swift Observatory (Gehrels
et al., 2004), (hereafter Swift). Together with the observations of other
instruments (NuSTAR, MAGIC, VERITAS, etc.) this resulted in the accumulation
of an extremely rich multi-frequency data set mapping both emission
components. The available data can be combined to build the broadband SED of
BL Lac in many different periods with (quasi) contemporaneous data. The
theoretical interpretation of these SEDs can help understanding the physical
processes that dominate in different periods. For example, a similar study of
the broadband emission of 3C 454.3 allowed to estimate the main parameters
describing the jet and emitting electrons as well as to investigate their
evolution in time (Sahakyan, 2021). Moreover, BL Lac was in active flaring
states in the optical and $\gamma$-ray bands in December 2020 and January 2021
(e.g., Marchini et al., 2021a, b; Cutini, 2021; D’Ammando, 2021b; Hosokawa et
al., 2021; D’Ammando, 2021c) when there was also observed the brightest
$\gamma$-ray flare from BL Lac (Mereu, 2020); on October 6 2020, the daily
averaged $\gamma$-ray flux of BL Lac was $(5\pm 1)\times 10^{-6}\>{\rm
photons\>cm^{-2}\>s^{-1}}$. The available multiwavelength data and the
extraordinary flaring activity of BL Lac in 2020/2021 motivated us to have a
new look on the origin of the broadband emission from it.

In this paper, analyzing the data observed by Fermi-LAT, Swift X-ray Telescope
(XRT) and Ultraviolet and Optical Telescope (UVOT) accumulated in the previous
thirteen years, we perform an intense broadband study of BL Lac. The paper is
organized as follows. The Fermi-LAT and Swift data collected for the analysis
and its reduction methods are described in Section 2. The spectral changes in
different bands and the broadband SED modeling is discussed in Section 3. The
discussion and conclusion are provided in Section 4.

##  2 Fermi-LAT Observations and Data Analysis

Since the launch of Fermi-LAT in 2008 BL Lac was constantly observed providing
unprecedented information on its emission in the $\gamma$-ray band. Fermi-LAT
is a pair conversion telescope sensitive to $\gamma$-rays in the energy range
from 100 MeV to 500 GeV. By default, it operates in all sky scanning mode,
mapping the entire $\gamma$-ray sky every three hours. Further details on
Fermi-LAT are given in Atwood et al. (2009).

For the current study, publicly available data accumulated between 04 August,
2008 and 01 March, 2021 are used (MET 239557417 - 636249605). The data have
been analysed by using Fermi ScienceTools version 1.2.1. The PASS8 “Source”
class events with a higher probability of being photons (evclass = 128,
evtype=3) in the energy range from 100 MeV to 500 GeV were analyzed using
P8R3_ SOURCE_ V2 instrument response function. The events were downloaded from
a region of interest (ROI) defined as a circular region with $12^{\circ}$
around the $\gamma$-ray position of BL Lac. The events are binned within a
$16.9^{\circ}\times 16.9^{\circ}$ square region into pixels of
$0.1^{\circ}\times 0.1^{\circ}$ and into 37 equal logarithmically spaced
energy bins. The model was created using the Fermi-LAT fourth source catalog
(4FGL-DR2; Ajello et al., 2020) where all sources within $17^{\circ}$ around
the target as well as the Galactic (gll_ iem_ v07) and the isotropic (iso_
P8R3_ SOURCE _ V2_ v1) diffuse emission components are included. The spectral
parameters of the background sources falling between RoI and RoI+$5^{\circ}$
were fixed to their catalog values, while the parameters of the other sources
and background models were left free. Binned likelihood analysis was applied
with gtlike tool to find the best matches between spectral models and the
data.

The source variability was investigated by dividing the entire period to three
day bins. During these short periods the source spectrum was modeled using a
power-law function, and the photon flux and index were estimated by applying
unbinned likelihood analysis with the appropriate quality cuts mentioned
above. The spectral indices of all sources (except BL Lac) are fixed to the
values obtained for the whole time period to reduce the uncertainties in the
flux estimations. Also, the normalization of both background components is
fixed to their average values obtained in $\sim 13$ years as no variability is
expected from them. In addition to the three day binned light curve, an
adaptively binned light curve was generated by adjusting the time bin widths
so as to attain 20% uncertainty in the flux estimation above an optimal energy
(see Lott et al., 2012, for details). This light curve with unequal time bins
has been proven to be particularly efficient for the identification of flaring
states (e.g., see Gasparyan et al., 2018; Sahakyan et al., 2018; Zargaryan et
al., 2017; Baghmanyan et al., 2017; Britto et al., 2016; Rani et al., 2013).

Figure 1: The multiwavelength light curve of BL Lac between August 04, 2008
and March 01, 2020. a) Adaptively binned $\gamma$-ray light curve ($>196.7$
MeV); b) 3-day binned $\gamma$-ray light curve ($>100$ MeV); c) Swift XRT
light curve in the 2.0-10 keV range; d) 0.3-10.0 keV X-ray photon index; e, f)
flux in V, B, U, W1, M2 and W2 filters.

The adaptively binned (E $>196.7$ MeV) and three day (E $>100$ MeV)
$\gamma$-ray light curves are shown in Fig. 1 panels a) and b) respectively.
Both light curves show the complex behaviour of BL Lac in the $\gamma$-ray
band; the mean $\gamma$-ray flux of the source is $4.46\times 10^{-7}\>{\rm
photon\>cm^{-2}\>s^{-1}}$ which increases up to $(4.39\pm 1.01)\times
10^{-6}\>{\rm photon\>cm^{-2}\>s^{-1}}$ (above $196.7$ MeV) observed on MJD
59231.34. In the considered thirteen years the source flux was above
$10^{-6}\>{\rm photon\>cm^{-2}\>s^{-1}}$ in total for $41.5$ days. The photon
index variation in time is investigated using a 3-day binned light curve. The
photon index is mostly soft with a mean value of $\Gamma_{\rm mean}=2.15$ but
occasionally it hardened to $\Gamma<2.0$. The hardest indexes of $1.48\pm
0.22$ and $1.61\pm 0.17$ were observed on MJD 57771.16 and 55782.16,
respectively.

###  2.1 Gamma-ray spectral analysis

The time evolution of the $\gamma$-ray emission is also investigated by
generating the SED at different times. When the SEDs are constructed for short
periods (e.g., three day bins or for the time intervals identified in the
adaptive bins lightcurve) the spectrum can be measured only up to the moderate
1 GeV range, not enough for a detailed study. Therefore, the Bayesian block
algorithm (Scargle et al., 2013) is used to divide the $\gamma$-ray light
curve into optimal intervals which are represented by an approximately
constant flux. By applying this algorithm, the points where the flux changes
from one state to another will be identified, providing the $\gamma$-ray
spectra of the source in different states. The Bayesian block algorithm
applied to the adaptively binned light curve divides the entire period into
218 intervals with a similar flux level. The shortest period is $5.81$ hours
during a flare while in the low emission state the longest period is $278.74$
days. The spectral analysis is applied by limiting the time for each interval
selected based on the Bayesian block. During the analysis, the spectrum of BL
Lac is assumed to be a power-law with spectral index and normalization left as
free parameters. The best matches between the spectral models and events are
obtained with an unbinned likelihood analysis implemented in gtlike. Depending
on source intensity the spectrum of BL Lac is obtained by separately running
the analysis for 4 or 7 energy bands of equal width in log scale.

Figure 2: The SED of BL Lac in different periods. Left panel: The periods with
soft X-ray spectrum as modeled within two-zone leptonic scenario. Right panel:
The usual hard X-ray spectrum (Swift XRT and Nustar) observed on MJD 59133.88.

###  2.2 Swift XRT

During the considered period, the Swift satellite (Gehrels et al., 2004)
observed BL Lac 610 times with single exposures ranging from 1.13 to 16.46 ks.
All the data were downloaded and processed using Swift_xrtproc automatic tool
for XRT data analysis developed within the Open Universe Initiative (Giommi et
al., 2021). This tool automatically downloads the raw data and processes it
using the XRTPIPELINE task adopting standard parameters and filtering
criteria. For each observation, it extracts the source events from a circle
with a radius of 20 pixels centered at the position of the source, while the
background counts are taken from an annular ring centered at the source. The
tool applies also pile-up correction when the source count rate is above
$0.5\>{\rm counts\>s^{-1}}$. Then it loads the ungrouped data in XSPEC
(version 12.11) for spectral fitting using Cash statistics (Cash, 1979),
modeling the source spectrum as a power-law and a log-parabola model with the
Galactic absorption column density fixed to $1.75\times 10^{21}\>{\rm
cm^{-2}}$.

The 2-10 keV X-ray flux variation is shown in Fig. 1 panel c). The baseline
flux is around $2\times 10^{-11}\>{\rm erg\>cm^{-2}\>s^{-1}}$ although small
amplitude changes are visible in different observations. In three periods, MJD
56300, MLD 59140 and MJD 59235, the flux substantially increased reaching the
maximum of $(1.54\pm 0.04)\times 10^{-10}\>{\rm erg\>cm^{-2}\>s^{-1}}$ on MJD
59128.18. This is the historical highest flux of BL Lac in the soft X-ray
band.

The X-ray photon index in different observations is shown in Fig. 1 panel d).
Most of the time, the photon index is hard and $\leq 2.0$ implying that the
X-ray emission is due to the rising part of the HE component in the SED of BL
Lac. However, the photon index undergoes interesting modifications reaching
$2.0$ which corresponds to a flat distribution in $\nu f\nu$ vs $\nu$
representation. For example, such tendency can be noticed after the X-ray
flare around MJD 56300. In 18 observations the X-ray photon index is $>2.2$
(considering only the observations when the number of counts was $>100$) which
is unusual for BL Lac and more typical of HBL blazars. The optical/UV and
X-ray spectrum of BL Lac during such changes is shown in Fig 2. The X-ray
component started to soften starting from MJD 59124.73 when an index of
$\Gamma_{\rm X}=2.37\pm 0.10$ was observed. Then, the softest index of
$\Gamma_{\rm X}=2.56\pm 0.03$ was observed on MJD 59128.18 during the
brightest X-ray emission state (light blue squares in Fig. 2). In this period
the optical/UV flux increased substantially as well, showing that the low-
energy component now extends to the X-ray band. Such soft X-ray emission with
$\Gamma_{\rm X}=2.48\pm 0.07$ and a slightly lower flux of $(4.01\pm
0.40)\times 10^{-11}\>{\rm erg\>cm^{-2}\>s^{-1}}$ (red circle in Fig.2) was
also observed on MJD 59128.91. In the next two observations (MJD 59129.90 and
59131.83), the X-ray flux was constantly decreasing and the photon index was
$\Gamma_{\rm X}=2.22-2.37$. The soft component with $\Gamma_{\rm X}=2.53\pm
0.11$ and a low X-ray flux of $(1.23\pm 0.18)\times 10^{-11}\>{\rm
erg\>cm^{-2}\>s^{-1}}$ was still observed on MJD 59132.88 (magenta triangles
in Fig. 2). However, this component fades in the next observations (on MJD
59133.81) and the usual component with a hard index ($\Gamma_{\rm X}=1.61\pm
0.11$) rises. There were additional periods when prolonged moderate softening
in the X-ray band was observed; for example, on MJD 58685.98 and 58686.90 and
between MJD 58740.42-58741.41 the X-ray photon index was $\Gamma_{\rm
X}=2.2-2.3$ with a flux between $(0.49\pm 1.11)\times 10^{-11}\>{\rm
erg\>cm^{-2}\>s^{-1}}$.
The X-ray flux evolution was further investigated by comparing it with the
photon index in different states. When considering the entire observational
period with diverse X-ray emission properties, any trend (if present) will be
smoothed out. For this reason, the X-ray photon index versus the flux was
investigated by selecting the periods around two major flares visible in Fig.
1; namely within MJD 56160-56350 and MJD 59000-59350. The results are shown in
Fig. 3. The linear-Pearson correlation test applied to the data during the
first flare (MJD 56160-56350) yields $-0.56$ implying a negative correlation
between the flux and photon index, i.e., when the source gets brighter, the
photon index decreases (hardens). This behaviour has already been observed for
many flaring blazars (e.g. Giommi et al., 1990). On the other hand, the
linear-Pearson test results in $0.44$ for the second flare. This implies that
during the X-ray flare the photon index softens, so a softer-when-brighter
trend is observed. This shows that two major flares observed in the X-ray band
for BL Lac are different by their nature and are caused by different
processes.

Figure 3: BL Lac X-ray photon index versus the flux during two major X-ray
flares.

###  2.3 Swift UVOT

UVOT observed BL Lac in all six filters, V (500-600 nm), B (380-500 nm), U
(300- 400 nm), W1 (220-400 nm), M2 (200-280 nm) and W2 (180–260 nm)
simultaneously with the XRT. All single observations of BL Lac were downloaded
and reduced using HEAsoft version 6.27 with the latest release of HEASARC
CALDB. The data are reduced using standard procedures, by selecting source
counts from a circular region of $5^{\prime\prime}$ around the source, while
the background counts were estimated from a $20^{\prime\prime}$ region away
from the source. uvotsource tool was used to derive the magnitudes which were
converted to fluxes using the conversion factors provided by Poole et al.
(2008) which then were corrected for extinction using the reddening
coefficient $E(B-V)$ from the Infrared Science Archive
111http://irsa.ipac.caltech.edu/applications/DUST/.

The optical/UV flux evolution in time is shown in Fig. 1 panel e) and f)
separating the flux in V, B, U and W1, M2 and W2 filters. The source flux is
relatively constant at the level of a few times $10^{-11}\>{\rm
erg\>cm^{-2}\>s^{-1}}$ up to MJD $56500$. A flaring activity occurred around
MJD 56617-56622 when the flux in all filters exceeded $10^{-10}\>{\rm
erg\>cm^{-2}\>s^{-1}}$. The major flaring activity started on MJD 59072
coinciding with those in the X-ray and $\gamma$-ray bands and the baseline
flux level was above $10^{-10}\>{\rm erg\>cm^{-2}\>s^{-1}}$. The maximum flux
of $(5.80\pm 0.14)\times 10^{-10}\>{\rm erg\>cm^{-2}\>s^{-1}}$ was observed in
the V band on MJD $59249.26$. The maximum flux in B, U and M2 filters was also
above $5\times 10^{-10}\>{\rm erg\>cm^{-2}\>s^{-1}}$ while that in W1 and W2
was around $4\times 10^{-10}\>{\rm erg\>cm^{-2}\>s^{-1}}$.

###  2.4 Archival data

To achieve as much as possible a complete view of the broad-band emission from
BL Lac we have also considered all the available multi-frequency archival
measurements alongside with the data from Fermi-LAT, Swift XRT and UVOT. These
include a) optical data monitoring from the ASAS-SN Sky Patrol web site
222https://asas-sn.osu.edu/ (Kochanek et al., 2017), b) NuStar data from the
observations of BL Lac on 11 December 2012, 14 September 2019 and 11 October
2020 (from Middei et al. 2021, submitted), and c) any other multi-frequency
measurements available from the VOU_BLazars tool (Chang et al., 2020) and the
ASI Space Science Data Center (SSDC) archive 333https://www.ssdc.asi.it. In
addition we also considered the observations of BL Lac carried out by VERITAS
on June 11, 2011 (MJD 55740 Arlen et al., 2013) and on October 5, 2016 (MJD
57697 Abeysekara et al., 2018) and by MAGIC between June 15 and 28, 2015 (MJD
57188-57201 MAGIC Collaboration et al., 2019). The combination of all these
data sets results in an unprecedented amount of observations of BL Lac
covering the spectrum from radio frequencies to HE and VHE $\gamma$-ray bands
over a period of nearly 13 years, from 2008 to 2021.

##  3 Modeling the SEDs

In this section we use the data assembled as described above to investigate
the evolution of the broadband spectrum of BL Lac between 2008 and 2021. To
this end we have generated a large number of quasi contemporaneous SEDs by
plotting the computed $\gamma$-ray spectra together with the data available in
all other energy bands in each of the Bayesian intervals defined in Sec. 2.1.
To illustrate the temporal evolution of the broad band emission from BL Lac in
a visually effective way we have combined these SEDs to form an animation that
is available at the following link: youtube.com/L1yT105UGYM. Flux changes in
the optical/UV, X-ray and $\gamma$-ray bands are evident with variations of
almost two orders of magnitude. In the $\gamma$-ray band, the spectrum hardens
together with the flux amplification, resulting in a shift of the peak of the
second component to higher energy values. During the brightest and softest
X-ray state (on MJD 59128.91), the low-energy SED component extended to the
X-ray band as a consequence of a significant change of the location of the
synchrotron peak from the usual $\sim 10^{14}$ Hz to $\sim 10^{16}$ Hz, well
into the HBL regime (Padovani & Giommi, 1995; Abdo et al., 2010). Such a large
modification, never observed before in BL Lac , marks the extraordinary nature
of this flare, which has been studied also by D’Ammando (2021a) and Prince
(2021).

The classical double-humped SED of BL Lac is usually interpreted within
leptonic scenarios. The EGRET observations of BL Lac (Sambruna et al., 1999;
Madejski et al., 1999) revealed that modeling of the HE data requires a
component that extends beyond the SSC radiation generated in a single emission
zone: one-zone leptonic modeling requires a very high $\delta$ or an extended
emission region. Since then the SED of BL Lac has been conventionally modelled
within two-zone scenarios (e.g. Abdo et al., 2011; MAGIC Collaboration et al.,
2019) or assuming inverse Compton scattering of external (EIC) photons (e.g.,
Madejski et al., 1999; Böttcher et al., 2013). External Compton scenarios are
favoured in BL Lac, considering the detection, although weak, of the $H\alpha$
line (Corbett et al., 1996; Capetti et al., 2010), which points to the
presence of a broad-line region (BLR). Even if this BLR is not large enough to
absorb VHE $\gamma$-rays through $\gamma-\gamma$ interaction (e.g., Donea &
Protheroe, 2003), it can provide targets for inverse Compton up-scattering.
For example, by modeling the SED of BL Lac, Abdo et al. (2011) showed that the
SSC+ERC scenario provides reasonable modeling of the data also during the low
state and the inverse Compton scattering of the BLR-reprocessed radiation
strongly dominates over that directly from the disc.

In an effort to understand the processes dominating in the jet in different
physical conditions we investigated the broadband emission from BL Lac by
modeling the SEDs observed in different periods. From the SED/light curve
animation discussed above we have selected all periods with sufficient multi-
wavelength data, typically those with flux measurements in at least the
optical/UV, X-ray and $\gamma$-ray bands. This allowed us to assemble 520
high-quality and quasi-simultaneous SEDs representing BL Lac in a variety of
emission states. All these SEDs are modeled assuming that the emission region
(’the blob’) is a sphere with radius $R$, including a magnetic field of
intensity $B$ and a population of relativistic electrons following an energy
distribution defined by a power law with an exponential cutoff, as expected
from shock accelerations:

$N(\gamma^{\prime})=N^{\prime}_{e}\gamma^{\prime-p}\>Exp(-\gamma^{\prime}/\gamma^{\prime}_{\rm
cut})$ (1)

for $\gamma^{\prime}>\gamma^{\prime}_{\rm min}$ where $\gamma^{\prime}_{\rm
min}$ and $\gamma^{\prime}_{\rm cut}$ are the minimum and cut-off energy of
the electrons, respectively. It is assumed that the emission region is located
inside the BLR and the low energy SED component is interpreted as synchrotron
emission of relativistic electrons, while the second SED component is due to
inverse Compton up-scattering of photon fields from the jet itself (SSC model
e.g., Maraschi et al., 1992; Bloom & Marscher, 1996) and those reflected from
the BLR clouds (EIC BLR; Sikora et al., 1994). The BLR is assumed to be
spherical shell with an average radius of $R_{\rm BLR}=7\times 10^{16}$ cm and
lower and upper boundaries of $0.9\times R_{\rm BLR}$ and $1.2\times R_{\rm
BLR}$, respectively (Donea & Protheroe, 2003). The BLR reflects 10% of the
disc luminosity whose emission is approximated as a mono-temperature black
body with a luminosity of $L_{\rm d}=3\times 10^{43}\>{\rm erg\>s^{-1}}$. This
luminosity was estimated with a requirement that the disc component does not
overproduce the optical/UV data in any period. We note that $R_{\rm BLR}$ and
$L_{\rm d}$ define the density of the external photon fields, so their small
changes do not affect the results and only will result in moderate changes in
the normalization of the electrons.
Our 520 SEDs represent an ample variety of different states of BL Lac and in
some periods the simple one-zone model described above cannot explain the
observed data. For example, when the X-ray spectrum softens neither the
synchrotron component (defined by the optical/UV data) nor SSC component which
has a rising shape cannot account for the X-ray flux. In these cases the SEDs
were modeled considering two-zone scenarios, assuming that one region is
within the BLR and the other is outside (e.g., Tavecchio et al., 2011).
The broadband SEDs model fitting was carried out using the open source package
JetSet (Massaro et al., 2006; Tramacere et al., 2009, 2011; Tramacere, 2020).
The free model parameters ($p$, $\gamma_{\rm cut}$, $\gamma_{\rm min}$,
$\delta$, $R$ and $B$) are constrained by using the Minuit optimizer and then
improved by Markov Monte Carlo Chain sampling of their distributions. In
principle, $R$ can be constrained using the variability time scales from the
relation $R\leq\delta\>t\>c/(1+z)$, but instead of using fixed $R$ for all the
periods we considered that it varies between $10^{15}$ and $10^{16}$ cm, which
corresponds to the day/sub-day variability that is often observed in BL Lac.
When two-zone modeling was considered, to reduce the number of free
parameters, $R$ was fixed for both emitting regions. Since the TeV or X-ray
observations in the bright states reveal that the flux varies on minute
scales, implying that the emitting region outside BLR should be very compact,
$R=10^{15}$ cm was assumed. On the other hand, the optical/UV and $\gamma$-ray
fluxes vary albeit not on such short scales, so for the blob within the BLR
$R=10^{16}$ cm was used. Also, both emitting regions were assumed to have the
same Doppler boosting factor. In principle, because of orientations those
regions can have different Doppler boosting factors which, however, will
introduce an additional free parameter.

###  3.1 SEDs modeling results

Figure 4: The broadband SEDs of BL Lac observed in different emission states.
a) The emission components obtained from the modeling of all SEDs with
contemporaneous data collected during 2008-2021. b) The SED of BL Lac in the
normal emission state. c) BL Lac SED modeling with VHE $\gamma$-ray data from
MAGIC observations (MAGIC Collaboration et al., 2019). d) The SED of BL Lac
during the brightest X-ray period modeled within two-zone scenario.

The animation of our 520 high-quality and quasi-contemporaneous SEDs of BL Lac
together with the corresponding modeling is available at the following link
youtube.com/watch?v=f3a5CGukbbE. In this animation, the sum of all model
components is plotted as a solid blue line, while the SSC and EIC components
appear as green and orange lines, respectively. The disc emission,
approximated as a black body with intensity that is always below the
synchrotron component, is shown in magenta. Fig. 4 shows the emission
components in all 520 SEDs (panel a) and some frames representing special
states (panels b-d). The optical/UV data constrain the tail of the synchrotron
component which peaks at $\sim 10^{14}$ Hz and, despite large flux
variability, it remains almost unchanged as can be seen from Fig. 4 panel a)
(blue curves; the bright and soft X-ray periods were not considered). The SSC
emission of the synchrotron emitting electrons starts to dominate around
$10^{17}$ Hz extending up to $\sim 10^{23}$ Hz (green dot-dashed lines in Fig.
4 panel a) while at higher frequencies EIC of BLR photons dominates (orange
dot-dot-dashed lines in Fig. 4 panel a). The change of intensity of these
components show high-amplitude variability of BL Lac emission in the
optical/UV, X-ray and $\gamma$-ray bands. The variability in the radio band
cannot be tested, as the data are missing for most of the cases. Moreover, the
radio emission can be produced, with significant time-lags, by the low-energy
electrons in extended regions which is not associated with the emission in
other bands.

Figure 5: The evolution of electron parameters obtained from the modeling of
520 quasi-contemporaneous broadband SEDs of BL Lac.

The modeling provides estimates of the physical parameters describing the
emission from BL Lac and allows us to investigate their changes in time. The
evolution of the $p$, $\gamma_{\rm min}$, $\gamma_{\rm cut}$, $\delta$ and $B$
parameters is shown in Fig. 5 panels a) to e). The minimum energy of electrons
$\gamma_{\rm min}$ is mostly below $10$ (Fig. 5 panel a) implying that even
lower-energy electrons are efficiently accelerated. The power-law index ($p$)
is mostly within $1.3-2.3$ (Fig. 5 panel b) and is defined by fitting the
X-ray data with the SSC component. Its variation is in accordance with the
changes of the X-ray photon index shown in Fig. 1 (panel d). The variation of
the cut-off energy is shown in Fig. 5 panel c) which is defined mostly by the
optical/UV and sometimes by $\gamma$-ray data; it is within $301.4-2457.5$.
Despite such change in the cut-off energy, the low- and high- energy peaks in
the SED do not deviate to higher frequencies; the peak positions are
determined by the product of $\gamma_{\rm cut}$ and $B$. In this case, the
magnetic field (see Fig. 5 panel c) varies in the relatively narrow range of
$B=1.08-12.75$ G but, in particular, when a high $\gamma_{\rm cut}$ was
estimated, $B$ was around its lower level (see Fig. 5). It is also interesting
to investigate the evolution of the Doppler factor $\delta$ (Fig. 5 panel d).
This parameter remained mostly below 20 but it reached the maximum value of
$\delta=63$ on MJD 59117.4 when the major $\gamma$-ray flares were observed.

An example of SED modeling during the period of MJD 54989.3-55014.9, when the
source was not flaring in any band, is shown in Fig. 4 panel b). The
synchrotron emission of the accelerated electrons explains the archival (gray)
and optical/UV (light blue) data, while the X-ray emission (blue) is due to
SSC emission. This SSC component extends up to $10^{20}$ Hz and the emission
in the $\gamma$-ray band is dominated by the IC scattering of BLR photons.
This shows that even in the quiescent state of the source, the external photon
field (BLR) is necessary to explain HE $\gamma$-ray data. The power-law index
of the emitting electrons is $p=1.92\pm 0.02$ and $\gamma_{\rm cut}=678.83\pm
22.9$, while the emitting regions size is $R=4.04\times 10^{15}$ cm which
moves with a Doppler factor of $\delta=12.94$. The system is close to the
equipartition $U_{\rm e}/U_{\rm B}\simeq 1.5$ with $B=2.92\pm 0.09$ G.

An example of the BL Lac SED modeling with VHE $\gamma$-ray data from MAGIC
observations is shown in Fig. 4 panel c). When considering a one-zone
scenario, the optical/UV data with a decreasing trend constrains the cut-off
energy of the emitting electrons and the IC scattering of $1-5$ eV synchrotron
or BLR photons on them will only reach the MeV/GeV band, unable to account for
the VHE $\gamma$-ray data. These periods are modeled within two-zone
scenarios, considering the emitting regions inside and outside the BLR. The
emission observed in the radio to HE $\gamma$-ray bands is dominated by that
from the region within the BLR, and the VHE $\gamma$-ray data are explained by
the SSC emission from the compact region. In the extended region ($R=10^{16}$
cm), $B=1.87\pm 0.08$ G and the emitting electrons have
$\gamma^{\prime-1.29\pm 0.02}$ distribution with $\gamma^{\prime}_{\rm
cut}=1082.77\pm 36.1$. Instead, in the compact region, requiring that its
synchrotron emission is lower than that from the other region, $B=3.32\times
10^{-2}$ G is estimated, implying that the electrons can be accelerated to
higher energies, i.e., $\gamma^{\prime}_{\rm cut}=1.39\times 10^{5}$ was
estimated in this case. The contribution of these electrons with $p=2.17$
starts to dominate above $\sim 30$ GeV when the spectrum of IC scattering of
BLR photons decreases, explaining the data observed by MAGIC.

Similar two-zone models are also required when the soft component in the X-ray
band is observed. The synchrotron component defined by the available
optical/UV data effectively extends up to $10^{17}$ Hz ($\sim 400$ eV), unable
to account for the observed X-ray data. Thus, in the X-ray band an additional
component is dominating. Among the selected periods, 13 SEDs with a soft X-ray
spectrum were modeled within the two-zone scenario; an example of SED modeling
when the brightest X-ray emission was observed is shown in Fig. 4 panel d).
The Doppler boosting factor of both emitting regions is $\delta=16.24$ but
they are filled with different distributions of electrons. For example, the
X-ray emitting electrons (the region outside BLR) have $\gamma^{-1.50}$
distribution above $\gamma_{\rm min}=760.92$ with a cut-off energy of
$1.13\times 10^{4}$. Instead, the electrons in the other region have a softer
distribution with $1.97\pm 0.16$ and are accelerated only up to moderate
energies of $\gamma_{\rm cut}=1554.84\pm 214.0$. The magnetic field in the
region outside the BLR is stronger ($B=7.12$ G) than that in the other region
($B=1.67\pm 0.22$ G) which is because i) the first region has a smaller radius
and ii) the synchrotron emission should be at higher frequencies, reaching up
to the X-ray band.

##  4 Discussion and Conclusion

We performed a comprehensive investigation of the large and complex luminosity
and spectral variability of BL Lac using the data taken by Fermi-LAT, Swift-
XRT and Swift-UVOT between 2008 and 2021. In all the energy bands considered,
the source shows multiple periods when the flux exceeds its average level by
substantial amounts. The maximum $\gamma$-ray flux of
$F_{\rm\gamma}(>196.7\>{\rm MeV})=(4.39\pm 1.01)\times 10^{-6}\>{\rm
photon\>cm^{-2}\>s^{-1}}$ was observed on MJD 59231.3, associated with a flat
spectral slope with photon index of $2.03\pm 0.21$. This implies an energy
flux, ${\rm\epsilon F_{\epsilon}}$, of $9.39\times 10^{-9}\>{\rm
erg\>cm^{-2}\>s^{-1}}$ in the 0.1-300 GeV energy range, corresponding to an
isotropic $\gamma$-ray luminosity of $L_{\rm\gamma}=4\>\pi\>d_{\rm
L}^{2}\epsilon F_{\epsilon}\simeq 10^{47}\>{\rm erg\>s^{-1}}$ for a $\sim 307$
Mpc distance. Assuming a Doppler factor of $\delta=20$ this corresponds to
$L_{\rm\gamma}/\delta^{2}\simeq 2\times 10^{44}\>{\rm erg\>s^{-1}}$ in the
proper frame of the jet. This value largely exceeds the disc luminosity
estimated under any reasonable assumption (e.g., $3\times 10^{43}\;{\rm
erg\>s^{-1}}$ to not overproduce the observed optical/UV data) implying
extreme energetics during the $\gamma$-ray flares (e.g., Ghisellini et al.,
2014). Unlike the flux, the photon index does not usually vary significantly,
although it occasionally hardens to values $<2.0$. Such hardening, for
example, was noticed after the brightening observed on MJD 59247.4 when the
$\gamma$-ray flux measured within three days was $(3.32\pm 0.15)\times
10^{-6}\>{\rm photon\>cm^{-2}\>s^{-1}}$; the photon index was within
$1.79-1.89$ during MJD 59247.4-59256.4.

An interesting behaviour of BL Lac was observed in the X-ray band where it
showed two major flaring activities on MJD 56268.65 and 59128.18 with
different characteristics. During the first flare, when the 2-10 keV flux
reached it maximum value of $(7.59\pm 0.52)\times 10^{-11}\>{\rm
erg\>cm^{-2}\>s^{-1}}$ and the traditional harder-when-brighter trend was
observed. In the X-ray band, this is a known behaviour for blazars (Giommi et
al., 1990; Zhang et al., 1999; Giommi et al., 2005). On the other hand, during
the second flare the linear-Pearson test resulted in $0.44$ clearly showing a
softer-when-brighter trend. Such a pattern was also observed in the X-ray
emission of OJ 287 (e.g., Kapanadze et al., 2018; Prince et al., 2021; Giommi
et al., 2021). During this flare, the spectral index of $2.56\pm 0.03$ was the
softest ever observed, very different from the values normally observed in BL
Lac (typically $\Gamma_{\rm X}\leq 2.0$). Significant softening of the X-ray
spectrum, with $\Gamma_{\rm X}\geq 2.2$, has been observed on 18 additional
occasions. When this soft component is associated with a high X-ray flux
(i.e., $F_{\rm X(2-10\>keV)}\geq 10^{-11}\>{\rm erg\>cm^{-2}\>s^{-1}}$) the
peak of the SED low-energy component reaches frequencies of ${\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}}10^{15}$ \- $10^{16}$ Hz, instead
of the usual $\sim 10^{14}$ Hz, temporarily placing BL Lac into the domain of
HBL blazars. During these flares the source X-ray emission is driven by a new
HBL-like component, which significantly differs from the usual X-ray spectrum
of BL Lac. Our modeling shows that this component may come from a separate
emission zone with specific properties, like the size of the emission region,
the population of electrons, etc. This component is present during the period
MJD 59124.73-59132.88 when the flux rises from $F_{\rm X(2-10\>keV)}=(2.72\pm
0.31)\times 10^{-11}\>{\rm erg\>cm^{-2}\>s^{-1}}$ with $\Gamma_{\rm X}=2.37\pm
0.10$, and reaches $F_{\rm X(2-10\>keV)}=(1.54\pm 0.04)\times 10^{-10}\>{\rm
erg\>cm^{-2}\>s^{-1}}$ ($\Gamma_{\rm X}=2.56\pm 0.03$) on MJD 59128.18. This
component fades on MJD 59132.88 when the flux decreases to $F_{\rm
X(2-10\>keV)}=(1.23\pm 0.18)\times 10^{-11}\>{\rm erg\>cm^{-2}\>s^{-1}}$
($\Gamma_{\rm X}=2.53\pm 0.11$). The SEDs observed in this period are shown in
Fig. 2. This new component clearly goes beyond the synchrotron radiation
constrained by the optical/UV data and is interpreted as synchrotron emission
from the second region (dashed line) containing much more energetic particles.
For example, the electrons should be accelerated up to $\gamma_{\rm
cut}=1.76\times 10^{4}$ with $p=1.15$ to explain the data observed on MJD
59124.73 (red line). Similar parameters obtained from the modeling of the SED
on MJD 59128.18 are $\gamma_{\rm cut}=1.13\times 10^{4}$ and $p=1.50$ (blue
line) but the magnetic field is $7.12$ G, significantly higher compared to the
previous case ($1.65$ G). Such a large magnetic field is required because of
the increase in the X-ray flux ($\sim 5.7$ times) which cannot be explained by
changing $\gamma_{\rm cut}$; $p$ and $\gamma_{\rm cut}$ are also constrained
by the $\gamma$-ray data. The X-ray flux variation impacts the magnetic field
which decreases to $B=0.27$ G on MJD 59132.88 when the soft X-ray component
was with a low state. The electrons in the emitting region are still energetic
with $\gamma_{\rm cut}=5.60\times 10^{4}$ but their contribution starts to be
subdominant. In fact, the usual hard X-ray component which is interpreted as
SSC radiation from the blob inside the BLR dominates already on MJD 59133.81
(Fig. 2 right panel). This implies that either the acceleration/injection of
the energetic electrons is not efficient anymore or due to the drop of the
magnetic field these electrons cool down on longer time scales; for example,
when $B=0.1$ G the cooling of $\gamma=5\times 10^{4}$ electrons in the
observer frame is $t_{\rm cool}=6\>\pi\>m_{e}\>c\>(1+z)/\sigma_{\rm
T}\>B^{2}\>\delta\gamma=1.3\>(B/0.1\>G)^{-2}\>(\delta/15)^{-1}\>(\gamma/5\times
10^{4})^{-1}$ day.

The composition of the emitting region ($U_{\rm e}/U_{\rm B}$) changes for the
considered periods. A slightly particle dominated region ($U_{\rm e}/U_{\rm
B}\simeq 12$) is necessary to explain the SED on MJD 59124.73 while it should
be magnetically dominated with $U_{\rm e}/U_{\rm B}\simeq 0.43$ to explain the
data observed on MJD 59128.18 and it is strongly particle dominated on MJD
59132.88 with $U_{\rm e}/U_{\rm B}\simeq 2.8\times 10^{3}$. This indicates
that the X-ray flare was caused by the injection of new energetic particles
and by sudden increase of the magnetic field. This is in agreement with the
observation of the softer photon index during the bright X-ray state; due to
the high magnetic field the electrons cool faster forming a soft spectrum with
the increase of intensity.

The emission from a second region is also required when BL Lac was flaring in
the VHE $\gamma$-ray band, but under different conditions for the emitting
region. It still contains energetic electrons with $\gamma_{\rm cut}\simeq
10^{5}$ but the magnetic field is low ($0.02-0.1$ G), so the region is
strongly particle dominated with $U_{\rm e}/U_{\rm B}\geq 10^{3}$. This makes
the emission from these electrons significant in the VHE $\gamma$-ray band
with no significant contribution at lower energies.

An emission from the second more energetic population of electrons is required
when modeling the flaring components in the X-ray and VHE $\gamma$-ray bands.
An excess of a new component in the X-ray band was already noticed in previous
observations of BL Lac in 2007-2008: the XMM Newton observations showed that
the X-ray spectrum was flat/concave producing a mild soft-X-ray excess,
suggesting that two components are contributing in this band (Raiteri et al.,
2009). This excess was interpreted by the helical jet model of (Villata &
Raiteri, 1999). Also BeppoSAX observations indicated presence of two
synchrotron components in the broadband SED of BL Lac (Ravasio et al., 2003).

The formation of the second emission region can be explained in the framework
of other hypotheses as well. For example, it can be a local reconnection
outflow in the “jet in a jet” scenario (Giannios et al., 2009, 2010). The
second population of energetic electrons can be formed also when the energetic
protons interact in the jet; the electrons are produced from the decay of
muons and are more energetic than the initial cooled electrons. The radiative
signature of these electrons initially appears at HE, which in time is shifted
to lower energies. In the scenario considered by Mastichiadis & Petropoulou
(2021) the second emission region could be where occasionally accelerated
protons radiate X-ray photons via the synchrotron mechanism and interact with
them via photomeson processes producing detectable high-energy neutrinos
(Stathopoulos et al., 2021).

Except for the cases when BL Lac was in a soft X-ray emission or a flaring
$\gamma$-ray state, a one-zone leptonic model considering the inverse Compton
scattering of synchrotron and BLR-reflected photons gives a reasonable
modeling of the data. The parameters values obtained from our modeling are
similar to those typically estimated for blazars; their evolution in time is
shown in Fig. 5. In some cases, the power-law index of the emitting electrons
(panel a) on Fig. 5) should be hard ($<2.0$) to account for the observed X-ray
data. Such a hard injection index is hardly possible within standard shock
accelerations, but it can be achieved when the particles are accelerated via
magnetic reconnection (Sironi & Spitkovsky, 2014). For example, Guo et al.
(2014) using fully kinetic simulations demonstrated that in highly magnetized
environments ($\sigma>>1$) the spectral index of the particles approaches
$p=1.0$; the $\sigma>>1$ condition is required so that the time scale over
which particles are injected into the acceleration region is longer than the
first-order Fermi acceleration time. The cut-off energy of the emitting
electrons takes value between $(3.01-24.57)\times 10^{2}$ and shows that the
electrons are efficiently accelerated up to $\sim 1$ GeV during the flares.
The cut-off energy is naturally formed when the acceleration is limited by the
cooling or dynamical time scale (e.g., Yan et al., 2013; Zheng et al., 2018;
Sahakyan, 2020b; Baring et al., 2017). For example, when the particles are
abruptly injected into the emitting region, they start to loose energy or
escape the region, so the HE tail of the particle distribution steepens and a
cut-off is formed. In principle, a cut-off energy within $(3.01-24.57)\times
10^{2}$ can be obtained under a reasonable assumption for the injection and
escape times. It should be noted that a similar cut-off feature in the
electron spectrum will be formed also in the case of an episodic injection
with an energy-dependent escape.

The Doppler boosting factor substantially increases when the source is bright
in the $\gamma$-ray band (Fig. 5). This is a consequence of the current
interpretation the $\gamma$-ray data as inverse Compton scattering of the
external photon field whose density transferred to the jet frame is $u_{\rm
ph}^{{}^{\prime}}\sim\delta^{2}\>u$. Therefore, any increase in the
$\gamma$-ray luminosity would require a larger $\delta$. Although the values
estimated for BL Lac and shown in Fig. 5 are not physically unrealistic, in
other interpretations, e.g., in two zone emission scenarios, a lower value of
$\delta$ would be acceptable.

The modeling allows us to assess the luminosity of the jet in various periods.
The evolution of the power (luminosity) carried by the jet in the form of
electrons and magnetic field computed as $L_{e}=\pi cR_{b}^{2}\Gamma^{2}U_{e}$
and $L_{B}=\pi cR_{b}^{2}\Gamma^{2}U_{B}$ respectively, is shown in Fig. 5
panels f) and g), respectively. Both are relatively constant with a mean
luminosity of $L_{\rm e,mean}=1.13\times 10^{44}\>{\rm erg\>s^{-1}}$ and
$L_{\rm B,mean}=2.81\times 10^{44}\>{\rm erg\>s^{-1}}$, but they slightly
increased during flaring periods. The jet is magnetically dominated in the
periods when the optical/UV data (defined by the synchrotron component)
exceeds the X-ray data (defined by SSC), for the other cases $L_{\rm e}/L_{\rm
B}>1$ was estimated. The total luminosity ($L_{\rm tot}=L_{\rm e}+L_{\rm B}$)
varies within $(0.07-5.55)\times 10^{45}\>{\rm erg\>s^{-1}}$ which is lower
than the Eddington luminosity $4.75\times 10^{46}\>{\rm erg\>s^{-1}}$ for a
black hole mass of $3.8\times 10^{8}\>{\rm M_{\rm BH}}$ (Wu et al., 2009).
Moreover, this condition will be still satisfied when considering the second
emission region with $L_{\rm tot}=(0.1-7.3)\times 10^{44}\>{\rm erg\>s^{-1}}$,
comparable to the luminosity of the other region.

## Acknowledgements

We acknowledge the use of data, analysis tools and services from the Open
Universe platform, the ASI Space Science Data Center (SSDC), the Astrophysics
Science Archive Research Center (HEASARC), the Fermi Science Tools, the All-
Sky Automated Survey for Supernovae (ASAS-SN), the Astrophysics Data System
(ADS), and the National Extra-galactic Database (NED).

This work was supported by the Science Committee of RA, in the frames of the
research project No 20TTCG-1C015.

This work used resources from the ASNET cloud and the EGI infrastructure with
the dedicated support of CESGA (Spain).

## Data availability

All data used in this paper is public and is available from the Swift, Fermi
and NuSTAR archives, from the Open Universe tools and on-line services or from
the ASAS-SN web page. The SED data used for model fitting can be shared on
reasonable request to the corresponding author.

## References

  * Abdo et al. (2010) Abdo A. A., et al., 2010, ApJ, 716, 30
  * Abdo et al. (2011) Abdo A. A., et al., 2011, ApJ, 730, 101
  * Abeysekara et al. (2018) Abeysekara A. U., et al., 2018, ApJ, 856, 95
  * Ackermann et al. (2011) Ackermann M., et al., 2011, ApJ, 743, 171
  * Ackermann et al. (2017) Ackermann M., et al., 2017, ApJ, 837, L5
  * Ajello et al. (2020) Ajello M., et al., 2020, ApJ, 892, 105
  * Albert et al. (2007) Albert J., et al., 2007, ApJ, 666, L17
  * Ansoldi et al. (2018) Ansoldi S., et al., 2018, ApJ, 863, L10
  * Arlen et al. (2013) Arlen T., et al., 2013, ApJ, 762, 92
  * Atwood et al. (2009) Atwood W. B., et al., 2009, ApJ, 697, 1071
  * Baghmanyan et al. (2017) Baghmanyan V., Gasparyan S., Sahakyan N., 2017, ApJ, 848, 111
  * Baring et al. (2017) Baring M. G., Böttcher M., Summerlin E. J., 2017, MNRAS, 464, 4875
  * Biteau et al. (2020) Biteau J., et al., 2020, Nature Astronomy, 4, 124
  * Błażejowski et al. (2000) Błażejowski M., Sikora M., Moderski R., Madejski G. M., 2000, ApJ, 545, 107
  * Bloom & Marscher (1996) Bloom S. D., Marscher A. P., 1996, ApJ, 461, 657
  * Bloom et al. (1997) Bloom S. D., et al., 1997, ApJ, 490, L145
  * Böttcher & Bloom (2000) Böttcher M., Bloom S. D., 2000, AJ, 119, 469
  * Böttcher et al. (2013) Böttcher M., Reimer A., Sweeney K., Prakash A., 2013, ApJ, 768, 54
  * Britto et al. (2016) Britto R. J., Bottacini E., Lott B., Razzaque S., Buson S., 2016, ApJ, 830, 162
  * Bulgarelli et al. (2019) Bulgarelli A., et al., 2019, A&A, 627, A13
  * Capetti et al. (2010) Capetti A., Raiteri C. M., Buttiglione S., 2010, A&A, 516, A59
  * Cash (1979) Cash W., 1979, ApJ, 228, 939
  * Catanese et al. (1997) Catanese M., et al., 1997, ApJ, 480, 562
  * Cerruti et al. (2019) Cerruti M., Zech A., Boisson C., Emery G., Inoue S., Lenain J.-P., 2019, MNRAS, 483, L12
  * Chang et al. (2020) Chang Y. L., Brandt C. H., Giommi P., 2020, Astronomy and Computing, 30, 100350
  * Corbett et al. (1996) Corbett E. A., Robinson A., Axon D. J., Hough J. H., Jeffries R. D., Thurston M. R., Young S., 1996, MNRAS, 281, 737
  * Costamante et al. (2001) Costamante L., et al., 2001, A&A, 371, 512
  * Costamante et al. (2018) Costamante L., Bonnoli G., Tavecchio F., Ghisellini G., Tagliaferri G., Khangulyan D., 2018, MNRAS, 477, 4257
  * Cutini (2012) Cutini S., 2012, The Astronomer’s Telegram, 4028, 1
  * Cutini (2021) Cutini S., 2021, The Astronomer’s Telegram, 14330, 1
  * D’Ammando (2021a) D’Ammando F., 2021a, arXiv e-prints, p. arXiv:2107.05660
  * D’Ammando (2021b) D’Ammando F., 2021b, The Astronomer’s Telegram, 14342, 1
  * D’Ammando (2021c) D’Ammando F., 2021c, The Astronomer’s Telegram, 14350, 1
  * Donea & Protheroe (2003) Donea A.-C., Protheroe R. J., 2003, Astroparticle Physics, 18, 377
  * Gao et al. (2019) Gao S., Fedynitch A., Winter W., Pohl M., 2019, Nature Astronomy, 3, 88
  * Gasparyan et al. (2018) Gasparyan S., Sahakyan N., Baghmanyan V., Zargaryan D., 2018, ApJ, 863, 114
  * Gehrels et al. (2004) Gehrels N., et al., 2004, ApJ, 611, 1005
  * Ghisellini & Tavecchio (2009) Ghisellini G., Tavecchio F., 2009, MNRAS, 397, 985
  * Ghisellini et al. (1985) Ghisellini G., Maraschi L., Treves A., 1985, A&A, 146, 204
  * Ghisellini et al. (2014) Ghisellini G., Tavecchio F., Maraschi L., Celotti A., Sbarrato T., 2014, Nature, 515, 376
  * Giannios et al. (2009) Giannios D., Uzdensky D. A., Begelman M. C., 2009, MNRAS, 395, L29
  * Giannios et al. (2010) Giannios D., Uzdensky D. A., Begelman M. C., 2010, MNRAS, 402, 1649
  * Giommi et al. (1990) Giommi P., Barr P., Garilli B., Maccagni D., Pollock A. M. T., 1990, ApJ, 356, 432
  * Giommi et al. (1999) Giommi P., Menna M. T., Padovani P., 1999, MNRAS, 310, 465
  * Giommi et al. (2005) Giommi P., Piranomonte S., Perri M., Padovani P., 2005, A&A, 434, 385
  * Giommi et al. (2021) Giommi P., et al., 2021, arXiv e-prints, p. arXiv:2108.07255
  * Guo et al. (2014) Guo F., Li H., Daughton W., Liu Y.-H., 2014, Phys. Rev. Lett., 113, 155005
  * Hosokawa et al. (2021) Hosokawa R., et al., 2021, The Astronomer’s Telegram, 14334, 1
  * IceCube Collaboration et al. (2018a) IceCube Collaboration et al., 2018a, Science, 361, 147
  * IceCube Collaboration et al. (2018b) IceCube Collaboration et al., 2018b, Science, 361, eaat1378
  * Kapanadze et al. (2018) Kapanadze B., Vercellone S., Romano P., Hughes P., Aller M., Aller H., Kapanadze S., Tabagari L., 2018, MNRAS, 480, 407
  * Keivani et al. (2018) Keivani A., et al., 2018, ApJ, 864, 84
  * Kochanek et al. (2017) Kochanek C. S., et al., 2017, PASP, 129, 104502
  * Larionov et al. (2010) Larionov V. M., Villata M., Raiteri C. M., 2010, A&A, 510, A93
  * Lott et al. (2012) Lott B., Escande L., Larsson S., Ballet J., 2012, A&A, 544, A6
  * MAGIC Collaboration et al. (2019) MAGIC Collaboration et al., 2019, A&A, 623, A175
  * Madejski et al. (1999) Madejski G. M., Sikora M., Jaffe T., BŁażejowski M., Jahoda K., Moderski R., 1999, ApJ, 521, 145
  * Mannheim (1993) Mannheim K., 1993, A&A, 269, 67
  * Mannheim & Biermann (1989) Mannheim K., Biermann P. L., 1989, A&A, 221, 211
  * Maraschi et al. (1992) Maraschi L., Ghisellini G., Celotti A., 1992, ApJ, 397, L5
  * Marchini et al. (2021a) Marchini A., et al., 2021a, The Astronomer’s Telegram, 14318, 1
  * Marchini et al. (2021b) Marchini A., et al., 2021b, The Astronomer’s Telegram, 14328, 1
  * Marscher et al. (2008) Marscher A. P., et al., 2008, Nature, 452, 966
  * Massaro et al. (2006) Massaro E., Tramacere A., Perri M., Giommi P., Tosti G., 2006, A&A, 448, 861
  * Mastichiadis & Petropoulou (2021) Mastichiadis A., Petropoulou M., 2021, ApJ, 906, 131
  * Mereu (2020) Mereu I., 2020, The Astronomer’s Telegram, 14072, 1
  * Mücke & Protheroe (2001) Mücke A., Protheroe R. J., 2001, Astroparticle Physics, 15, 121
  * Mücke et al. (2003) Mücke A., Protheroe R. J., Engel R., Rachen J. P., Stanev T., 2003, Astroparticle Physics, 18, 593
  * Murase et al. (2018) Murase K., Oikonomou F., Petropoulou M., 2018, ApJ, 865, 124
  * Neshpor et al. (2001) Neshpor Y. I., Chalenko N. N., Stepanian A. A., Kalekin O. R., Jogolev N. A., Fomin V. P., Shitov V. G., 2001, Astronomy Reports, 45, 249
  * Nilsson et al. (2018) Nilsson K., et al., 2018, A&A, 620, A185
  * Padovani & Giommi (1995) Padovani P., Giommi P., 1995, ApJ, 444, 567
  * Padovani et al. (2018) Padovani P., Giommi P., Resconi E., Glauch T., Arsioli B., Sahakyan N., Huber M., 2018, MNRAS, 480, 192
  * Pian et al. (1998) Pian E., et al., 1998, ApJ, 492, L17
  * Poole et al. (2008) Poole T. S., et al., 2008, MNRAS, 383, 627
  * Prince (2021) Prince R., 2021, arXiv e-prints, p. arXiv:2105.00221
  * Prince et al. (2021) Prince R., Raman G., Khatoon R., Agarwal A., Varun Gupta N., Czerny B., Majumdar P., 2021, arXiv e-prints, p. arXiv:2105.03937
  * Raiteri et al. (2009) Raiteri C. M., et al., 2009, A&A, 507, 769
  * Raiteri et al. (2013) Raiteri C. M., et al., 2013, MNRAS, 436, 1530
  * Rani et al. (2013) Rani B., Lott B., Krichbaum T. P., Fuhrmann L., Zensus J. A., 2013, A&A, 557, A71
  * Ravasio et al. (2002) Ravasio M., et al., 2002, A&A, 383, 763
  * Ravasio et al. (2003) Ravasio M., Tagliaferri G., Ghisellini G., Tavecchio F., Böttcher M., Sikora M., 2003, A&A, 408, 479
  * Righi et al. (2019) Righi C., Tavecchio F., Pacciani L., 2019, MNRAS, 484, 2067
  * Sahakyan (2018) Sahakyan N., 2018, ApJ, 866, 109
  * Sahakyan (2019) Sahakyan N., 2019, A&A, 622, A144
  * Sahakyan (2020a) Sahakyan N., 2020a, MNRAS, 496, 5518
  * Sahakyan (2020b) Sahakyan N., 2020b, A&A, 635, A25
  * Sahakyan (2021) Sahakyan N., 2021, MNRAS, 504, 5074
  * Sahakyan & Giommi (2021) Sahakyan N., Giommi P., 2021, MNRAS, 502, 836
  * Sahakyan et al. (2018) Sahakyan N., Baghmanyan V., Zargaryan D., 2018, A&A, 614, A6
  * Sahakyan et al. (2020) Sahakyan N., Israyelyan D., Harutyunyan G., Khachatryan M., Gasparyan S., 2020, MNRAS, 498, 2594
  * Sambruna et al. (1999) Sambruna R. M., Ghisellini G., Hooper E., Kollgaard R. I., Pesce J. E., Urry C. M., 1999, ApJ, 515, 140
  * Scargle et al. (2013) Scargle J. D., Norris J. P., Jackson B., Chiang J., 2013, ApJ, 764, 167
  * Sikora et al. (1994) Sikora M., Begelman M. C., Rees M. J., 1994, ApJ, 421, 153
  * Sironi & Spitkovsky (2014) Sironi L., Spitkovsky A., 2014, ApJ, 783, L21
  * Stathopoulos et al. (2021) Stathopoulos S. I., Petropoulou M., Giommi P., Vasilopoulos G., Padovani P., Mastichiadis A., 2021, arXiv e-prints, p. arXiv:2107.14632
  * Tavecchio et al. (2011) Tavecchio F., Becerra-Gonzalez J., Ghisellini G., Stamerra A., Bonnoli G., Foschini L., Maraschi L., 2011, A&A, 534, A86
  * Tramacere (2020) Tramacere A., 2020, JetSeT: Numerical modeling and SED fitting tool for relativistic jets (ascl:2009.001)
  * Tramacere et al. (2009) Tramacere A., Giommi P., Perri M., Verrecchia F., Tosti G., 2009, A&A, 501, 879
  * Tramacere et al. (2011) Tramacere A., Massaro E., Taylor A. M., 2011, ApJ, 739, 66
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, PASP, 107, 803
  * Villata & Raiteri (1999) Villata M., Raiteri C. M., 1999, A&A, 347, 30
  * Weaver et al. (2020) Weaver Z. R., et al., 2020, ApJ, 900, 137
  * Wehrle et al. (2016) Wehrle A. E., et al., 2016, ApJ, 816, 53
  * Wu et al. (2009) Wu Z.-Z., Gu M.-F., Jiang D.-R., 2009, Research in Astronomy and Astrophysics, 9, 168
  * Yan et al. (2013) Yan D., Zhang L., Yuan Q., Fan Z., Zeng H., 2013, ApJ, 765, 122
  * Zargaryan et al. (2017) Zargaryan D., Gasparyan S., Baghmanyan V., Sahakyan N., 2017, A&A, 608, A37
  * Zhang et al. (1999) Zhang Y. H., et al., 1999, ApJ, 527, 719
  * Zheng et al. (2018) Zheng Y. G., Long G. B., Yang C. Y., Bai J. M., 2018, PASP, 130, 083001