# Catching profound optical flares in blazars

Gopal Bhatta1, Staszek Zola2, M. Drozdz3, Daniel Reichart4, Joshua Haislip4,
Vladimir Kouprianov4, Katsura Matsumoto5, Eda Sonbas6,7, D. Caton8, Urszula
Pajdosz-Śmierciak2, A. Simon9, J. Provencal10,11, Dariusz Góra1, and Grzegorz
Stachowski2,
1Institute of Nuclear Physics Polish Academy of Sciences, PL-31342 Kraków,
Poland
2Astronomical Observatory of the Jagiellonian University, ul. Orla 171, 30-244
Kraków, Poland
3Mt. Suhora Observatory, Pedagogical University, ul. Podchorazych 2, 30-084
Kraków, Poland
4Dept. of Physics and Astronomy, University of North Carolina at Chapel Hill,
Chapel Hill, NC 27599, USA
5Astronomical Institute, Osaka Kyoiku University, 4-698 Asahigaoka, Kashiwara,
Osaka 582-8582, Japan
6University of Adiyaman, Department of Physics, 02040 Adiyaman, Turkey
7Astrophysics Application and Research Center, Adiyaman University, Adiyaman
02040, Turkey
8Dark Sky Observatory, Dept. of Physics and Astronomy, Appalachian State
University, Boone, NC 28608, USA
9Astronomy and Space Physics Department, Taras Shevshenko National University
of Kyiv, Volodymyrska str. 60, 01033 Kyiv, Ukraine
10University of Delaware, Department of Physics and Astronomy Newark, DE
19716, USA
11Delaware Asteroseismic Research Center, Mt. Cuba Observatory, Greenville, DE
19807, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Flaring episodes in blazars represent one of the most violent processes
observed in extra-galactic objects. Studies of such events shed light on the
energetics of the physical processes occurring in the innermost regions of
blazars, which cannot otherwise be resolved by any current instruments. In
this work, we present some of the largest and most rapid flares captured in
the optical band in the blazars 3C 279, OJ 49, S4 0954+658, TXS 1156+295 and
PG 1553+113. The source flux was observed to increase by nearly ten times
within a timescale of a few weeks. We applied several methods of time series
analysis and symmetry analysis. Moreover, we also performed searches for
periodicity in the light curves of 3C 279, OJ 49 and PG 1553+113 using the
Lomb-Scargle method and found plausible indications of quasi-periodic
oscillations (QPOs). In particular, the 33- and 22-day periods found in 3C
279, i.e. a 3:2 ratio, are intriguing. These violent events might originate
from magnetohydrodynamical instabilities near the base of the jets, triggered
by processes modulated by the magnetic field of the accretion disc. We present
a qualitative treatment as the possible explanation for the observed large
amplitude flux changes in both the source-intrinsic and source-extrinsic
scenarios.

###### keywords:

radiation mechanisms: non-thermal, optical — galaxies: active — blazars: jets
— method: time series analysis

††pubyear: 2015

##  1 Introduction

Blazars are a sub-class of radio-loud active galactic nuclei (AGN) featuring
relativistic jets which are closely aligned to the line of sight (Urry &
Padovani, 1995). The blazar continuum emission is non-thermal in nature, and
it is Doppler boosted and highly variable over a wide range of spatial and
temporal frequencies. Blazars consist of two kinds of sources: flat-spectrum
radio quasars (FSRQ) and BL Lacertae (BL Lac) objects. Of the two types of
sources, FSRQs show broad emission lines, while BL Lacs exhibit either weak
emission lines or their absence over the continuum. Nonetheless, the objects
are visible in the TeV energy range and constitute the dominant population of
discrete gamma-ray sources in the sky. Historically, FSRQs are considered to
be more luminous than BL Lacs, however new reports suggest that some of the BL
Lacs can be more luminous than FSRQs (see Sheng et al., 2022). The broadband
non-thermal spectrum of blazars, which extends from radio to the highest
energy $\gamma$-rays such as TeV emission, exhibits two distinct low- and
high-energy components which respectively peak between radio and soft X-rays
and between X-rays and $\gamma$-rays. The low-energy component is well
explained in terms of synchrotron emission by relativistic plasma in the
magnetized jets. However, various models, mainly leptonic or hadronic
scenarios, have been put forward to explain the origin of high-energy
emission. According to the leptonic scenario, ultra-relativistic electrons up-
scatter low-energy seed photons into X-rays and $\gamma$-rays via the inverse
Compton mechanism. Further, in the synchrotron self-Compton (SSC; e.g.
Marscher & Gear 1985; Maraschi, Ghisellini, & Celotti 1992) model, the
synchrotron photons produced by the electrons constitute a radiation field
which is inverse-Compton scattered by the co-spatial leptons. In the external
Compton (EC) models, AGN components such as the accretion disc, broad-line
region and dusty torus (see Ghisellini & Madau, 1996; Sikora, Begelman, &
Rees, 1994; Dermer, Schlickeiser, & Mastichiadis, 1992) could contribute the
low-energy seed photons required for the inverse-Compton scattering.
Conversely, the hadronic models suggest that the protons may be accelerated to
very high energies, which then produce the high-energy spectral component via
direct proton-synchrotron and/or photon-initiated cascades (M’́ucke et al.,
2003; Aharonian, 2000; Mannheim, 1993). Although both leptonic and hadronic
models can account for the origin of high-energy emission from blazars, the
hadronic models require a stronger magnetic field to radiatively cool more
massive protons. Similarly, the proton-synchrotron models also require the
acceleration of ultra high energy cosmic rays (UHECRs) in the blazar jets,
which trigger photo-pionic interactions resulting in gamma-ray emission along
with some secondary particles. In such interactions, neutrino production in
the jets becomes a natural outcome. Indeed, the IceCube (IceCube Collaboration
et al., 2018) experiment detected high energy neutrinos potentially associated
with gamma-ray flaring in the blazar TXS 0506+056.

Blazars are characterized by multi-wavelength (MWL) flux variability over
diverse timescales (e.g. see Bhatta & Webb, 2018; Bhatta et al., 2016b;
Bhatta, Mohorian, & Bilinsky, 2018; Bhatta & Dhital, 2020; Bhatta, 2021).
Although the statistical variability properties can be largely represented by
a single power-law spectral density (see optical; Nilsson et al. 2018; and
$\gamma$-ray; Bhatta & Dhital 2020), blazars often display complex variability
patterns, such as red-noise like variability superimposed by occasional sharp
rises in the flux. In majority of the cases, such rises in flux possess a
well-defined shape and last for a definite duration and therefore can be
identified as distinct flaring events, which can either be observed
simultaneously in MWL observations (e.g. Acciari et al., 2020; Abeysekara et
al., 2018; Baloković et al., 2016; Hayashida et al., 2015; Aleksić et al.,
2015), or only in a specific wavelength band, often termed as “orphan” flares,
(e.g. see Chatterjee et al., 2013, for optical orphan flares). Some of the
dominant flaring events, which show a large rise in the flux, typically last
from a few weeks to months. These events suggest the presence of extreme
physical conditions prevalent around the central engine as well as in the
jets, which may drive the most efficient particle acceleration and cooling
processes. In the gamma-ray domain, blazar flaring episodes are often
accompanied by a large swing in the rotation of the plane of optical
polarization (Blinov et al., 2018), consistent with violent collision between
the relativistic shock waves and stationary structures such as Mach disc.
Moreover, such events are also seen to be accompanied by ejection of radio
knots seen in Very Long Baseline Array images (Park et al., 2019; Marscher et
al., 2010). Studies of blazar radio jets seem to indicate that, in most cases,
flaring events in blazars can be linked with disturbances propagating along
the jet which lead to the ejection of radio knots showing apparent
superluminal motion(Jorstad & Marscher, 2016; Kellermann et al., 2004; Jorstad
et al., 2001). In the case of weakly magnetized jets, the shock waves
propagating along the jet can energize the particles to induce the flaring
events (Lind & Blandford, 1985). On the other hand, if the jets are highly
magnetized, magnetic reconnection might play a dominant role in particle
acceleration(Giannios, 2013; Nalewajko et al., 2011) before energy
dissipation. Also the re-collimation of the shock waves can lead to the
formation of rapid flares (Bromberg & Levinson, 2009)

Table 1: General information about the sample of blazar targets  Source name | Source class | R.A. (J2000) | Dec. (J2000) | Redshift (z)
---|---|---|---|---
OJ 49 | BL Lac, LSP | $08^{h}31^{m}48.88^{s}$ | $+04^{d}29^{m}39.086^{s}$ | 0.17386
S4 0954+658 | BL Lac, LSP | $09^{h}58^{m}47.2^{s}$ | $+65^{d}33^{m}55^{s}$ | 0.368
TXS 1156+295 | FSRQ, LSP | $11^{h}59^{m}032.07^{s}$ | $+29^{d}14^{m}42.0^{s}$ | 0.729
3C 279 | FSRQ, LSP | $12^{h}56^{m}11.1665^{s}$ | $-05^{d}47^{m}21.523^{s}$ | 0.536
PG 1553+113 | BL Lac, HSP | $15^{h}55^{m}43.044^{s}$ | $+11^{d}11^{m}24.365^{s}$ | 0.36

In this work, we present the results of our analysis of flaring observations
of five blazars obtained by our group in the course of long-term optical
monitoring of AGN. In Section 2, the observations of the source sample and the
relevant optical data processing are described. A brief description of each of
the sample sources is presented in Section 3. In Section 4, several analytical
approaches using various methods including fractional variability, flux
distribution, PSD and QPOs are introduced, and the results of the analyses of
the light curves are presented. The results along with their possible
implications are discussed in Section 5, and we summarize our conclusions in
Section 6.

##  2 Observations and data processing

Long term monitoring of a sample of quasars was primarily carried out using
small telescopes operated by the Skynet Robotic Telescope Network (Zola et
al., 2021). Additional data were collected with the 60 cm telescope located at
Adiyaman University Astrophysics Application and Research Center, Turkey, the
50 cm telescope at the Osaka University Observatory in Japan, and with two
telescopes in Poland: a 60 cm at the Mt. Suhora Observatory and a 50 cm at the
Jagiellonian University Observatory in Krakow. The wide band R filter (Bessell
prescription) was most often used for the monitoring. Longer runs were also
carried out, mostly at the Krakow and Mt. Suhora sites.

The data taken by Skynet consists of several scientific images of a target
taken each night, subsequently reduced for bias, dark and flatfield by the
network pipeline. Other sites provided raw images accompanied by calibration
frames. We performed reduction of raw images with the standard procedure:
calibration for bias, dark and flatfield (usually taken on the sky) with the
IRAF package, while extraction of magnitudes was done using aperture
photometry with the CMunipack program, which implements the DAOPHOT algorithm.
As a result, differential magnitudes were derived with comparison stars for
each object chosen so as to be visible in the field of view of all telescopes.
Their constancy was verified with check stars which were similarly chosen. The
sample of sources along with their classes, positions and red-shifts are
presented in Table 1, which also lists the blazar source classification based
on synchrotron peak frequency, that is, high synchrotron peaked blazars (HSP;
$\nu_{\mathrm{peak}}^{S}>10^{15}$ Hz), intermediate synchrotron peaked blazars
(ISP; $10^{14}<\nu_{\mathrm{peak}}^{S}<10^{15}$ Hz), and low synchrotron
peaked blazars (LSP; $\nu_{\mathrm{peak}}^{S}<10^{14}$Hz) (Abdo et al., 2010).
Similarly, for a given source, the total observation duration, the number of
observations and the mean magnitude are listed in the 2nd, 3rd and 4th
columns, respectively, of Table 2.

Figure 1: Optical (R-band) observations of the sample of blazars

##  3 Source description

### OJ 49

Blazar OJ 49 is a BL Lacertae object located at a red-shift of 0.17386. The
Very Large Array (VLA) radio image of the source at 20 cm (Antonucci, &
Ulvestad, 1985) shows a sharply curved extending jet, while a 43 GHz Very Long
Baseline Array (VLBA) image shows a highly polarized jet extending about 0.6
mas from the core (Lister et al., 1998). Similarly, VLBA images at 22 GHz
revealed a prominent jet ejecting knotty components at apparently superluminal
speeds (Jorstad et al., 2001). During VLBI Space Observatory Program mission,
the core size of the 5 GHz radio emission was estimated to be 0.5 mas (Dodson
et al., 2008). In the optical and near-IR bands, the blazar exhibits strong
intraday variability of polarization and total flux (see Sitko et al., 1985;
Smith et al., 1987).

### S4 0954+658

BL Lac S4 0954+658 is known to display significant flux and polarization
variability on both intra-day and longer timescales. Long term optical
variability of the source by Raiteri et al. (1999) showed large amplitude-flux
modulations. Similarly, Papadakis et al. (2004) made multi-band optical
observations of the source for a few nights and found the source to be
variable by $<$ 5% within the observational period. Morozova et al. (2014)
studied the multi-wavelength behavior of S4 0954+658 during a powerful
outburst in March-April 2011 using optical (R-band) photometric and
polarimetric monitoring and VLBA observations. The authors reported an
increase of the flux by 2.8 mag within a period of two months, and a steep
intra-night flux increase of $\sim$0.7 mag. Similarly, in a study involving
multi-color photometric and polarization observations of the blazar during
2008–2012, the source revealed a power-law spectrum and high degree of
polarization, confirming the synchrotron nature of the emission (Hagen-Thorn
et al., 2015). The source also displayed significant color variability with a
trend of bluer-when-brighter (BWB). Morozova et al. (2016) investigated the
behavior of the source during the outburst in early 2015 using observations
from cm-wave to $\gamma$-ray energies. It was found that the optical flaring
coincided with a similar flare observed in $\gamma$-rays. During the same
period, very-high-energy $\gamma$-ray emission was detected and the ejection
of a new, bright polarized superluminal knot was observed by the VLBA at 43
GHz. More recently, Vlasyuk et al. (2022) reported the fastest fall in the
source’s optical flux, by $\sim$ 0.25 mag. within 15 minutes, accompanied by a
minute-timescale QPO.

### TXS 1156+295

Blazar TXS 1156+295, also known as Ton 599 and as 4FGL J1159.5 + 2914 in the
Fermi-LAT 4th catalog, is an FSRQ situated at the position of R.A. =
$11^{h}59^{m}31.8^{s}$ and Dec. = +29${}^{o}14^{\prime}43^{\prime\prime}.8$
and lies at a redshift $z=0.725$ (Hewett & Wild, 2010). The source was first
detected in the $\gamma$-ray band by the Energetic Gamma Ray Experiment
Telescope (EGRET), and later was also detected in very high-energy emission (
$>$ 100 GeV) by VERITAS (Mukherjee & VERITAS Collaboration, 2017). In the
optical band, the blazar is highly variable in all timescales. Fan et al.
(2006) presented a study of the sources using photometric observations which
showed a large variation ($\Delta m\sim$ 5.8 mag) in the optical flux on
timescales of a few years. In 2017 the source was reported to have undergone
optical flaring (Pursimo et al., 2017). In the $\gamma$-ray band, power
spectral density and flux distribution analysis of decade-long Fermi/LAT
observations was carried out by (Bhatta & Dhital, 2020). More recently, Rajput
& Pandey (2021) studied the flux and spectral variability of the source during
its $\gamma$-ray flaring in 2021. Also, variability of emission-line during a
non-thermal outburst was reported by Hallum et al. (2022).

### 3C 279

Blazar 3C 279 is a FSRQ source profusely emitting in hard X-ray and
$\gamma$-rays. Highly variable across a wide range of spectral bands (see
Hayashida et al., 2015; Paliya et al., 2016, and the references therein), it
is one of the few FSRQs detected above 100 GeV (MAGIC Collaboration et al.,
2008). The source reveals a compact, milliarcsecond-scale radio core ejecting
radio knots with a bulk Lorentz factor $\Gamma=15.5\pm 2.5,$ in a direction
making an angle $\theta_{obs}=2.1\pm 1.1^{\circ}$ to the line of sight
(Jorstad et al., 2005, 2004). During Whole Earth Blazar Telescope (WEBT)
campaigns, the source flux in the optical was reported to undergo an
exponential-type decay on a timescale of $\sim 10$ d (B’́ottcher et al.,
2007). Similarly, Larionov et al. (2008) in another WEBT campaign in 2006-2007
observed a slower but large flux decline, $\sim 3$ mag on a timescale of
$\sim$100 d, in the optical and near-IR band. Bhatta, Mohorian, & Bilinsky
(2018) discussed the flux and spectral variability properties of the source in
the hard X-ray band during intra-day timescales. More recently, Agarwal et al.
(2019) presented multi-band optical variability of the source and found that
large amplitude variability within a timescale of a few months and a mild BWB
trend on shorter time-scales. The optical observations of this source are
presented in the bottom panels of Figure 1. A sharp fall in the flux ($\Delta
m\sim 2.0$) within a timescale of 100 days is separately shown in the right
panel.

### PG 1553+113

Blazar PG 1553+113 is a BL Lac source which has been studied from radio to
$\gamma$-rays in different observation campaigns (e.g. Osterman et al., 2006;
Ackermann et al., 2015; Raiteri et al., 2015, 2017). The source is famous for
its $\sim$ 2.18 year periodicity, first revealed in Fermi/LAT observations
(Ackermann et al., 2015). It was classified as a BL Lac object based on its
featureless spectrum (Falomo & Treves 1990) and was further sub-classified as
a high-peaked BL Lac (HBL) object (see Beckmann et al., 2002, and references
therein). Evidence of very-high-energy $\gamma$-ray emission from this source
was first reported by H.E.S.S. in 2005 (Aharonian et al., 2006) and was later
confirmed by observations above 200 GeV with the MAGIC telescope at a
significance level of 8.8 $\sigma$ (Albert et al., 2007). Due to its
featureless optical spectrum, the redshift of PG 1553+113 remains highly
uncertain. Measurements using the Cosmic Origins Spectrograph onboard the
Hubble Space Telescope yielded a lower limit of 0.395 (Danforth et al., 2010).
The small statistical uncertainties of the VERITAS energy spectrum help
constrain the upper limit of the redshift to a value of $\sim$ 0.62 (Aliu et
al., 2015). During multi-frequency WEBT campaign in 2013 April–August, the
source was found to display a general BWB trend in the optical regime (Raiteri
et al., 2015). Meng et al. (2018) studied the source using the optical multi-
band observations in the yearly timescale. During the observation period the
source exhibited moderately varying multi-band emission, without any inter-
band lag. Pandey et al. (2019) reported variable emission in the optical V and
R bands with a mean optical spectral index of $\sim$ 0.83$\pm$ 0.21.

Figure 2: Fitting of the flares observed in optical (R-band) observations of
the sample blazars with a curve (shown in blue color) parameterized with an
exponential rise and decay are presented in the upper panels of the figures.
The lower panels show the ratio between the observations and the models. The
corresponding rise and decay times from the best-fitting models are presented
in the 8th and 9th columns of Table 2.

##  4 Analysis and Results

In order to characterize the properties of the flaring episodes in the target
sources, we analyzed the optical observations through multiple analysis
methods. The methods and results of the analyses for the individual sources
are presented below.

###  4.1 Variability measures and timescales

In order to obtain a quantified measure of the observed variability in the
sources, the following three measures of variability are estimated.
Variability amplitude (VA) as given in Heidt & Wagner (1996) provides an
measure of the net flux change during the observation period and is written as

$\rm VA=\sqrt{\left(A_{max}-A_{min}\right)^{2}-2\sigma^{2}}$ (1)

where $A_{max}$ and $A_{min}$ are the maximum and minimum of the source
magnitude, respectively, and $\sigma$ represents mean error in the magnitude
measurements. The VA of the sample sources are listed in Table 2.
Subsequently, the VA in magnitudes then can be directly converted into flux
ratio using the relation $f_{1}/f_{2}=10^{-0.4(m_{1}-m_{2})}$, where $m_{1}$
and $m_{2}$ are the initial and final magnitudes of the variable source. This
allows us to compute the factor by which the flux changed during the
observational period.

We observed high-amplitude rapid variability in the blazars included in our
study. As seen in Figure 1, and also indicated by the VA (in magnitudes)
listed in the 5th column of Table 2, the optical flux of blazar 3C 279 was
observed to change by $\sim 13$ times within a timescale of a few hundred
days. Similarly, in OJ 49 we measured a flux change by a factor of $\sim$12
within the period of $\sim 12$ days as shown in the top panel of Figure 1, and
in TXS 1156+295 by a factor of $\sim$16 within a timescale of $\sim$22 days as
seen in the middle panel of Figure 1. Likewise, the optical light curve
showing flux modulation in the blazar PG 1553+113 is presented in the top
right panel of Figure 1.

The VA considers only extreme values and therefore provides a measure for the
peak-to-peak magnitude change. Average variability during the entire period
can be quantified by estimating their fractional variability (FV) given as

$F_{var}=\sqrt{\frac{S^{2}-\left\langle\sigma_{err}^{2}\right\rangle}{\left\langle
F\right\rangle^{2}}},$ (2)

where $S^{2}$ and $\left\langle\sigma_{err}^{2}\right\rangle$ represent the
variance and the mean of the squared measurement errors, respectively; and the
uncertainty in FV can be expressed as

$\centering\sigma_{F_{var}}=\sqrt{F_{var}^{2}+\sqrt{\frac{2}{N}\frac{\left\langle\sigma_{err}^{2}\right\rangle^{2}}{\left\langle
F\right\rangle^{4}}+\frac{4}{N}\frac{\left\langle\sigma_{err}^{2}\right\rangle}{\left\langle
F\right\rangle^{2}}F_{var}^{2}}}-F_{var}\@add@centering$ (3)

(Vaughan et al. 2003, see also Bhatta & Webb 2018). The FV values for the
sample sources are listed in the 6th column of Table 2, which show large large
average flux variability during the period.

The variability timescale ($\tau_{var}$) can be taken as the e-folding
timescale of flux change given by

$\tau_{var}=\left|\frac{\Delta t}{\Delta lnF}\right|,$ (4)

Burbidge et al. 1974, see also Bhatta, Mohorian, & Bilinsky 2018), where
$\Delta t$ is the time interval corresponding to the change in natural
logarithm of flux measurements. Using the above relation, $\tau_{var}$ of the
sample light curves are calculated and listed in the 7th column of Table 2. Of
the sample sources, a shortest $\tau_{var}$ of 11 minutes is observed blazar
3C 279, whereas the source PG 1553+113 appears to show relatively slower
variability with a $\tau_{var}$ of 65 minutes.

Rapid variability in shorter timescales can be associated with the synchrotron
cooling timescales, which can be expressed as,

$t_{cool}\sim 7.74\times 10^{8}\gamma^{-1}B^{-2}\ s,$ (5)

where $\gamma$ and B represent electron Lorentz factor and ambient magnetic
field, respectively. Also, we use $\beta\sim 1$ considering ultra-relativistic
electrons (see also Bhatta, Mohorian, & Bilinsky, 2018). Using a typical jet
magnetic field of 1 Gauss and a minimum variability timescale of 30 minutes,
the energy of the relativistic electrons emitting optical synchrotron emission
can be estimated to $\sim 3\times 10^{5}$ Lorentz factors. Moreover, following
the causality argument, the timescale $\tau_{var}$ can be used to estimate the
upper limit for the minimum size of the emitting region ($R$) as given by
$R\geq\frac{\delta}{\left(1+z\right)}c\tau_{var}$; where $\delta$, Doppler
factor, is defined as $\delta=(\Gamma\left(1-\beta cos\theta\right))^{-1}$,
and for the velocity $\beta=v/c$ the bulk Lorentz factor can be written as
$\Gamma=1/\sqrt{1-\beta^{2}}$. If we use a typical value of Doppler factor of
10 and $z=0.5$, a 30-minute minimum variability timescale could have arisen
from a region that is 1/1000-th of a parsec in scale.

Table 2: Optical observations of a sample of blazars and their variability properties Source name | Duration (d) | Npt. | mean mag | VA (mag) | Fvar (%) | tvar (min.) | tr (d) | td (d)
---|---|---|---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9)
OJ 49 | 164.28 | 620 | 16.28 | 2.75 | 66.23$\pm$0.24 | 38.24$\pm$11.80 | 0.57 | 0.44
TXS 1156+295 | 79.65 | 650 | 17.45 | 3.30 | 62.78$\pm$0.13 | 19.06$\pm$13.73 | - | 4.10
PG 1553+113 | 97.43 | 1071 | 15.57 | 0.91 | 23.61$\pm$0.21 | 64.76$\pm$11.18 | - | -
3C 279 | 461.58 | 1831 | 17.59 | 2.78 | 46.16$\pm$0.50 | 11.73$\pm$7.80 | - | -
S4 0954+658 | 242.70 | 2988 | 14.80 | 2.61 | 47.60$\pm$0.17 | 17.10$\pm$6.18 | 3.05 | 5.15

###  4.2 Symmetry Analysis: Rise and Decay Profiles

The rise and decay profiles of flares in blazar light curves can be associated
with the particle acceleration and cooling timescales, respectively, and
thereby can be linked to the physical processes leading to the flaring
episodes. To characterize the flaring properties, we performed symmetry
analysis of the flares. For that purpose, we took parts of the light curves of
the sources OJ 49, S4 0954+658, and TXS 1156+295 (see Figure 1) which showed
well resolved flares, with a well-defined generally monotonic rise and decay,
and the three flares were fitted by a functional form representing the
temporal structure of the exponential rise and decay described by

$F\left(t\right)=F_{c}+F_{0}\left[e^{\frac{t_{0}-t}{t_{r}}}+e^{\frac{t-t_{0}}{t_{d}}}\right]^{-1},$
(6)

where $F_{c}$ is the constant flux level, $F_{0}$ is the amplitude of the
flaring structure and $t_{0}$ the center of the flare, and $t_{r}$ and $t_{d}$
are the rise and decay times of the flares (Abdo et al., 2010). The parts of
the light curves of the sources which show distinct flares and the
corresponding functional fits are shown in Figure 2. The rise and decay
timescales determined from the fitting are listed in the 8th and 9th columns,
respectively, of Table 2. Using the obtained timescales, a symmetry parameter
written as $\xi=\left(t_{d}-t_{r}\right)/\left(t_{d}+t_{r}\right)$ can be
defined within [-1,1] such that $\xi=-1,+1$ represent completely right- and
left-asymmetric flares, whereas $\xi=0$ represents completely symmetric
flares. The asymmetry parameters for OJ 49 and S4 0954+658 $\xi=0.13$ and
$\xi=-0.25$, respectively. The result indicates that the flares in these
sources are, respectively, slightly left- and right-asymmetric. Since the
source TXS 1156+295 was only observed during its decay phase its $\xi$ could
not be estimated. Also, although the amplitude of the flares in OJ 49 and S4
0954+658 are similar, with normalized amplitudes ($F_{0}/F_{c}$)=8.5 and 8.2
respectively, the rise and decay times are considerably shorter in OJ 49,
indicating fast flaring events. In the case of TXS 1156+295 an even faster
decay of flux, by 11.3 normalized amplitudes, is observed within 4.3 days.

###  4.3  Periodicity Analysis

QPOs in blazars with characteristic timescales of a few years have been
frequently reported (Bhatta et al. 2016c, see also Zola et al. 2016). However,
not many blazar QPOs are observed on timescales of a few days or weeks. We
searched for the possible periodic flux modulations in the optical bands using
Lomb-Scargle method (Lomb, 1976; Scargle, 1982). The method modifies the
conventional discrete Fourier periodogram such that the least-square fitting
of sine waves of the form $X_{f}(t)=A\cos\omega t+B\sin\omega t$ to the data
is minimized. The periodogram is given as

$P=\frac{1}{2}\left\\{\frac{\left[\sum_{i}x_{i}\cos\omega\left(t_{i}-\tau\right)\right]^{2}}{\sum_{i}\cos^{2}\omega\left(t_{i}-\tau\right)}+\frac{\left[\sum_{i}x_{i}\sin\omega\left(t_{i}-\tau\right)\right]^{2}}{\sum_{i}\sin^{2}\omega\left(t_{i}-\tau\right)}\right\\}\,,$
(7)

where $\tau$ is given by $\tan\left(2\omega\tau\right)=\sum_{i}\sin 2\omega
t_{i}/\sum_{i}\cos 2\omega t_{i}$ .

The Lomb-Scargle periodogram (LSP) of the source light curves of PG 1553+113
and 3C 279 are presented in the right panels of Figure 3 and similarly LSP of
the source OJ 049 is presented in Figure 4. The light curve of blazar PG
1553+113 shows peaks at the timescales of $7\pm 0.8$ ,$25\pm 3$, and $59\pm 7$
d, which possibly may be the signs of quasi-periodic oscillations. Similarly
the LSP of the blazar 3C 279 shows two prominent peaks around the timescales
of 22 and 33 days. Furthermore, in the LSP of the source OJ 049 a prominent
peak around 12 day time scale is observed.

Figure 3:  Lomb-Scargle periodogram and the R band light curve of the blazars
PG 1553+113 and 3C 279 are shown in the top and bottom panels, respectively.
The 90 and 99% significance contours from simulation are shown on the LSP
diagram by magenta and red curves, respectively. Figure 4:  Lomb-Scargle
periodogram of the blazar OJ 49 showing possible $12\pm 1.5$ d periodicity.
The 90 and 99% significance contours from simulation are shown by the magenta
and red curves, respectively

To compute the significance of the observed periodogram features in the source
LSP, a large number of light curves were generated by Monte Carlo simulations
and their statistical properties were utilized. In particular, the source
periodograms based on Discrete Fourier Transform (DFT) were linearly fit in
the logarithmic frequency space to obtain a best representative model power
spectral density (PSD) using the method described in Vaughan (2005). To avoid
artifacts associated with DFT, some of the data were linearly interpolated to
create a evenly spaced sampling. Using the best-fit PSD model, 10 000 light
curves mimicking the observations in duration and sampling rate were simulated
and the distribution of these simulated LSP periodograms was used to determine
90 and 99% significance contours (for detail see Bhatta et al., 2016c; Bhatta
& Dhital, 2020), shown in magenta and red curves in the corresponding LSP
figures. It can be seen that only the 12 d period in the source OJ 049 appears
at over the 99% significance level, although the 7 d period in PG 1552+113,
and 22 and 33 d periods in the source 3C 279 are clearly observed. We note
that in another blazar Mrk 501 a QPO of a similar timescale of 23 d has been
reported previously (see Rieger & Mannheim, 2000).

##  5  Discussion

Blazars are found to be violently variable over all timescales from a few
minutes to decades. Apart from stochastic variability often represented by
power-law PSD, blazar light curves frequently show MWL flaring events
characterized by well-resolved trends of monotonic rise or decline of flux
which can last from a few weeks to months. Studies by several blazar
monitoring groups are particularly focused on flares in gamma-rays, which are
frequently associated with the ejection of radio knots visible in VLBA images
(Agudo et al., 2011; Wehrle et al., 2012; Jorstad et al., 2013; Morozova et
al., 2014). There is growing evidence that flaring events observed in the
gamma-ray wavelength range could be associated with superluminally-moving
features crossing stationary features along the jet; it has been found that a
large number of flares occur in coincidence with the passage of superluminal
knots through the millimeter-wave core(see Jorstad et al., 2017, and the
references therein).

Flares can be explained in terms of the sudden enhancement of the blazar flux
followed by a strong energy dissipation event. The particles are accelerated
to high energies via a number of particle acceleration mechanisms. In the
shock diffusive acceleration model, shock waves compress and order the
upstream magnetic field, causing the particles to accelerate through the Fermi
acceleration mechanism and subsequently leading to the formation of outbursts
of synchrotron emission (Blandford & Eichler, 1987; Hughes et al., 1998).
Similarly, in the turbulent jet models the main jet can be thought of being
divided into a large number of sub-volumes moving relativistically in random
directions (see e.g. Marscher, 2014; Narayan & Piran, 2012). In this scenario,
the passage of shock waves in the turbulent flow of relativistic plasma can
heat the particles by compressing the plasma, and subsequently depending upon
the size and direction of the motion one single turbulent cell can beam
dominantly to appear as a flare. In the case that the jet is highly
magnetized, the turbulence can trigger intermittent magnetic instabilities,
such as kink instabilities (Spruit et al., 2001) or inversions of magnetic
field near the base of the jet (Giannios & Uzdensky, 2019), which lead to
magnetic reconnection events and consequently acceleration of the particles to
high energies. (Guo et al., 2014a, b). Similarly, magnetic reconnection events
lead to the formation of plasmoids, which can grow to produce rapid flares
with distinct and resolved envelopes. In this study, the flare are found to be
slightly asymmetrical, which might suggest that the observed skewness in the
flare profiles can be linked to both disturbance passing through the emission
region and/or geometric effects such as those resulting from light crossing
time. The origin of these flares could be extrinsic and intrinsic in nature,
as discussed qualitatively in the following sections.

###  5.1  Source intrinsic scenario

The flares in blazars can be explained in the context of internal shocks
propagating along the blazar relativistic jets. As an illustrative example of
blazar flares which last about a week, here we attempt to produce a flare
using the method worked out by Kirk et al. (1998). This provides a framework
for time-dependent analysis of the homogeneous single zone leptonic model, in
which variable emission leading to large flares in the source flux is expected
owing to gradual particle acceleration at the shock front and subsequent
radiative cooling at the emission region. The model assumes that the shocks
are propagating along a cylindrical jet aligned near the line of sight, and
that particles injected at a constant rate at the shock wavefront lose energy
primarily via synchrotron emission.

Figure 5: Top: Time evolution of maximum electron Lorentz factor as it
accelerates (left panel). Cooling length behind the shock wave front as a
function of electron energy for three different times. Bottom: The normalized
integrated particle (electron) density in the acceleration zone for three
different times 5 tacc, 50 tacc and 500 tacc as shown by the black, blue and
red curves, respectively (left panel). Rise and decay of the blazar emission
due to particles accelerated at the relativistic shocks front and thereby
cooling, respectively (right panel).

The evolution of the particles in the acceleration zone is given by the
diffusion equation,

$\frac{\partial N}{\partial
t}+\frac{\partial}{\partial\gamma}\left[\left(\frac{\gamma}{t_{acc}}-\beta_{s}\gamma^{2}\right)\right]N+\frac{N}{t_{acc}}=Q\delta(\gamma-\gamma_{0}).$
(8)

where $\beta_{s}\gamma^{2}$ represents the loss of the energy by the
synchrotron radiation with

$\beta_{s}=\frac{4}{3}\frac{\sigma_{T}}{m_{e}c^{2}}\left(\frac{B^{2}}{2\mu_{0}}\right),$
(9)

where $\sigma_{T}$, B and $\mu_{0}$ are the Thompson-scattering cross-section,
the magnetic field and the permeability of free space, respectively.

Particle acceleration reaching up to $\gamma_{max}=(\beta_{s}t_{acc})^{-1}$ is
described by the equation

$\gamma
1\left(t\right)=\left(\frac{1}{\gamma_{max}}+\left[\frac{1}{\gamma_{0}}-\frac{1}{\gamma_{max}}\right]e^{-t/t_{acc}}\right)^{-1},$
(10)

which shown in the top panel of Figure 5. The distributions of the
(normalized) particle densities over the particle energies for the three
different acceleration times, i. e., 5, 50 and 500 tacc are presented in the
bottom left panel of Figure 5 (see similar Figure 1 in Kirk et al., 1998). The
breaks which naturally appear in the curves divide the population into
particles with energies greater than those which cool within the source and
low energy particles which do not cool within the emission region.

The particle enhancement in the regions is described by the equations,
$Q\left(t\right)=Q_{0}$ for t$<$0 and t$>$$t_{f}$ and
$Q\left(t\right)=\left(1+\eta_{f}\right)Q_{0}$ for 0$<$t$<$ $t_{f}$.

We have then,

$I\left(\nu,t\right)=I_{1}\left(\nu,t\right)+\eta_{f}\left[I_{1}\left(\nu,t\right)-I_{1}\left(\nu,\left(1-u_{s}/c\right)t_{f}\right)\right].$
(11)

For small angles $cos\theta\rm\sim 1$, and $\delta\sim 2\Gamma$ such that the
relations $I(\nu,t)=\delta^{3}I(\nu^{\prime},t^{\prime})\sim
8\Gamma^{3}I(\nu^{\prime},t^{\prime})$, and
$\nu=\Gamma\left(1+\beta\right)\approx 2\Gamma\nu^{\prime}$ are used to
transform the source rest-frame quantities (primed) to the observer’s frame
(unprimed). We take $\nu=4.55\times 10^{14}$ Hz, the frequency corresponding
to the mean effective wavelength of R-band (658 nm) and the shocks traveling
down the blazar jet with relativistic speeds ($\beta_{s}=0.1$). The resulting
normalized intensity profile mimicking the flaring behavior observed in the
blazars is presented in the bottom right panel of Figure 5.

###  5.2  Source extrinsic scenario

A large flare in the blazar flux can also appear owing to the Doppler boosted
emission when the emission regions travel along a curve trajectory in the
jets. In such cases, the rest frame flux (${F}^{\prime}_{{\nu^{\prime}}}$) is
related to the observed flux ($F_{\nu}$) through the equations

$\frac{F_{\nu}(\nu)}{{F}^{\prime}_{{\nu}^{\prime}}(\nu)}=\delta^{3+\alpha}\quad\text{and}\quad\delta(t)=\frac{1}{\Gamma\left(1-\beta
cos\theta\right)}$ (12)

Figure 6: Left: As the angle between the emission region and the line of sight
decreases (top panel), the flux appears to flare as a result of relativistic
beaming (bottom panel). The three curves correspond to the three different
values of the bulk Lorentz factors. Right: Similarly, flaring of flux (bottom
panel) owing to the increase in the bulk Lorentz factor (top panel). The three
curves correspond to the three different values of the angles of sight.

The optical spectral slope is, in general, larger than 1 in LSPs, (in our
case, OJ 49, S4 0954+658, TXS 1156+295, and 3C 279) and less than 1 in case of
HSPs (here, PG1553+113). For illustrative purposes, however, here the spectral
index is taken as $\alpha\sim 1$. If we assume the apparent flux rise happens
purely due to the changes in the Doppler factor ($\delta$), which further can
be related to the changes in the angle with the line of sight or/and the
changes in the bulk Lorentz factor. Here as an illustration we treat two cases
considering change in $\theta$ and $\Gamma$ separately.

  * •

Change in the angle of the line of sight:
The angle between the emission region and the line of sight is allowed to
gradually decrease, as approximated by $\rm{\theta=\theta_{0}-Asin^{2}\omega
t}$, where $\theta_{0}$=5.7 and A=3.7 as shown in the top left panel of Figure
6. The resulting flux rise profiles for the three values of the bulk Lorentz
factor, i.e. $\Gamma$=10, 15 and 20 represented by the blue, magenta and red
curves, respectively, are shown in the lower left panel of Figure 6.
Furthermore to approximate the flaring behavior that last about a week in the
observer’s frame, $\omega$ was chosen to be $15\times 10^{-6}$ rad/sec.

  * •

Change in $\Gamma$ :
The bulk Lorentz factor of the dominant emission region is allowed to
gradually increase, as approximated by $\rm{\Gamma=\Gamma_{0}+Asin^{2}\omega
t}$, with $\Gamma_{0}$=10 and A=20, so that the plasma blob traveling at a
speed of $\Gamma=10$ accelerates to attain $\Gamma=30$ and subsequently
decelerates back to the previous speed. The evolution of $\Gamma$ in time can
be seen in the top right panel of Figure 6. As in the previous case, the value
of $\omega$ was chosen to be $15\times 10^{-6}$ rad/sec such that the event
lasts about a week in the observer’s frame. The resulting flux rise profiles
for the three values of the angles of the line of sights, i.e. $\theta$=1.5,
3.5 and 6.5o are represented by the blue, magenta and red curves,
respectively, shown in the lower right panel of Figure 6.

###  5.3 Stochastic flux variability and possible periodicity

Blazars exhibit complex flux variability patterns with a possible mixture of
several components such as general stochastic fluctuations, occasional large
amplitude flares and possible quasi-periodic oscillations arising from various
instabilities both in the disc and the jet. The general aperiodic and
stochastic variability observed in multi-wavelength observations of AGN is
largely represented by red-noise (e.g. see Isobe et al., 2015; Kelly et al.,
2011; Bhatta & Dhital, 2020). The spectral power density of such variability
is most consistent with a single power-law model, with negative spectral index
ranging from $\sim$1–2. A negative power-law implies that the variability
power grows towards longer timescales, meaning the source fluxes vary by
larger amplitudes over yearly timescales compared to shorter ones, i.e.
intraday, daily, weekly and monthly timescales. However it is interesting to
note that in the particular case of the 3C 279 light curve in the bottom left
panel of Figure 3, the mean flux nearly remains stationary while the shorter
term flux fluctuates rapidly. Similar observations can be made about PG
1153+113 when considering only the light curve before the vertical dashed
line. It is only when we consider the full-length light curve that the longer
term variability appears dominant. It is possible that the processes driving
variability which can be characterized by a negative power-law index are
distinct in origin compared to the commonly-observed variability characterized
by positive power-law index. They could be signatures of the quasi-stationary
isolated events which gets mixed with the general red-noise like variability.

Apart from the general aperiodic variability, Lomb-Scargle periodogram
analysis of the light curves revealed hints of QPOs in some of the sources.
However, the estimated significance of the peaks in the periodogram is
moderate at about 90 %. Nevertheless, it is important to note that in blazar
light curves heavily dominated by red noise the actual QPO signals could
appear relatively weak (see discussion in Bhatta & Dhital, 2020). The
periodicity analysis indicates that in blazar 3C 279 the two periods, well
above 90% significance level, appear to be in a 3:2 ratio. In X-ray binaries
such QPOs with periods in a 3:2 ratio are interpreted in terms of resonant
oscillations of accretion flow (e.g. Kluźniak & Abramowicz, 2005). Blazar QPOs
on timescales of a few months are relatively rare, although yearly timescales
QPOs have been reported in several blazars . Note that although QPOs in a
timescales of a few days are more likely to originate at the accretion disc,
we detect them through the jet emission. Possible interpretations of the
blazar QPO are discussed in detail in our previous works (see Bhatta & Dhital,
2020; Bhatta, 2019, 2017, and references therein).

##  6 Conclusion

A variability study was carried out using optical observations of the blazars
3C 279, OJ 49, S4 0954+658, TXS 1156+295 and PG 1553+113 spanning several
weeks. The light curves using observations acquired through our ground based
telescope networks showed some of the extraordinary flaring events in which a
flux change of nearly 10 times was observed within a timescale of $\sim$10
days. Flares in the sources OJ 49, S4 0954+658, and TXS 1156+295 were studied
using a functional form of exponential rise and decay. Of these sources,
slightly asymmetric flares of comparable normalized amplitude of $\sim 8$ were
observed in the blazars OJ 49, S4 0954+658; whereas light curves of the source
TXS 1156+295 revealed a fasted flux decay by $\sim 11$ normalized amplitudes
within the timescale of four days. Such a rapid, large amplitude flux change
reflects the most violent processes in the jet, e.g. shock waves and magnetic
re-connection events. To explain the observed flares, qualitative descriptions
of possible intrinsic and extrinsic scenarios were presented. In the source-
intrinsic scenario, particle injection at the shock wave front can result in a
flare in the source flux; whereas in the extrinsic scenario turbulent flow in
the jet might lead to a change in the Doppler factor of a single, dominant
energized cell such that a slight change in the speed or angle to the line of
sight can yield a large observed change in flux.

Furthermore, periodicity analysis was performed on the sources 3C 279, OJ 49,
TXS 1156+295 and PG 1553+113 using Lomb-Scargle periodograms supplemented by a
large number of simulated light curves generated using Monte Carlo methods.
The result indicated that a strong spectral peak was detected at the
characteristic timescale of $\sim$ 12 days to above 99% significance level
against power-law noise. In addition, two potentially significant periods 33
and 22 day, at a ratio of 3:2, were found to be detected at a significance
level of $\sim 95\%$.

## Acknowledgments

We are grateful to the anonymous reviewer for their constructive comments,
which helped improve the quality of the paper significantly. We acknowledge
the support of the Polish National Science Centre through the grants
UMO-2017/26/D/ST9/01178 (GB), 2018/09/B/ST9/02004 (SZ), and 2020/39 / B / ST9
/ 01398 (DG). KM acknowledges JSPS KAKENHI grant number 19K03930.

## Data availability

The data used in this work can be shared on reasonable request to the
corresponding author.

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 722, 520
  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, ApJ, 716, 30. doi:10.1088/0004-637X/716/1/30
  * Abramowicz et al. (2003) Abramowicz, M. A., Bulik, T., Bursa, M., et al. 2003, A&A, 404, L21
  * Abeysekara et al. (2018) Abeysekara, A. U., Benbow, W., Bird, R., et al. 2018, ApJ, 856, 95. doi:10.3847/1538-4357/aab35c
  * Acciari et al. (2020) Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2020, ApJS, 248, 29. doi:10.3847/1538-4365/ab89b5
  * Ackermann et al. (2015) Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJL, 813, L41
  * Agarwal et al. (2019) Agarwal, A., Cellone, S. A., Andruchow, I., et al. 2019, MNRAS, 488, 4093
  * Agudo et al. (2011) Agudo, I., Marscher, A. P., Jorstad, S. G., et al. 2011, ApJ, 735, L10. doi:10.1088/2041-8205/735/1/L10
  * Aharonian et al. (2007) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, ApJL, 664, L71z
  * Aharonian et al. (2006) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 448, L19
  * Aharonian (2000) Aharonian F. A., 2000, NewA, 5, 377. doi:10.1016/S1384-1076(00)00039-7
  * Albert et al. (2007) Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 654, L119
  * Aleksić et al. (2015) Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2015, A&A, 573, A50. doi:10.1051/0004-6361/201322906
  * Aliu et al. (2015) Aliu, E., Archer, A., Aune, T., et al. 2015, ApJ, 799, 7
  * Antonucci, & Ulvestad (1985) Antonucci, R. R. J., & Ulvestad, J. S. 1985, ApJ, 294, 158
  * Beckmann et al. (2002) Beckmann, V., Wolter, A., Celotti, A., et al. 2002, A&A, 383, 410
  * Bhatta (2021) Bhatta, G. 2021, ApJ, 923, 7. doi:10.3847/1538-4357/ac2819
  * Bhatta & Dhital (2020) Bhatta, G. & Dhital, N. 2020, ApJ, 891, 120. doi:10.3847/1538-4357/ab7455
  * Bhatta (2019) Bhatta G., 2019, MNRAS, 487, 3990. doi:10.1093/mnras/stz1482
  * Bhatta (2018) Bhatta, G. 2018, Galaxies, 6, 136
  * Bhatta, Mohorian, & Bilinsky (2018) Bhatta, G., Mohorian M., and Bilinsky I. 2018 A&A, 619, A93
  * Bhatta & Webb (2018) Bhatta, G., & Webb, J. 2018, Galaxies, 6, 2
  * Bhatta (2017) Bhatta, G. 2017, ApJ, 847, 7
  * Bhatta et al. (2016c) Bhatta, G., Zola S., Stawarz, Ł., et al. 2016c, ApJ, 832, 47
  * Bhatta et al. (2016b) Bhatta, G., Stawarz, Ł., Ostrowski, M., et al. 2016b, ApJ, 831, 92
  * Blinov et al. (2018) Blinov, D., Pavlidou, V., Papadakis, I., et al. 2018, MNRAS, 474, 1296. doi:10.1093/mnras/stx2786
  * Baloković et al. (2016) Baloković, M., Paneque, D., Madejski, G., et al. 2016, ApJ, 819, 156. doi:10.3847/0004-637X/819/2/156
  * Bromberg & Levinson (2009) Bromberg, O. & Levinson, A. 2009, ApJ, 699, 1274. doi:10.1088/0004-637X/699/2/1274
  * B’́ottcher et al. (2007) B’́ottcher, M., Basu, S., Joshi, M., et al. 2007, ApJ, 670, 968
  * Burbidge et al. (1974) Burbidge, G. R., Jones, T. W., & Odell, S. L. 1974, ApJ, 193, 43
  * Camenzind & Krockenberger (1992) Camenzind M., Krockenberger M., 1992, A&A, 255, 59
  * Chatterjee et al. (2013) Chatterjee, R., Fossati, G., Urry, C. M., et al. 2013, ApJ, 763, L11. doi:10.1088/2041-8205/763/1/L11
  * Czerny (2006) Czerny, B. 2006, Astronomical Society of the Pacific Conference Series, 360, 265
  * Danforth et al. (2010) Danforth, C. W., Keeney, B. A., Stocke, J. T., Shull, J. M., & Yao, Y. 2010, ApJ, 720, 976
  * Dermer, Schlickeiser, & Mastichiadis (1992) Dermer C. D., Schlickeiser R., Mastichiadis A., 1992, A&A, 256, L27
  * Dodson et al. (2008) Dodson, R., Fomalont, E. B., Wiik, K., et al. 2008, ApJS, 175, 314
  * Marscher (2014) Marscher A. P., 2014, ApJ, 780, 87. doi:10.1088/0004-637X/780/1/87
  * Morozova et al. (2014) Morozova, D. A., Larionov, V. M., Troitsky, I. S., et al. 2014, AJ, 148, 42.
  * Morozova et al. (2016) Morozova, D., Larionov, V., Jorstad, S., et al. 2016, Galaxies, 4, 24.
  * Edelson et al. (2013) Edelson, R., Mushotzky, R., Vaughan, S., et al. 2013, ApJ, 766, 16
  * Emmanoulopoulos et al. (2013) Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013, MNRAS, 433, 907
  * Fan et al. (2006) Fan, J. H., Tao, J., Qian, B. C., et al. 2006, PASJ, 58, 797. doi:10.1093/pasj/58.5.797
  * Foster (1996) Foster G., 1996, AJ, 112, 1709
  * Fu & Lai (2012) Fu, W., & Lai, D. 2012, MNRAS, 423, 831
  * Ghisellini & Madau (1996) Ghisellini G., Madau P., 1996, MNRAS, 280, 67. doi:10.1093/mnras/280.1.67
  * Giannios (2013) Giannios, D. 2013, MNRAS, 431, 355. doi:10.1093/mnras/stt167
  * Giannios & Uzdensky (2019) Giannios, D. & Uzdensky, D. A. 2019, MNRAS, 484, 1378. doi:10.1093/mnras/stz082
  * Guo et al. (2014a) Guo, X., Sironi, L., & Narayan, R. 2014, ApJ, 797, 47. doi:10.1088/0004-637X/797/1/47
  * Guo et al. (2014b) Guo, X., Sironi, L., & Narayan, R. 2014, ApJ, 794, 153. doi:10.1088/0004-637X/794/2/153
  * Heidt & Wagner (1996) Heidt J., Wagner S. J., 1996, A&A, 305, 42
  * Hagen-Thorn et al. (2015) Hagen-Thorn, V. A., Larionov, V. M., Arkharov, A. A., et al. 2015, Astronomy Reports, 59, 551.
  * Hallum et al. (2022) Hallum, M. K., Jorstad, S. G., Larionov, V. M., et al. 2022, ApJ, 926, 180. doi:10.3847/1538-4357/ac4710
  * Hayashida et al. (2015) Hayashida, M., Nalewajko, K., Madejski, G. M., et al. 2015, ApJ, 807, 79
  * Hewett & Wild (2010) Hewett, P. C. & Wild, V. 2010, MNRAS, 405, 2302. doi:10.1111/j.1365-2966.2010.16648.x
  * Hughes et al. (1998) Hughes P. A., Aller H. D., Aller M. F., 1998, ApJ, 503, 662
  * IceCube Collaboration et al. (2018) IceCube Collaboration, Aartsen M. G., Ackermann M., Adams J., Aguilar J. A., Ahlers M., Ahrens M., et al., 2018, Sci, 361, eaat1378. doi:10.1126/science.aat1378
  * Isobe et al. (2015) Isobe, N., Sato, R., Ueda, Y., et al. 2015, ApJ, 798, 27
  * Jorstad et al. (2013) Jorstad, S. G., Marscher, A. P., Smith, P. S., et al. 2013, ApJ, 773, 147. doi:10.1088/0004-637X/773/2/147
  * Jorstad & Marscher (2016) Jorstad, S. & Marscher, A. 2016, Galaxies, 4, 47. doi:10.3390/galaxies4040047
  * Jorstad et al. (2017) Jorstad, S. G., Marscher, A. P., Morozova, D. A., et al. 2017, ApJ, 846, 98. doi:10.3847/1538-4357/aa8407
  * Jorstad et al. (2004) Jorstad, S. G., Marscher, A. P., Lister, M. L., et al. 2004, AJ, 127, 3115
  * Jorstad et al. (2005) Jorstad, S. G., Marscher, A. P., Lister, M. L., et al. 2005, AJ, 130, 1418
  * Jorstad et al. (2001) Jorstad, S. G., Marscher, A. P., Mattox, J. R., et al. 2001, ApJS, 134, 181
  * Kastendieck et al. (2011) Kastendieck, M. A., Ashley, M. C. B., & Horns, D. 2011, aap, 531, A123
  * Kellermann et al. (2004) Kellermann, K. I., Lister, M. L., Homan, D. C., et al. 2004, ApJ, 609, 539. doi:10.1086/421289
  * Kelly et al. (2011) Kelly, B. C., Sobolewska, M., & Siemiginowska, A. 2011, ApJ, 730, 52. doi:10.1088/0004-637X/730/1/52
  * Kirk et al. (1998) Kirk, J. G., Rieger, F. M., & Mastichiadis, A. 1998, A&A, 333, 452
  * Kluźniak & Abramowicz (2005) Kluźniak, W. & Abramowicz, M. A. 2005, Ap&SS, 300, 143. doi:10.1007/s10509-005-0020-6
  * Larsson (1996) Larsson S. 1996, A&As, 117, 197
  * Lachowicz et al. (2009) Lachowicz, P., Gupta, A. C., Gaur, H., & Wiita, P. J. 2009, A&A, 506, L17
  * Larionov et al. (2008) Larionov, V. M., Jorstad, S. G., Marscher, A. P., et al. 2008, A&A, 492, 389
  * Li & Narayan (2004) Li, L.-X., & Narayan, R. 2004, ApJ, 601, 414
  * Blandford & Eichler (1987) Blandford, R. & Eichler, D. 1987, Phys. Rep., 154, 1. doi:10.1016/0370-1573(87)90134-7
  * Lind & Blandford (1985) Lind, K. R. & Blandford, R. D. 1985, ApJ, 295, 358. doi:10.1086/163380
  * Lister et al. (1998) Lister, M. L., Marscher, A. P., & Gear, W. K. 1998, ApJ, 504, 702
  * Lomb (1976) Lomb, N. R. 1976, Ap& SS, 39, 447
  * Mannheim (1993) Mannheim K., 1993, A&A, 269, 67
  * Maraschi, Ghisellini, & Celotti (1992) Maraschi L., Ghisellini G., Celotti A., 1992, ApJL, 397, L5. doi:10.1086/186531
  * Marscher (2014) Marscher, A. P. 2014, ApJ, 780, 87
  * Marscher et al. (2010) Marscher, A. P., Jorstad, S. G., Larionov, V. M., et al. 2010, ApJ, 710, L126. doi:10.1088/2041-8205/710/2/L126
  * Marscher & Gear (1985) Marscher A. P., Gear W. K., 1985, ApJ, 298, 114. doi:10.1086/163592
  * MAGIC Collaboration et al. (2008) MAGIC Collaboration, Albert, J., Aliu, E., et al. 2008, Science, 320, 1752
  * McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083
  * Meng et al. (2018) Meng, N., Zhang, X., Wu, J., Ma, J., & Zhou, X. 2018, ApJS, 237, 30
  * Morozova et al. (2014) Morozova, D. A., Larionov, V. M., Troitsky, I. S., et al. 2014, AJ, 148, 42. doi:10.1088/0004-6256/148/3/42
  * M’́ucke et al. (2003) M’́ucke A., Protheroe R. J., Engel R., Rachen J. P., Stanev T., 2003, APh, 18, 593. doi:10.1016/S0927-6505(02)00185-8
  * Mukherjee & VERITAS Collaboration (2017) Mukherjee, R. & VERITAS Collaboration 2017, The Astronomer’s Telegram, 11075
  * Nalewajko et al. (2011) Nalewajko, K., Giannios, D., Begelman, M. C., et al. 2011, MNRAS, 413, 333. doi:10.1111/j.1365-2966.2010.18140.x
  * Narayan & Piran (2012) Narayan, R. & Piran, T. 2012, MNRAS, 420, 604. doi:10.1111/j.1365-2966.2011.20069.x
  * Nilsson et al. (2018) Nilsson, K., Lindfors, E., Takalo, L. O., et al. 2018, A&A, 620, A185
  * Nilsson et al. (2009) Nilsson, K., Pursimo, T., Villforth, C., et al. 2009, A&A, 505, 601
  * Osterman et al. (2006) Osterman, M. A., Miller, H. R., Campbell, A. M., et al. 2006, AJ, 132, 873
  * Paliya et al. (2016) Paliya, V. S., Diltz, C., B’́ottcher, M., et al. 2016, ApJ, 817, 61
  * Pandey et al. (2019) Pandey, A., Gupta, A. C., Wiita, P. J., & Tiwari, S. N. 2019, ApJ, 871, 192
  * Park et al. (2019) Park, J., Lee, S.-S., Kim, J.-Y., et al. 2019, ApJ, 877, 106. doi:10.3847/1538-4357/ab1b27
  * Papadakis et al. (2004) Papadakis, I. E., Samaritakis, V., Boumis, P., et al. 2004, A&A, 426, 437
  * Raiteri et al. (1999) Raiteri, C. M., Villata, M., Tosti, G., et al. 1999, A&A, 352, 19.
  * Raiteri et al. (2015) Raiteri, C. M., Stamerra, A., Villata, M., et al. 2015, MNRAS, 454, 353
  * Raiteri et al. (2017) Raiteri, C. M., Villata, M., Acosta-Pulido, J. A., et al. 2017, Nature, 552, 374
  * Rajput & Pandey (2021) Rajput, B. & Pandey, A. 2021, Galaxies, 9, 118. doi:10.3390/galaxies9040118
  * Rieger & Mannheim (2000) Rieger F. M., Mannheim K., 2000, A&A, 359, 948
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
  * Shah et al. (2018) Shah, Z.; Mankuzhiyil, N.; Sinha, A. et al. 2003, PASJ, 55L3
  * Sikora, Begelman, & Rees (1994) Sikora M., Begelman M. C., Rees M. J., 1994, ApJ, 421, 153. doi:10.1086/173633
  * Sikora & Begelman (2013) Sikora, M., & Begelman, M. C. 2013, ApJL, 764, L24
  * Sitko et al. (1985) Sitko, M. L., Schmidt, G. D., & Stein, W. A. 1985, ApJS, 59, 323
  * Smith et al. (1987) Smith, P. S., Balonek, T. J., Elston, R., et al. 1987, ApJS, 64, 459
  * Sheng et al. (2022) Sheng, Y., Rajagopal, M., Kaur, A., et al. 2022, arXiv:2209.09877
  * Spruit et al. (2001) Spruit, H. C., Daigne, F., & Drenkhahn, G. 2001, A&A, 369, 694. doi:10.1051/0004-6361:20010131
  * Pursimo et al. (2017) Pursimo, T., Sagues, A., Telting, J., et al. 2017, The Astronomer’s Telegram, 10948
  * Timmer & Koenig (1995) Timmer, J., & Koenig, M. 1995, A& A, 300, 707
  * Urry & Padovani (1995) Urry, C. M. & Padovani, P. 1995, PASP, 107, 803. doi:10.1086/133630
  * Vaughan (2005) Vaughan, S. 2005, A&A, 431, 391. doi:10.1051/0004-6361:20041453
  * Vaughan et al. (2003) Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P., 2003, MNRAS, 345, 1271
  * Vlasyuk et al. (2022) Vlasyuk, V. V., Spiridonova, O. I., Moskvitin, A. S., et al. 2022, The Astronomer’s Telegram, 15376
  * Wang et al. (2014) Wang, J.-Y., An, T., Baan, W. A., & Lu, X.-L. 2014, MNRAS, 443, 58
  * Wehrle et al. (2012) Wehrle, A. E., Marscher, A. P., Jorstad, S. G., et al. 2012, ApJ, 758, 72. doi:10.1088/0004-637X/758/2/72
  * Wu et al. (2006) Wu, J., Zhou, X., Wu, X.-B., et al. 2006, AJ, 132, 1256
  * Xie et al. (2008) Xie, G. Z., Yi, T. F., Li, H. Z., Zhou, S. B., & Chen, L. E. 2008, AJ, 135, 2212
  * Zhang et al. (2014) Zhang, B.-K., Zhao, X.-Y., Wang, C.-X., & Dai, B.-Z. 2014, Research in Astronomy and Astrophysics, 14, 933-941
  * Zhang et al. (2017) Zhang, P.-f., Yan, D.-h., Liao, N.-h., & Wang, J.-c. 2017, ApJ, 835, 260
  * Zheng et al. (2016) Zheng Z.-Y., Butler N. R., Shen Y. et al., 2016, ApJ, 827, 56
  * Zola et al. (2021) Zola, S., Kouprianov, V., Reichart, D. E., et al. 2021, Revista Mexicana de Astronomia y Astrofisica Conference Series, 53, 206. doi:10.22201/ia.14052059p.2021.53.40
  * Zola et al. (2016) Zola, S., Valtonen, M., Bhatta, G., et al. 2016, Galaxies, 4, 41