# Investigation of the correlation between optical and $\gamma$-ray flux
variation in the blazar Ton 599

Bhoomika Rajput1,2, Amit Kumar Mandal3, Ashwani Pandey4, C. S. Stalin5, Walter
Max-Moerbeck1, and Blesson Mathew2
1Departamento de Astronomiá, Universidad de Chile, Camino El Observatorio
1515, Las Condes, Santiago, Chile
2Department of Physics, CHRIST (Deemed to be University), Hosur Road,
Bangalore 560 029, India
3Astronomy Program, Department of Physics and Astronomy, Seoul National
University, Seoul 151-742, Republic of Korea
4Center for Theoretical Physics, Polish Academy of Sciences, Al. Lotników
32/46, 02-668 Warsaw, Poland
5Indian Institute of Astrophysics, Block II, Koramangala, Bangalore 560034,
India
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The correlation between optical and $\gamma$-ray flux variations in blazars
reveals a complex behaviour. In this study, we present our analysis of the
connection between changes in optical and $\gamma$-ray emissions in the blazar
Ton 599 over a span of approximately 15 years, from August 2008 to March 2023.
Ton 599 reached its highest flux state across the entire electromagnetic
spectrum during the second week of January 2023. To investigate the connection
between changes in optical and $\gamma$-ray flux, we have designated five
specific time periods, labeled as epochs A, B, C, D, and E. During periods B,
C, D, and E, the source exhibited optical flares, while it was in its
quiescent state during period A. The $\gamma$-ray counterparts to these
optical flares are present during periods B, C, and E, however during period
D, the $\gamma$-ray counterpart is either weak or absent. We conducted a
broadband spectral energy distribution (SED) fitting by employing a one-zone
leptonic emission model for these epochs. The SED analysis unveiled that the
optical-UV emission primarily emanated from the accretion disk in quiescent
period A, whereas synchrotron radiation from the jet dominated during periods
B, C, D, and E. Diverse correlated patterns in the variations of optical and
$\gamma$-ray emissions, like correlated optical and $\gamma$-ray flares, could
be accounted for by changes in factors such as the magnetic field, bulk
Lorentz factor, and electron density. On the other hand, an orphan optical
flare could result from increased magnetic field and bulk Lorentz factor.

###### keywords:

galaxies: active - galaxies: nuclei - galaxies:jets - individual: Ton 599 -
$\gamma$-ray:galaxies

††pubyear: 2015††pagerange: Investigation of the correlation between optical
and $\gamma$-ray flux variation in the blazar Ton 599–Investigation of the
correlation between optical and $\gamma$-ray flux variation in the blazar Ton
599

##  1 Introduction

Blazars, the peculiar category of active galactic nuclei (AGN) known for
emitting non-thermal variable radiation covering the whole range of the
electromagnetic spectrum, are among the most luminous objects in the Universe
(luminosity $\approx$ 1042$-$1048 erg s-1). These objects emit relativistic
jets that are oriented in the same direction as the observer’s line of sight
(Antonucci, 1993; Urry & Padovani, 1995), resulting in strong Doppler
boosting. They reside within elliptical galaxies (Urry et al., 2000) and their
energy source is derived from the process of material accreting onto a
supermassive black hole. (106-1010 M⊙)(Shakura & Sunyaev, 1973; Shaw et al.,
2012).

Blazars are categorized into two groups: flat spectrum radio quasars (FSRQs)
and BL Lacs. This classification relies on the rest-frame equivalent width
(EW) of the emission lines observed in their optical spectra. FSRQs exhibit
strong emission lines (EW $>$ 5 Å), while BL Lacs exhibit either weak or no
emission lines (EW $<$ 5 Å) (Stocke et al., 1991; Scarpa & Falomo, 1997). The
absence of broad emission lines in the optical spectra of BL Lac objects may
be caused by Doppler-boosted non-thermal continuum, which could swamp out the
spectral emission lines (Blandford & Rees, 1978). Additionally, the presence
of EW $>$ 5 Å observed in BL Lac objects might be the outcome of an especially
low state of the beamed continuum (Vermeulen et al., 1995). Consequently, the
EW by itself is not a reliable indicator of the difference between the two
categories of blazars. Ghisellini et al. (2011) and Sbarrato et al. (2012)
proposed a categorization scheme for FSRQs and BL Lacs that hinges on the
luminosity of the broad-line region (BLR), which is measured in Eddington
units. In accordance with this scheme, the luminosity of the broad line region
(LBLR) to the Eddington luminosity ratio for FSRQs is greater than 5 $\times$
10-4 while it is less than 5 $\times$ 10-4 for BL Lacs.

The broadband spectral energy distribution (SED) of blazars has a two-hump
structure, where one hump peaks at lower energies in the optical/IR/X-ray
bands and the second hump peaks at higher energies in the MeV/GeV bands
(Fossati et al., 1998; Mao et al., 2016). The low-energy emission is
attributed to the synchrotron emission process of the relativistic electrons
in the jets (Urry & Mushotzky, 1982; Impey & Neugebauer, 1988). Padovani &
Giommi (1995) and Abdo et al. (2010c) proposed an alternative method of
categorizing blazars by considering the position of the synchrotron peak
frequency ($\nu^{\rm syn}_{\rm peak}$) in their SED. These classes include
low-synchrotron peaked blazars (LSP, $\nu^{\rm syn}_{\rm peak}$ $<$ 1014 Hz),
intermediate-synchrotron peaked blazars (ISP, 1014 Hz $<$ $\nu^{\rm syn}_{\rm
peak}$ $<$ 1015 Hz), and high-synchrotron peaked blazars (HSP, 1015 Hz $<$
$\nu^{\rm syn}_{\rm peak}$ $<$ 1017 Hz). Most FSRQs are LSP blazars, while BL
Lac sources exhibit all three blazar behaviors. Some BL Lacs are extremely
high-synchrotron peaked blazars (EHBL or extreme blazars) with $\nu^{\rm
syn}_{\rm peak}$ $>$ 1017 Hz (Foffano et al., 2019).

It is renowned that the synchrotron emission process results in low-energy
emission. However, the origin of the higher energy peak remains partially
understood. To explain this numerous models have been developed. Leptonic
model often provides a satisfactory explanation for the observed broadband SED
of the majority of blazars. In the leptonic scenario, the high-energy emission
is produced by the inverse compton (IC) process (Abdo et al., 2010c). The seed
photons that cause the IC process because of the scattering of the
relativistic electron in the jets can originate from a variety of locations.
These could originate from within the jets by synchrotron emission process
(synchrotron self-Compton or SSC; Konigl 1981; Marscher & Gear 1985;
Ghisellini & Maraschi 1989; Tramacere et al. 2009) or from outside the jets
(external Compton or EC) mainly from the accretion disk (Boettcher et al.,
1997; Dermer & Schlickeiser, 2002), from the BLR (Ghisellini & Madau, 1996;
Dermer et al., 2009) and from the dusty torus (Błażejowski et al., 2000;
Sikora et al., 2008).

In addition to the leptonic model, the high-energy peak in the broadband SED
of blazars can also be accounted for by the hadronic model (Böttcher et al.,
2013). Hadronic model for blazar postulates that protons are accelerated to
ultra-relativistic energies and dominate the high-energy emission through
proton-synchrotron radiation (Aharonian, 2000; Mücke et al., 2003) or photo-
pion production via proton-proton or proton-photon interaction (Mannheim,
1993), which results in the production of ultra-high-energy photons and
neutrinos. The hybrid lepto-hadronic model has also been proposed in some
cases to explain high-energy emission in blazar (Diltz & Böttcher, 2016;
Paliya et al., 2016; Cerruti et al., 2019).

In addition to the broadband SED, the study of the temporal correlation
between the two major components of the broadband emission i.e. optical and
$\gamma$-ray emission can give major insights to the high-energy emission in
blazars. In the leptonic model of blazar emission, since the relativistic
electrons in the jets are responsible for both the optical and $\gamma$-ray
emission (Böttcher, 2007), a close correlation between the optical and
$\gamma$-ray flux changes is anticipated. The launch of the Fermi Gamma-ray
space telescope111https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ (referred
to as Fermi hereafter; Atwood et al. 2009) in 2008 gave us unusual opportunity
to characterise the $\gamma$-ray for a large number of sources (Abdollahi et
al., 2020). To investigate the connection between the optical and $\gamma$-ray
bands, Fermi data was combined with numerous ground-based optical and IR
observations. Several studies have been conducted on this topic, yielding
diverse conclusions. For example, correlations between optical and
$\gamma$-ray changes, with or without lag, have been found in blazars (Bonning
et al., 2009; Chatterjee et al., 2012; Liao et al., 2014; Liodakis et al.,
2018). These findings firmly support the leptonic model of blazar emission.
However, many studies disfavour the one-zone leptonic model due to the absence
of a correlation between changes in optical and $\gamma$-ray flux. In these
studies, either $\gamma$-ray flares are present without their corresponding
optical flares (orphan $\gamma$-ray flares) (Dutka et al., 2013; Cohen et al.,
2014; MacDonald et al., 2015; Rajput et al., 2020) or optical flares are
present without their $\gamma$-ray counterparts (orphan optical flares)
(Chatterjee et al. 2013; Rajput et al. 2019, 2020 and references therein).

The existing observations suggest that the correlation between variations in
optical flux and $\gamma$-ray flux is intricate. In this work, we conducted a
thorough examination over the long-term duration to understand the diverse
behaviour between changes in optical and $\gamma$-ray flux through broadband
SED modeling in the blazar Ton 599. This blazar was recently observed at its
brightest state (Prince, 2023; Garrappa & Valverd, 2023; Tripathi et al.,
2023). The availability of multi-wavelength data over the long-term period
renders this source appropriate candidate for conducting study on the optical-
GeV connection. Ton 599, also known as 4C $+$29.45, is named as 4FGL
J1159.5$+$2914 in the 4th catalog of Fermi Large Area Telescope (4FGL;
Abdollahi et al. 2022). It is located at a redshift z = 0.725 (Hewett & Wild,
2010), with RA = 179∘.882641 and Dec = 29∘.245507. This source is Classified
as an FSRQ source and it falls under the category of LSP source (Abdo et al.,
2010a). It displays high optical variability and exhibits strong polarization
(Fan et al., 2006). In the $\gamma$-ray band this source was first detected by
the Energetic Gamma Ray Experiment Telescope (EGRET) (Thompson et al., 1995)
and later by the Fermi (Abdo et al., 2010b). Moreover, this source was
detected in the very high-energy band (VHE $>$ 100 GeV) by VERITAS (Mukherjee
& VERITAS Collaboration, 2017). This source is observed with significant
variability in the optical and $\gamma$-ray, as well as in other bands of the
electromagnetic spectrum (Fan et al., 2006; Patel et al., 2018; Prince, 2019;
Rajput & Pandey, 2021; Bhatta et al., 2023). In 2017, Ton 599 was detected in
an unprecedented massive flaring state across the entire electromagnetic
spectrum (Cheung et al., 2017; Carrasco et al., 2017; Acosta Pulido et al.,
2017). The optical and GeV connection for this source during its 2017 outburst
was investigated by Prince (2019), and optical and $\gamma$-ray emission
regions are found to be co-spatial with the time lag of a few days. This
source was modelled using a two-zone leptonic emission model during the same
flare by Patel & Chitnis (2020), and it was found that the GeV emission was
from EC process, where the seed photons originated from the dusty torus. The
largest flare ever recorded by Fermi for this source in a period of around 15
years occurred in January 2023 when Ton 599 was once again found to be
flaring, with a flux value of 3.60 $\pm$ 0.27 $\times$ 10-6 ph cm-2 s-1. Here,
we present the findings of a multi-wavelength analysis performed on the blazar
Ton 599 using data collected over a period of about 15 years, from August 2008
to March 2023. The main objective of this analysis is to determine any
potential correlation between variations in the optical V-band flux and
$\gamma$-ray flux and, as a result, putting constraints on the emission
process in blazars. We provide the information about the data utilized in this
work in section 2. Section 3 provides a detail of the analysis, and section 4
presents the findings and a discussion. The summary of the work is given in
the final section.

Figure 1: Multi-wavelength light curves of the source Ton 599 for the period
from MJD 54686 to MJD 60008. The top panel refers to the one day binned
$\gamma$-ray (0.1-300 GeV) light curve, the second panel refers to the X-ray
(0.3-10 keV) light curve, the third panel refers to the UV light curve, the
fourth and fifth panels refer to the optical light curves. The $\gamma$-ray
fluxes are measured in the units of 10-6 ph cm-2 s-1 and the UV and optical
fluxes are measured in the units of 10-11 erg cm-2 s-1. The two vertical black
dashed lines represent the quiescent epoch A. The vertical blue dashed lines
refer to the flaring epochs B, C and E, when both optical and $\gamma$-ray
flares are present. The two vertical green dashed lines represent the flaring
epoch D, when optical flare is present without its corresponding $\gamma$-ray
counterpart. The peaks of optical and $\gamma$-ray flares are denoted by the
vertical cyan dashed lines. The horizontal black dashed lines in the
$\gamma$-ray and optical V-band light curves represents the average flux from
August 2008 to March 2023. Table 1: Details of the epochs studied in this work
for optical-GeV correlation and broadband SED modelling. In this table the
optical flux values are in units of 10-11 erg cm-2 s-1, and the $\gamma$-ray
flux values are in units of 10-6 ph cm-2 s-1.

|  | Time-period (MJD) | Calender date | Average Flux | Peak Flux
---|---|---|---|---|---
Epoch | Status | Start | End | Start | End | Optical | $\gamma$ | Optical | $\gamma$
A | Quiescent | 55650 | 55750 | 2011-03-30 | 2011-07-08 | 0.37$\pm$0.02 | 0.11$\pm$0.02 | 0.56$\pm$0.03 | 0.21$\pm$0.09
B | Optical flare with $\gamma$-ray counterpart | 58007 | 58107 | 2017-09-11 | 2017-12-20 | 3.89$\pm$0.02 | 0.99$\pm$0.02 | 5.47$\pm$0.13 | 2.49$\pm$0.28
C | Optical flare with $\gamma$-ray counterpart | 59355 | 59455 | 2021-05-21 | 2021-08-29 | 3.62$\pm$0.03 | 0.63$\pm$0.02 | 5.60$\pm$0.14 | 2.42$\pm$0.28
D | Optical flare without $\gamma$-ray counterpart | 59500 | 59600 | 2021-10-13 | 2022-01-21 | 3.03$\pm$0.02 | 0.45$\pm$0.01 | 4.67$\pm$0.12 | 0.81$\pm$0.05
E | Optical flare with $\gamma$-ray counterpart | 59905 | 60005 | 2022-11-22 | 2023-03-02 | 5.55$\pm$0.03 | 1.32$\pm$0.03 | 10.30$\pm$0.14 | 3.60$\pm$0.27

$\begin{array}[]{rrr}\includegraphics[width=193.47873pt,height=256.0748pt]{lightcurve_epoch_a.pdf}\includegraphics[width=193.47873pt,height=256.0748pt]{lightcurve_epoch_b.pdf}\includegraphics[width=193.47873pt,height=256.0748pt]{lightcurve_epoch_c.pdf}\end{array}$

$\begin{array}[]{rr}\includegraphics[width=193.47873pt,height=256.0748pt]{lightcurve_epoch_d.pdf}\includegraphics[width=193.47873pt,height=256.0748pt]{lightcurve_epoch_e.pdf}\end{array}$

Figure 2: Multi-wavelength light curves for the selected epochs A, B, C, D,
and E. The light curves for epochs A, B, and C are displayed in the upper
left, middle, and right panels, respectively. The bottom left and right panels
represent the light curves for epochs D and E, respectively. The detail of the
panels are same as in Figure 1.

##  2 Multi-wavelength data reduction

To carry out a multi-wavelength analysis of the source Ton 599, we used
publicly accessible data in the $\gamma$-ray, X-ray, UV, and optical bands
that span around 15 years, from August 2008 to March 2023.

###  2.1 $\gamma$-ray data

We analyzed about 15 years of data collected from the Large Area Telescope
(LAT) aboard the Fermi. The Fermi-LAT is a pair-conversion telescope and has a
very large effective area ($>$ 8000 cm2) (Atwood et al., 2009). It usually
operates in the scanning mode and scans the entire sky once every $\sim$ 3
hours covering an extensive spectral energy range from 20 MeV to TeV energies.
For this study, we performed binned likelihood analysis by using the Fermi-LAT
data in the energy range of 100 MeV-300 GeV from 2008 August 08 to 2023 March
05, a period of approximately 15 years (MJD 54686$-$60008). To analyze the
Fermi data we used the current version of Fermitools 2.2.0. We used PASS 8
data where the photon-like events are categorized as ’evclass=128, evtype=3’.
All the events were chosen within a 10∘ circular region of interest (ROI)
centered on the source and a 90∘ zenith angle cut was made to eliminate
background $\gamma$-rays from the Earth’s limb. After that, we generated a
good time interval using the suggested criteria "(DATA QUAL $>$ 0)&&(LAT
CONFIG==1)". For this analysis, we used the latest instrument response
function (IRF) "P8R3_SOURCE_V3", Galactic diffuse emission model "gll_iem_v07"
and the extragalactic isotropic diffuse emission model
"iso_P8R3_SOURCE_V3_v1". Our source model includes all the point sources from
the 4FGL catalog that fall within 10∘ of the source, as well as the Galactic
and the extragalactic isotropic diffuse emission components. The spectral
shapes of all the sources are adopted from the 4FGL catalog. For the sources
lying between 0∘ to 5∘, the spectral parameters are left free except the
scaling factor and for the source lying between 5∘ to 10∘, all the spectral
parameters are kept fixed to their 4FGL catalog value. To determine the
significance of the detection, we used the maximum likelihood (ML) test
statistic, defined as TS = 2$\Delta$ log(L), where L represents the likelihood
function comparing models with and without a $\gamma$-ray point source at the
source’s location. During the likelihood analysis, the sources with TS $<$ 25,
which corresponds to a 5$\sigma$ detection (Mattox et al., 1996), were
removed. This modified model is then employed for the analysis of the light
curve and spectra. To generate the one-day binned $\gamma$-ray light curve
over the time period of interest, we regarded the source as detectable when
the value of TS $>$ 9, equivalent to a 3$\sigma$ detection.

###  2.2 X-ray data

####  2.2.1 Swift-XRT data

The X-ray telescope (XRT), as detailed in Burrows et al. (2005), mounted on
the Swift satellite (Gehrels et al., 2004), operates within the 0.3 to 10 KeV
energy range, observed Ton 599 concurrently with Fermi. The X-ray data
utilized in this work was taken from the archives at
HEASARC222https://heasarc.gsfc.nasa.gov/docs/archive.html. To analyze the
data, we adhered the standard procedure given by the instrument pipeline. We
used the data obtained only in the photon counting (PC) mode for the light
curve and spectrum analysis. We analyzed the data with the xrtpipeline task,
utilizing the latest CALDB and response files provided in HEASOFT version
6.29. We employed the standard grade selection 0-12. To generate the energy
spectra the calibrated and cleaned events files were added. The source spectra
were extracted from a circular area with a radius of 60′′, positioned at the
center of the source, whereas, the background spectra were extracted from a
region with a radius of 80′′, located away from the source. The exposure maps
were combined using the ximage task, and the ancillary response files (ARFs)
were created using the xrtmkarf task. The corresponding response matrix files
(RMFs) were obtained from CALDB. The RMFs and ARFs files were loaded into the
grppha tool together with the source and background spectra in order to
combine and group them. To perform the fitting of the final X-ray spectra
within XSPEC (Arnaud, 1996), we employed an absorbed simple power law model
with a Galactic neutral hydrogen column density NH = 1.63$\times$1020 cm-2
(Dickey & Lockman, 1990; Kalberla et al., 2005).

####  2.2.2 NuSTAR data

NuSTAR, the first focusing hard X-ray telescope, was launched in June 2012. It
consists of two co-aligned X-ray detectors pair with their corresponding focal
plane modules, called FPMA and FPMB and it operates over a wide energy range
from 3 to 79 keV (Harrison et al., 2013). NuSTAR recorded two observations of
Ton 599 in May 2019 and June 2021. For this study, we used one observation
(ObsID 60463037004) during the 2021 flare of Ton 599 with an exposure time of
17.5 ks. The data was reduced using the NuSTAR Data Analysis Software package
NuSTARDAS within the
HEASOFT333http://heasarc.gsfc.nasa.gov/docs/nustar/analysis/ version 6.29. The
calibrated and cleaned level 2 event files were produced using the nupipeline
task and the CALDB version 20211202. We extracted the source and background
spectra using a circular region of 30′′ with the nuproducts script. We grouped
the spectra using grppha with at least 20 counts per bin and further processed
it into XSPEC.

###  2.3 UV and optical data

We obtained the data in the UV and optical bands from Swift-UV-Optical
telescope (Swift-UVOT). We analyzed this data using the online
tool444https://www.ssdc.asi.it/cgi-bin/swiftuvarchint. To generate the light
curves in the UV and optical bands, the magnitudes obtained using the online
tool were converted into fluxes using the zero point method described in
Breeveld et al. (2011) without correcting for galactic reddening. However, to
examine the spectral characteristics of the source, the magnitudes were
corrected for galactic extinction (Aλ), whose values were obtained from
NED555https://ned.ipac.caltech.edu/. In addition to using optical data
obtained from Swift-UVOT, we also incorporated optical V-band data from the
ground-based survey ASAS-SN666https://asas-sn.osu.edu/ and from the Steward
observatory (SO)777http://james.as.arizona.edu/$\sim$psmith/Fermi/. The ASAS-
SN ground-based survey operates with five stations positioned in both the
northern and southern hemispheres. It has the capability to conduct daily
observations of the entire visible sky, reaching a depth of g = 18.5 mag
(Shappee et al., 2014; Kochanek et al., 2017). The Steward observatory, which
is located in Arizona, employs the 2.3 m Bok Telescope and the 1.54 m Kuiper
Telescope for photometric observations.

##  3 Analysis

###  3.1 Multi-wavelength light curves

We have generated multi-wavelength light curves of the source Ton 599 for a
period of about 15 years from August 08, 2008 to March 05, 2023 (MJD
54686$-$60008). The multi-wavelength light curves, which include 1-day binned
$\gamma$-ray, X-ray (PC mode), UV and optical light curves, are shown in
Figure 1. From Figure 1 it is conspicuous that Ton 599 has numerous quiescent
and flaring stages during this period. We have identified specific time
periods within this 15-year duration. These include one quiescent period
labeled as epoch A and four flaring periods (when an optical or $\gamma$-ray
flare is present), marked as epochs B, C, D, and E. These quiescent and
flaring epochs have a 100-day time span and flaring epochs are centered on the
peak of either the $\gamma$-ray or optical light curves. The time periods of
these epochs are given in Table 1. These epochs were chosen based on the
following criteria: 1) The accessibility of multi-wavelength data 2) optical
and $\gamma$-ray light-curve data during the quiescent phase should be below
the average level 3) optical or $\gamma$-ray flare should gradually increase
by more than two times from the average level. The values of average optical
and $\gamma$-ray fluxes are 1.24$\pm$0.003 $\times$ $10^{-11}$ erg cm-2 s-1
and 0.42$\pm$0.003 $\times$ $10^{-6}$ ph cm-2 s-1, respectively, over the
period from August 2008 to March 2023. The average and peak values of optical
and $\gamma$-ray flux for each epoch are given in Table 1. The details of each
of the epochs are given below.

####  3.1.1 Epoch A (MJD 55650$-$55750)

The source was in a quiescent state in the $\gamma$-ray, X-ray, UV and optical
bands during this epoch. The multi-wavelength light curve of this epoch is
shown in the upper left panel of Figure 2.

####  3.1.2 Epoch B (MJD 58007$-$58107)

During this epoch, the source showed major variation in $\gamma$-ray flux. The
$\gamma$-ray flux increased by about six times from the average level during
this epoch. Along with the large changes in the $\gamma$-ray flux, changes in
the X-ray, UV and optical fluxes were also observed (refer to upper middle
panel of Figure 2). The optical flux also increased by nearly 5 times from the
average level. During this epoch, two short-term intense $\gamma$-ray flares
were observed superimposed on the large $\gamma$-ray flare at around MJD 58057
and MJD 58077. Visual inspection of the Figure 2 upper middle panel) also
gives the hint of corresponding short-term optical flare at around MJD 58057
but there is no hint of an optical flare at around MJD 58077 due to lack of
data points. We deduce that during this time period, we noted the same
patterns of flux variations in both the optical and GeV energy ranges.

####  3.1.3 Epoch C (MJD 59355$-$59455)

The source experienced a significant $\gamma$-ray flare during this epoch,
this one being around six times greater than the $\gamma$-ray flux on average.
This substantial flaring event was subsequently accompanied by flares in the
X-ray, UV, and optical bands. Notably, the optical flux experienced an
increase of approximately five times compared to its average level. Thus, in
this epoch too, we found that the variability pattern in the optical and
$\gamma$-ray bands appeared to be correlated. The multi-wavelength light curve
for this epoch is shown in the upper right panel of Figure 2.

####  3.1.4 Epoch D (MJD 59500$-$59600)

During this epoch, the $\gamma$-ray flux variations are moderate but the
changes in the X-ray, UV and optical flux are significant (see the lower left
panel of the Figure 2). In comparison to the average flux level, the
$\gamma$-ray flux increased by less than two times, whereas the optical flux
increased by about four times compared to the average flux level. Since there
was no clear large amplitude variation in the $\gamma$-ray light curve at this
epoch, we can infer that in this epoch we observed an orphan optical flare.

####  3.1.5 Epoch E (MJD 59905$-$60005)

The source underwent the major $\gamma$-ray flare during this epoch, which is
the largest $\gamma$-ray flare over the whole time period examined in this
work. Such large flux variations were also noticed in the X-ray, UV and
optical bands. The multi-wavelength light curve for this epoch is shown in
lower right panel of Figure 2. During this epoch the $\gamma$-ray flux
increased about nine times from the average flux level and the optical flux
increased nearly eight times compared to the average flux. We conclude that in
this epoch the optical and $\gamma$-ray flux changes are structurally
correlated.

$\begin{array}[]{rr}\includegraphics[width=213.39566pt,height=213.39566pt]{dcf_epoch_b.pdf}\includegraphics[width=213.39566pt,height=213.39566pt]{dcf_epoch_c.pdf}\end{array}$

$\begin{array}[]{rr}\includegraphics[width=213.39566pt,height=213.39566pt]{dcf_epoch_d.pdf}\includegraphics[width=213.39566pt,height=213.39566pt]{dcf_epoch_e.pdf}\end{array}$

Figure 3: The outcomes of the discrete cross-correlation analysis conducted on
the optical V-band and $\gamma$-ray light curves during epochs B, C, D, and E.
The DCF is indicated by the black dashed line. The significance levels are
represented by the yellow dashed-dotted line (68.27%, 1$\sigma$), green dotted
line (95.45%, 2$\sigma$), and red dashed line (99.73%, 3$\sigma$). A negative
lag indicates the $\gamma$-ray flux leads the optical flux variations.

###  3.2 Cross-correlation analysis

The initial investigation described in Section 3.1 involved visually examining
the correlation between changes in optical and $\gamma$-ray emissions over
different time periods. Furthermore, we explored this correlation in depth by
employing advanced cross-correlation analysis methods. Analyzing the cross-
correlation between variations in flux across various energy bands allows for
the determination of whether emission in different bands are co-spatial or if
they emanate from different regions in the jet. In this study, we conduct a
cross-correlation analysis using the discrete cross-correlation function (DCF)
as outlined in Edelson & Krolik (1988)) to measure the time lag between
unevenly sampled optical (V-band) light curve and $\gamma$-ray light curve
binned over one day. The significance of the cross-correlation is evaluated
through Monte Carlo simulations, with the assumption of a simple power-law
power spectral density (PSD) model. This model is defined as P($\nu$)
$\propto$ 1/$\nu^{\beta}$, where $\nu$ represents the temporal frequency, and
$\beta$ denotes the slope. A comprehensive explanation of this methodology is
provided by Max-Moerbeck et al. (2014b). We have adopted the value of
$\beta_{\gamma}$ = 1.5 (Max-Moerbeck et al., 2014a) and $\beta_{\rm{optical}}$
= 2 (Goyal, 2021). The DCF between optical and $\gamma$-ray light curves for
epochs B, C, D, and E are shown in Figure 3. During these epochs, we detected
no significant correlations between variations in optical and $\gamma$-ray
flux. The results of the DCF analysis are provided in Table 2.

Table 2: The outcomes of the DCF analysis between optical flux and $\gamma$-ray flux variations during epochs B, C, D, and E. Here $\tau$ represents the time lag, with negative values indicating that the $\gamma$-ray is leading the optical flux variations. The DCF values provide estimates of the cross-correlation coefficient, the p-value represents the likelihood of observing the correlation under the assumption of the null hypothesis. Epoch |  $\tau$ (days) | DCF | p-value
---|---|---|---
B | -10.0$\pm$6.1 | 0.37 | 0.365
C | 0.0$\pm$3.7 | 0.56 | 0.139
D | -20.0$\pm$14.8 | 0.31 | 0.287
E | 0.0$\pm$3.7 | 0.59 | 0.103

###  3.3 Optical spectral variations

Apart from changes in flux, blazars are also known for exhibiting variations
in their spectral characteristics. In general, FSRQs exhibit redder-when-
brighter (RWB) trend (Gu et al., 2006; Bonning et al., 2012; Zhang et al.,
2015; Rajput et al., 2019; Safna et al., 2020; Negi et al., 2022). However
bluer-when-brighter (BWB) pattern is also observed in FSRQ sources (Wu et al.,
2011; Rajput et al., 2019; Safna et al., 2020; Rajput et al., 2020). The
optical and UV spectral regions (V, B, U, W1, M2, and W2 bands) investigated
in this paper include contributions from both non-thermal emission from the
relativistic jet and thermal emission from the accretion disk. Hence, by
examining spectral variations, one can identify the many factors causing the
observed flux variations. To describe the spectral behaviour of the epochs
studied in this work we looked for the variations in the B-V colour against
the V-band magnitude. The colour-magnitude diagrams of these epochs are shown
in Figure 4. Due to lack of data, we have not included the color-magnitude
diagram for epoch A. The trend of spectral behaviour during epochs B, C, D and
E was identified by weighted linear least squares fit, which take into account
errors in both colour and V-band magnitude. We determined the Spearman rank
correlation coefficient (R) and probability of obtaining this correlation
under the assumption of the null hypothesis (p-value). If R $>$ 0.5 or R $<$
$-$0.5 with p-value $<$ 0.05, the source is considered to have displayed any
significant colour trend over the epochs. We found no spectral variations
during the epochs considered in this work. The results of the weighted linear
least squares fitting are given in Table 3.

$\begin{array}[]{rr}\includegraphics[width=213.39566pt,height=184.9429pt]{colour_magnitude_B.pdf}\includegraphics[width=213.39566pt,height=184.9429pt]{colour_magnitude_C.pdf}\end{array}$

$\begin{array}[]{rr}\includegraphics[width=213.39566pt,height=184.9429pt]{colour_magnitude_D.pdf}\includegraphics[width=213.39566pt,height=184.9429pt]{colour_magnitude_E.pdf}\end{array}$

Figure 4: Colour-magnitude diagram for Epochs B, C, D and E. The cyan line depicts the weighted linear least squares fit to the data. Table 3: Results of the weighted linear least squares fit and Spearman rank correlation test applied to the colour-magnitude diagram during epochs B, C, D, and E. None of the correlations are significant for a chosen significance level of 0.05. Epoch | Intercept | Slope | R | p-value
---|---|---|---|---
B | 2.80$\pm$1.64 | -0.15$\pm$0.12 | -0.18 | 0.52
C | 3.84$\pm$0.63 | -0.23$\pm$0.04 | -0.08 | 0.76
D | -0.47$\pm$1.57 | 0.07$\pm$0.11 | 0.15 | 0.56
E | 1.11$\pm$0.46 | -0.04$\pm$0.03 | -0.28 | 0.31

###  3.4 $\gamma$-ray spectral analysis

The inherent distribution of the emitting electrons in blazars may be
reflected in the curvature of the $\gamma$-ray spectra, which also provides
information on the potential acceleration and cooling processes in the jets.
The $\gamma$-ray spectra of the blazars is either fit by power law (PL) or log
parabola (LP) models. The PL model is defined as:

$dN(E)/dE=N_{\circ}(E/E_{\circ})^{{\Gamma_{\rm{P}}}}$ (1)

Here dN(E)/dE is the number of photons in cm-2 s-1 MeV-1, N∘ is the
normalization factor of the energy spectrum, E∘ is the scaling factor and
$\Gamma_{\rm{P}}$ is the photon index.
The LP model is defined as

$dN(E)/dE=N_{\circ}(E/E_{\circ})^{-\alpha-\beta ln(E/E_{\circ})}$ (2)

Here, $\alpha$ is the photon index at E and $\beta$ is the curvature index
that determines the curvature around the peak.
To evaluate the model that more accurately characterizes the $\gamma$-ray
spectra (PL against LP), we applied the maximum likelihood estimation
technique outlined in Mattox et al. (1996), utilizing the analysis tool
gtlike. The gtlike tool maximizes the likelihood by employing different
optimizers, which determine the optimal spectral parameter values for the best
fit. We also calculated the TScurve value, as specified in Nolan et al.
(2012), using the following equation:

$TS_{\rm{curve}}=2(logL_{\rm{LP}}-logL_{\rm{PL}})$ (3)

In the above equation, L is the likelihood function. In order to detect
significant curvature in the $\gamma$-ray spectra, we employed the threshold
value TScurve $>$ 16 (at 4$\sigma$ level; Mattox et al. 1996). The results of
the $\gamma$-ray spectral analysis are given in Table 4 and the $\gamma$-ray
spectra for the selected epochs are shown in Figure 5. The $\gamma$-ray
spectra are well fit with the LP function for all the epochs, except the
quiescent epoch A, where the value of TScurve is less than 16.

$\begin{array}[]{rrr}\includegraphics[width=170.71652pt,height=142.26378pt]{TON599_spec_a.pdf}\includegraphics[width=170.71652pt,height=142.26378pt]{TON599_spec_b.pdf}\includegraphics[width=170.71652pt,height=142.26378pt]{TON599_spec_c.pdf}\end{array}$

$\begin{array}[]{rr}\includegraphics[width=170.71652pt,height=142.26378pt]{TON599_spec_d.pdf}\includegraphics[width=170.71652pt,height=142.26378pt]{TON599_spec_e.pdf}\end{array}$

Figure 5: PL and LP fits to the $\gamma$-ray spectra for epochs A, B, C, D and
E. Among these, the PL fit provides a better description of the $\gamma$-ray
spectra in epoch A, while the LP fit better explains the $\gamma$-ray spectra
in epochs B, C, D, and E. Refer to Table 4 for detailed results. Table 4:
Spectral analysis results for the $\gamma$-ray observations. The $\gamma$-ray
flux values are expressed in units of 10-6 ph cm-2 s-1. The PL and LP fits to
the $\gamma$-ray spectra during these time periods are shown in Figure 5.

| PL | LP |
---|---|---|---
Epochs | $\Gamma_{\rm{P}}$ | Flux | TS | -Log LPL | $\alpha$ | $\beta$ | Flux | TS | -Log LLP | TScurve
A | 2.22$\pm$0.01 | 0.04$\pm$0.001 | 117.50 | 315166.67 | 2.16$\pm$0.16 | 0.42$\pm$0.78 | 0.04$\pm$0.008 | 117.82 | 315166.50 | 0.34
B | 1.95$\pm$0.01 | 0.92$\pm$0.01 | 19702.7 | 385467.61 | 1.88$\pm$0.002 | 0.49$\pm$0.01 | 0.88$\pm$0.003 | 19754.5 | 385446.44 | 42.34
C | 2.05$\pm$0.02 | 0.61$\pm$0.01 | 10498.9 | 406121.01 | 1.96$\pm$0.02 | 0.82$\pm$0.13 | 0.58$\pm$0.01 | 10552.7 | 406094.64 | 52.74
D | 2.04$\pm$0.001 | 0.49$\pm$0.001 | 8837.58 | 415533.94 | 1.97$\pm$0.02 | 0.56$\pm$0.12 | 0.47$\pm$0.01 | 8862.74 | 415521.70 | 24.48
E | 2.05$\pm$0.001 | 1.27$\pm$0.007 | 23812.1 | 378575.73 | 1.97$\pm$0.02 | 0.83$\pm$0.09 | 1.20$\pm$0.02 | 23812.4 | 378520.98 | 109.50

###  3.5 Broadband SED modelling

The source has shown varied flux variability behaviour over the last 15 years,
during which we have detected instances of correlated optical and GeV
variations as well as optical flare without $\gamma$-ray counterpart. To
understand the diverse behaviour of the source during different epochs we
constructed broadband spectral energy distributions and modelled them using
one-zone leptonic emission model. To model the broadband SED, we used the
publicly available code, JetSet (Tramacere et al., 2009, 2011; Tramacere,
2020). JetSet (version 1.2.2) fits numerical models to the data in order to
determine the optimal parameter values that best describe the observed data.

The photometric observations taken in both UV and optical wavelengths during
each epoch were averaged filter-wise to generate a single photometric data
point for each epoch. On the other hand, the average spectra for X-rays and
$\gamma$-rays were created using all of the data that was available for the
100-day span in each of the epochs. Since the EC process primarily accounts
for the high-energy $\gamma$-ray emission in the FSRQs (Paliya et al., 2015;
Shah et al., 2017; Rajput et al., 2019, 2020), we used synchrotron, SSC, and
EC processes to model our SEDs. In the one-zone leptonic emission model, the
spherical blob of radius R, located at a distance of RH form the central black
hole. It is postulated that this spherical region is injected with non-
thermally accelerated electrons that follow a distribution characterized by a
broken power law. The broken power law distribution of electrons is defined as
follows:

$N(\gamma)d\gamma=\begin{cases}K\gamma^{\rm{-p}}d\gamma,&\gamma_{\rm{min}}<\gamma<\gamma_{\rm{break}}\\\
K\gamma_{\rm{break}}^{\rm{p_{1}-p}}\gamma^{\rm{-p_{1}}}d\gamma,&\gamma_{\rm{break}}<\gamma<\gamma_{\rm{max}}\end{cases}$
(4)

Where, $\gamma$ is the electron Lorentz factor, N is the electron density in
units of cm-3 and K is the normalization constant. p and p1 are the low and
high energy spectral slopes, and $\gamma_{\rm{min}}$, $\gamma_{\rm{max}}$ and
$\gamma_{\rm{break}}$ are the Lorentz factors corresponding to the low-energy
cut off, high-energy cut off and turn over energy, respectively. The emission
region (spherical blob), which is moving down the jet with the bulk Lorentz
factor $\Gamma$, is permeated with the magnetic field B. The relativistic
electrons interact with the magnetic field and produce synchrotron radiation.
The SSC emission is produced by the same synchrotron photons generated within
the jet (Tramacere et al., 2009). The seed photons for the EC process are
external to the jet, which are produced (1) from the accretion disk (Dermer &
Schlickeiser, 2002; Ghisellini et al., 2009), (2) from broad line region
(Dermer et al., 2009), and (3) from dust torus (Sikora et al., 2008). The
external photons that directly come from the accretion disk and participate in
the EC process have luminosity Ldisk and temperature Tdisk. The reprocessed UV
photons from the BLR, which is located at a distance of RBLR from the black
hole, and IR photons from the dust torus, which is located at a distance of
RDT with a temperature of TDT, are involved in the EC process. The fraction of
Ldisk reprocessed by the BLR and dust torus are given as $\tau_{\rm{BLR}}$ and
$\tau_{\rm{DT}}$.

We have determined the size of the emission region R, which is defined as R
$\leq$ cTmin$\delta$/(1+z), where Tmin is the minimum variability timescale
and $\delta$ denotes the Doppler factor. We have calculated the minimum
variability timescale for the $\gamma$-ray light curve using the methodology
that models flares with an exponential profile, as outlined in Burbidge et al.
(1974):

$T_{\rm{min}}=\Delta t/ln(F_{2}/F_{1})$ (5)

Here, $\Delta$t = $\lvert$ t2-t1$\rvert$ and F2 and F1 are the flux values at
a time t2 and t1, respectively. The minimum variability timescale is
calculated to be $\sim$ 7.23 hours. We have adopted the Doppler factor value
for the source as $\delta$ = 18.2, which has been determined using the one-
zone leptonic emission model by Paliya et al. (2017). Using the values of Tmin
and $\delta$, we have approximated the size of emission region to be 8.24
$\times$ 1015 cm $\sim$ 8 $\times$ 1015 cm. Furthermore, assuming that this
blob covers the whole cross-section of the jets, we have estimated the
distance of this blob from the central black hole using the relation RH =
R/$\phi$, where $\phi$ represents the opening angle. The intrinsic opening
angle for the Ton 599 is estimated as 0.58∘ ($\sim$ 0.01 radian) based on
observations conducted with the 15.4 GHz Very Long Baseline Array (VLBA) from
the 2 cm VLBA MOJAVE program (Pushkarev et al., 2009). Utilizing the values of
R and $\phi$, we estimate that the location of the emission region is located
at roughly 8 $\times$ 1017 cm away from the central black hole.

We have considered a value for Ldisk equal to 4.5 $\times$ 1045 erg/sec, which
was determined through the application of a one-zone leptonic emission model
in studies by Ghisellini et al. (2010); Paliya et al. (2017). The inner radius
of the spherical shell like BLR (R${}_{\rm{BLR_{in}}}$) is considered to be
0.01 pc $\sim$ 3.1 $\times$ 1016 cm (Donea & Protheroe, 2003). The outer
radius of the BLR (R${}_{\rm{BLR_{out}}}$) is estimated to be around 2.1
$\times$ 1017 cm using the relationship RBLR = 1017 L${}^{1/2}_{\rm{disk,45}}$
cm, where Ldisk,45 represents the luminosity of the accretion disk in units of
1045 erg/sec. This estimation is based on the assumption that the BLR radius
scales in relation to the square root of the luminosity of the accretion disk
(Ghisellini & Tavecchio, 2008). Moreover, it is posited that the spherical
shell like torus is situated at a distance RDT = 2.5 $\times$ 1018
L${}^{\rm{1/2}}_{\rm{disk,45}}$ cm (Ghisellini & Tavecchio, 2008), which is
approximately equivalent to 5 $\times$ 1018 cm. We have assumed that the
fraction of the Ldisk is reprocessed by the BLR ($\tau_{\rm{BLR}}$) is 0.1 and
by the dust torus ($\tau_{\rm{DT}}$)is 0.5 (Ghisellini et al., 2010).
Furthermore, the disk temperature Tdisk and dust torus temperature TDT are
assumed to be at their typical values of 1.0 $\times$ 105 K (Patel & Chitnis,
2020) and 1.0 $\times$ 103 K (Donea & Protheroe, 2003), respectively, in our
study. And the ratio of cold protons to relativistic electrons is taken as 0.1
in our SED analysis (Ghisellini, 2012).

Table 5:  Values of the frozen parameters for the one-zone leptonic model fits
to the observed SEDs during the five epochs.

name of the parameters | values
---|---
$\gamma_{\rm{min}}$ | 40
$\gamma_{\rm{break}}$ | 1.1 $\times$ 103
$\gamma_{\rm{max}}$ | 1.0 $\times$ 104
TDT | 1000 K
RDT | 5 $\times$ 1018 cm
$\tau_{\rm{DT}}$ | 0.5
$\tau_{\rm{BLR}}$ | 0.1
R${}_{\rm{BLR_{in}}}$ | 3.1 $\times$ 1016 cm
R${}_{\rm{BLR_{out}}}$ | 2.1 $\times$ 1017 cm
Ldisk | 4.5 $\times$ 1045 erg/sec
Tdisk | 1.0 $\times$ 105 K
R | 8.0 $\times$ 1015 cm
RH | 8.0 $\times$ 1017 cm
cold proton to relativistic electron ratio | 0.1

In one-zone leptonic model, the characteristics of the observed SED are
primarily determined by nineteen parameters. Five parameters, namely, p, p1,
N, B, and $\Gamma$, are chosen to be free, while the remaining fourteen
parameters are chosen to be frozen. The purpose of parameter freezing is to
compare the changes in certain parameters, such as the magnetic field, bulk
Lorentz factor, and particle density, etc., between quiescent and flaring
epochs. Details of frozen parameters are given in Table 5. In order to
investigate the distinct correlation behavior between optical and $\gamma$-ray
flux variations during different flaring epochs, we initially fitted the
quiescent epoch to acquire the values for the parameters $\gamma_{\rm{min}}$,
$\gamma_{\rm{max}}$ and $\gamma_{\rm{break}}$. The fitting process using the
same frozen parameters and the previously acquired values of
$\gamma_{\rm{min}}$, $\gamma_{\rm{max}}$, and $\gamma_{\rm{break}}$ from the
fitting during the quiescent epoch, was repeated for the other epochs B, C, D,
and E with the selection of free parameters. Our SED analysis shows that the
leptonic scenario provides a good fit to the observed broadband SED across all
the epochs. The resultant best fit parameters for the selected epochs are
given in Table 6 and the fitted SEDs are shown in Figures 6.

$\begin{array}[]{r}\includegraphics[width=341.43306pt,height=213.39566pt]{epoch_a_new.pdf}\end{array}$

$\begin{array}[]{rr}\includegraphics[width=256.0748pt,height=213.39566pt]{epoch_b_new.pdf}\includegraphics[width=256.0748pt,height=213.39566pt]{epoch_c_new.pdf}\end{array}$

$\begin{array}[]{rr}\includegraphics[width=256.0748pt,height=213.39566pt]{epoch_d_new.pdf}\includegraphics[width=256.0748pt,height=213.39566pt]{epoch_e_new.pdf}\end{array}$

Figure 6: Observed broad band SED along with the one zone leptonic emission model fit during epochs A, B, C, D and E (top panel) and the residuals (bottom panel). Here the solid line is the best fit to the SED, orange dashed represents the emission from synchrotron process, green dashed line represents the emission from SSC process, red dashed line represents the EC emission from dusty torus and pink dashed line represents the EC emission from BLR. Purple and brown dashed lines represent the thermal emissions from the dusty torus and accretion disk, respectively. The fit to the quiescent state is applied for determining the values of the parameters $\gamma_{\rm{min}}$, $\gamma_{\rm{max}}$, and $\gamma_{\rm{break}}$ in the subsequent epochs modelled in this study. Table 6:  Results of the broadband SED fitting during the epochs A, B, C, D, and E.  Parameters | Epoch A | Epoch B | Epoch C | Epoch D | Epoch E
---|---|---|---|---|---
p | 2.29$\pm$0.31 | 2.40$\pm$0.11 | 2.04$\pm$0.06 | 2.55$\pm$0.11 | 2.17$\pm$0.13
p1 | 4.03$\pm$0.37 | 3.20$\pm$0.15 | 3.18$\pm$0.08 | 3.60$\pm$0.19 | 3.85$\pm$0.21
N ($\times$ 103 cm${}^{-3})$ | 0.98$\pm$0.20 | 1.19$\pm$0.15 | 1.31$\pm$0.09 | 0.89$\pm$0.11 | 0.84$\pm$0.11
B (Gauss) | 1.63$\pm$0.29 | 1.51$\pm$0.10 | 1.51$\pm$0.07 | 2.05$\pm$0.15 | 1.73$\pm$0.14
Bulk Lorentz factor | 22.48$\pm$2.78 | 37.31$\pm$1.75 | 28.64$\pm$0.83 | 37.37$\pm$1.87 | 38.54$\pm$2.01

##  4 Discussion

###  4.1 Connection between optical and $\gamma$-ray flux variations

Over a span of around 15 years of study, Ton 599 has shown varied correlated
behavior between optical V-band flux and $\gamma$-ray flux variations. We
found two cases, where, (1) both the optical flux and the $\gamma$-ray flux
showed flaring behavior, (2) an increase in optical flux without a
corresponding increase in $\gamma$-ray flux. We chose five epochs A, B, C, D,
and E based on the distinct behaviour between optical V-band and $\gamma$-ray
flux variations. During epoch A, the source was in its quiescent state in each
waveband. During epochs B, C and E, there was a structural correlation between
the variations in optical and $\gamma$-ray flux. However, during epoch D, the
source exhibited an optical flare without a corresponding $\gamma$-ray
response. Our investigation involving linear regression analysis on the
logarithmic values of $\gamma$-ray flux and optical flux suggests a similar
trend. The linear least square fits are shown in Figure 7, depicting the
logarithm of $\gamma$-ray flux plotted against the logarithm of optical flux.
The outcomes of the linear regression analysis are given in Table 7. Moreover,
we employed the sophisticated Discrete Cross Correlation Function to cross-
correlate the optical V-band and $\gamma$-ray light curves, aiming to
establish whether there is a correlation between variations in the optical
flux and $\gamma$-ray flux. We found that there are weakly significant
correlations between the optical and $\gamma$-ray flux variations. Given the
observed structural similarities between these two energy bands, this low
significance may result from the relatively low amount of optical data used.

A time lag of a few days within the period of October-December 2017, falling
within the timeframe of epoch B considered here, was previously identified by
Prince (2019). In the cross-correlation study of bright blazars, Liodakis et
al. (2018); Liodakis et al. (2019) found that the time lag between optical and
$\gamma$-ray bands are generally small. Moreover, de Jaeger et al. (2023)
found that the time lag distribution between optical and $\gamma$-ray bands is
consistent with zero lags. These results indicate that the optical and
$\gamma$-ray flux emission regions are co-spatial and favour leptonic emission
model of blazars, where the same relativistic electron population is
accountable for both the optical and $\gamma$-ray emission. Given the scenario
that both optical and $\gamma$-ray emissions are proportionate to the number
of emitting electrons, it is plausible that variations in the electron
population of the jet lead to the correlation between these two types of flux
variations, as proposed by D’Ammando et al. (2011).

Figure 7: $\gamma$-ray v/s Optical V-band fluxes during epochs B, C, D, and E.
Table 7: The results obtained by applying linear regression on the logarithm
of $\gamma$-ray flux and optical flux during epochs B, C, D, and E.

Epoch | Intercept | Slope | Correlation coefficient | p-value
---|---|---|---|---
B | -1.65$\pm$4.05 | 0.41$\pm$0.39 | 0.28 | 0.31
C | 2.31$\pm$2.76 | 0.81$\pm$0.26 | 0.64 | 0.008
D | -4.77$\pm$4.17 | 0.15$\pm$0.40 | 0.10 | 0.71
E | 0.71$\pm$1.53 | 0.64$\pm$0.15 | 0.77 | 0.001

###  4.2 $\gamma$-ray spectra

The $\gamma$-ray spectral shape of blazars can yield valuable information
about the inherent characteristics of the jets such as the distribution of
emitting particles. In particular, these $\gamma$-ray spectra can help us
determine if there is any absorption of $\gamma$-rays within the emission
region and can also reveal whether $\gamma$-ray emissions are attributed to a
single component of the jet or multiple components. The high-energy
$\gamma$-ray spectra of FSRQ sources are better explained by the broken power
law or log parabola model (Abdo et al., 2010a; Paliya et al., 2015; Rajput et
al., 2019; Sahakyan, 2020; Rajput et al., 2020; Rajput & Pandey, 2021). The
curvature observed in the $\gamma$-ray spectra of FSRQs may be attributed to
external factors, such as the absorption of $\gamma$-rays through photon-
photon pair production within the BLR (Poutanen & Stern, 2010), and the impact
of the Klein-Nishina effect on the inverse Compton scattering of BLR photons
by relativistic jet electrons with a curved distribution (Cerruti et al.,
2013). However, both the effects account for the IC scattering of the BLR
photons, whereas in the large number of FSRQs the $\gamma$-ray emission site
lies outside the BLR and $\gamma$-ray emission is not produced by the IC
scattering of the BLR photons (Costamante et al., 2018). Our SED analysis also
indicates that the $\gamma$-ray emission is produced by the IC scattering of
the IR photons from the dusty torus.

In our work, the LP model fits the $\gamma$-ray spectra in the flaring epochs
(and orphan optical flare epoch) adequately, whereas the PL model fits the
$\gamma$-ray spectra in the quiescent epoch well. The curvature observed in
the $\gamma$-ray spectra is probably caused by the inherent characteristics of
the electron distribution that produces emission. These electrons follow
either a cutoff energy distribution or a log parabola energy distribution.
However, during quiescent epoch the low quality of photons data available for
analysis makes it challenging for us to identify this curvature. In the LP
model fits to $\gamma$-ray spectra, the parameters $\alpha$ and $\beta$ convey
crucial information about the properties of the $\gamma$-ray spectra. Changes
in the $\alpha$ and $\beta$ values are reflected in the shape of $\gamma$-ray
spectra. We looked for the variations in the $\alpha$ and $\beta$ parameters
against $\gamma$-ray flux during epochs B, C, D and E, which are shown in the
Figure 8. We conducted an analysis to find a correlation between the $\alpha$
and $\beta$ parameters and the $\gamma$-ray flux, employing the Spearman rank
test. We considered the correlation is significant if the Spearman rank
correlation coefficient $>$ 0.5 or $<$ $-$0.5, along with a p-value $<$ 0.05.
We observed that the parameter $\alpha$ showed no significant pattern, and the
curvature parameter $\beta$ exhibited a mild negative correlation with
$\gamma$-ray flux throughout all epochs. The results of Spearman rank test are
given in Table 8. Decreasing trend of $\beta$ with increasing flux was also
noticed in the study of 3LAC (The Third Catalog of Active Galactic Nuclei
Detected by the Fermi) sources by Ackermann et al. (2015). According to Coogan
et al. (2016), the variations in the location of the $\gamma$-ray emission
area during various activity levels of the sources may be responsible for
these distinct spectral behaviors of the $\gamma$-ray spectra. Such varied
spectral behavior could also be caused by a change in the IC peak frequency,
where the external photon field is responsible for producing the $\gamma$-ray
emission in FSRQs (Rajput et al., 2020).

Table 8: Results of the Spearman rank test for the correlation of $\alpha$ and $\beta$ parameters with $\gamma$-ray flux. The plots of $\alpha$ and $\beta$ parameters v/s $\gamma$-ray flux are shown in Figure 8. | $\alpha$ | $\beta$
---|---|---
Epoch | R | p-value | R | p-value
B | -0.33 | 0.002 | -0.21 | 0.04
C | 0.06 | 0.58 | -0.34 | 0.001
D | 0.26 | 0.01 | -0.27 | 0.007
E | 0.02 | 0.84 | -0.30 | 0.003

$\begin{array}[]{rrr}\includegraphics[width=227.62204pt,height=213.39566pt]{alpha_beta_b.pdf}\includegraphics[width=227.62204pt,height=213.39566pt]{alpha_beta_c.pdf}\end{array}$

$\begin{array}[]{rr}\includegraphics[width=227.62204pt,height=213.39566pt]{alpha_beta_d.pdf}\includegraphics[width=227.62204pt,height=213.39566pt]{alpha_beta_e.pdf}\end{array}$

Figure 8: Correlation between the parameters $\alpha$ and $\beta$ and the
$\gamma$-ray flux. The blue line represents the least-square fit applied to
the data.

###  4.3 Broadband spectral energy distribution

In our study, we found instances of both optical and $\gamma$-ray flares as
well as optical flare without its corresponding $\gamma$-ray counterpart.
Liodakis et al. (2019) investigated for a correlation between optical and GeV
flux changes in 178 blazars and found that in nearly 50% of the blazars,
optical flares occurred without their $\gamma$-ray counterparts. A Few blazars
with complex behavior between changes in optical and $\gamma$-ray flux were
also studied by Chatterjee et al. (2013); Cohen et al. (2014); Rajput et al.
(2019, 2020). Hadronic processes could provide a possible explanation for
optical flux variations without $\gamma$-ray counterparts or vice versa (Mücke
& Protheroe, 2001). Furthermore, the leptonic scenario of blazars favours a
close correlation between changes in the optical and GeV flux. However,
results from Rajput et al. (2019, 2020) indicate that different correlation
behaviour between optical and $\gamma$-ray flux changes are well explained by
one-zone leptonic model scenario.

In this work, we performed broadband SED modeling using a one-zone leptonic
emission model to model the five epochs considered here. The leptonic model
provides a good fit during all the epochs with the seed photons for the EC
process from the torus. The source was in a quiescent condition during epoch
A, therefore we first fitted the broadband SED for epoch A and then looked for
changes in parameters in other epochs relative to the quiescent epoch. The
accretion disk component is clearly evident in the broadband SED during the
quiescent phase, but synchrotron radiation dominates the optical-UV emission
during epochs B, C, D, and E. The source showed both optical and $\gamma$-ray
flares during epochs B, C and E. During epochs B and E, the bulk Lorentz
factor was roughly 1.5 times higher than the quiescent epoch A. Furthermore,
there was an increase in electron density and a decrease in the magnetic field
during period B. In contrast, during period E, there was a decrease in
electron density and an increase in the magnetic field. The enhanced optical
and $\gamma$-ray fluxes in these two time intervals may be attributed to the
increase bulk Lorentz factor and changes in electron density and magnetic
field. In comparison to the quiescent period, the bulk Lorentz factor changed
slightly during epoch C, and there was also decrease in the magnetic field.
However, the electron density increased during epoch C, which might be
accountable for the optical and $\gamma$-ray flares during this epoch.

During Epoch D, when the source showed an orphan optical flare, the bulk
Lorentz factor increased almost 1.5 times, while the electron density
decreased when compared to Epoch A. Additionally, the magnetic field was at
its peak during this epoch. Consequently, the orphan optical flare can be the
result of the increased magnetic field and bulk Lorentz factor during epoch D.
According to Marscher & Gear (1985) the change in magnetic field can affect
the synchrotron emission in the jet, however the $\gamma$-ray emission will be
unperturbed by the change in the magnetic field. The presence of increased
optical polarization measurements and X-ray flare can support this argument. A
significant X-ray flare is observed during this time (see the Figure 1), which
is dominated by SSC process (the same photons emitted by the synchrotron
process are involved in the EC process), however, the optical polarization
measurements are not available at this time. From SED analysis, we draw the
conclusion that the leptonic model explains the optical-GeV variations during
the various epochs in Ton 599.

##  5 Summary

We investigated the correlation between the changes in optical V-band and GeV
flux in the FSRQ Ton 599 for the time period MJD 54686–60008. Our study
includes 1) the identification of the epochs exhibiting optical or
$\gamma$-ray flare, 2) cross-correlation analysis between optical flux and
$\gamma$-ray flux variations, 3) analysis of optical spectral variations, 4)
analysis of $\gamma$-ray spectra 5) broadband SED modeling with a one-zone
leptonic emission model for the chosen epochs. We summarize the results below

  1. 1.

In the period of around 15 years, we selected five epochs A, B, C, D, and E.
The source was in its quiescent state during epoch A. Our visual inspections
indicated that, in epochs B, C, and E, the source showed variations in both
optical and $\gamma$-ray emissions, showing structural similarities. However,
during epoch, D the source showed an optical flare, without a clear
$\gamma$-ray counterpart. The source was in its highest flaring state during
epoch E.

  2. 2.

Through cross-correlation analysis, we identified a weak linear correlation
between optical and $\gamma$-ray light curves.

  3. 3.

We found no optical spectral variations during the epochs.

  4. 4.

During the quiescent period, the $\gamma$-ray spectra of the source was well
described by the PL model, while during the flaring periods, the $\gamma$-ray
spectra was best described by the LP model.

  5. 5.

During all the epochs, from SED modelling, we found the $\gamma$-ray emission
to be produced by EC process, where the source of seed photons is the torus.
Our SED analysis showed that during the quiescent period, the optical-UV
emission was caused by the accretion disk component, but during the flaring
periods, the synchrotron emission was the cause for the optical-UV emission.
The varied behaviour observed in our work between the changes in optical and
GeV flux, could be explained by changes in the magnetic field, the bulk
Lorentz factor, and the electron density.

## Acknowledgements

We are thankful to the anonymous referee for providing valuable feedback,
which helped to improve the manuscript. B.R. gratefully acknowledges support
by ANID BASAL project FB210003 and FONDECYT Postdoctorado 3230631. W.M.
gratefully acknowledges support by the ANID BASAL project FB210003. Part of
this work was supported by the Polish Funding Agency National Science Centre,
project 2017/26/A/ST9/00756 (MAESTRO 9). This research has made use of data
from Fermi Science Support Center, which was made available by NASA’s Goddard
Space Flight Center (GSFC). The High Energy Astrophysics Science Archive
Research Center (HEASARC), a service of GSFC, is used to obtain the data,
software, and online tools. This study used NuSTAR data as well as NuSTAR Data
Analysis Software (NuSTARDAS), which was collaboratively developed by the
California Institute of Technology (Caltech, USA) and the ASI Science Data
Center (ASDC, Italy). We are thankful for photometry program from Paul Smiths
monitoring program at the Steward Observatory, which is supported by Fermi
Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G.
This study utilized ASAS-SN data, which was partially funded by the Alfred P.
Sloan Foundation grant G-2021-14192 and Gordon and Betty Moore Foundation
grants GBMF5490 and GBMF10501 to the Ohio State University. This research
utilized the publicly accessible numerical code
JetSet888https://github.com/andreatramacere/jetset, developed by Tramacera A.

## Data Availability

The multi-wavelength data utilized in this study are publicly accessible from
the Fermi-LAT999https://fermi.gsfc.nasa.gov/ssc/data/access/lat/, Swift-
XRT101010https://heasarc.gsfc.nasa.gov/docs/archive.html,
NuSTAR111111http://heasarc.gsfc.nasa.gov/docs/nustar/analysis/, Swift-
UVOT121212https://www.ssdc.asi.it/cgi-bin/swiftuvarchint, Steward
observatory131313http://james.as.arizona.edu/$\sim$psmith/Fermi/ and ASAS-
SN141414https://asas-sn.osu.edu/.

## References

  * Abdo et al. (2010a) Abdo A. A., et al., 2010a, ApJ, 710, 1271
  * Abdo et al. (2010b) Abdo A. A., et al., 2010b, ApJ, 715, 429
  * Abdo et al. (2010c) Abdo A. A., et al., 2010c, ApJ, 716, 30
  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Abdollahi et al. (2022) Abdollahi S., et al., 2022, ApJS, 260, 53
  * Ackermann et al. (2015) Ackermann M., et al., 2015, ApJ, 810, 14
  * Acosta Pulido et al. (2017) Acosta Pulido J. A., Redondo J., Carnerero M. I., Raiteri C. M., Villata M., Larionov V. M., Kopatskaya E. N., GASP-WEBT Collaboration 2017, The Astronomer’s Telegram, 10949, 1
  * Aharonian (2000) Aharonian F. A., 2000, New Astron., 5, 377
  * Antonucci (1993) Antonucci R., 1993, ARA&A, 31, 473
  * Arnaud (1996) Arnaud K. A., 1996, in Jacoby G. H., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 101, Astronomical Data Analysis Software and Systems V. p. 17
  * Atwood et al. (2009) Atwood W. B., et al., 2009, ApJ, 697, 1071
  * Bhatta et al. (2023) Bhatta G., et al., 2023, MNRAS,
  * Blandford & Rees (1978) Blandford R. D., Rees M. J., 1978, in Wolfe A. M., ed., BL Lac Objects. pp 328–341
  * Błażejowski et al. (2000) Błażejowski M., Sikora M., Moderski R., Madejski G. M., 2000, ApJ, 545, 107
  * Boettcher et al. (1997) Boettcher M., Mause H., Schlickeiser R., 1997, A&A, 324, 395
  * Bonning et al. (2009) Bonning E. W., et al., 2009, ApJ, 697, L81
  * Bonning et al. (2012) Bonning E., et al., 2012, ApJ, 756, 13
  * Böttcher (2007) Böttcher M., 2007, Ap&SS, 309, 95
  * Böttcher et al. (2013) Böttcher M., Reimer A., Sweeney K., Prakash A., 2013, ApJ, 768, 54
  * Breeveld et al. (2011) Breeveld A. A., Landsman W., Holland S. T., Roming P., Kuin N. P. M., Page M. J., 2011, in McEnery J. E., Racusin J. L., Gehrels N., eds, American Institute of Physics Conference Series Vol. 1358, Gamma Ray Bursts 2010. pp 373–376 (arXiv:1102.4717), doi:10.1063/1.3621807
  * Burbidge et al. (1974) Burbidge G. R., Jones T. W., O’Dell S. L., 1974, ApJ, 193, 43
  * Burrows et al. (2005) Burrows D. N., et al., 2005, Space Sci. Rev., 120, 165
  * Carrasco et al. (2017) Carrasco L., Escobedo G., Recillas E., Porras A., Chavushyan V., Mayya D. Y., 2017, The Astronomer’s Telegram, 10932, 1
  * Cerruti et al. (2013) Cerruti M., Dermer C. D., Lott B., Boisson C., Zech A., 2013, ApJ, 771, L4
  * Cerruti et al. (2019) Cerruti M., Zech A., Boisson C., Emery G., Inoue S., Lenain J. P., 2019, MNRAS, 483, L12
  * Chatterjee et al. (2012) Chatterjee R., et al., 2012, ApJ, 749, 191
  * Chatterjee et al. (2013) Chatterjee R., et al., 2013, ApJ, 763, L11
  * Cheung et al. (2017) Cheung C. C., Gasparrini D., Buson S., 2017, The Astronomer’s Telegram, 10931, 1
  * Cohen et al. (2014) Cohen D. P., Romani R. W., Filippenko A. V., Cenko S. B., Lott B., Zheng W., Li W., 2014, ApJ, 797, 137
  * Coogan et al. (2016) Coogan R. T., Brown A. M., Chadwick P. M., 2016, MNRAS, 458, 354
  * Costamante et al. (2018) Costamante L., Cutini S., Tosti G., Antolini E., Tramacere A., 2018, MNRAS, 477, 4749
  * D’Ammando et al. (2011) D’Ammando F., et al., 2011, A&A, 529, A145
  * Dermer & Schlickeiser (2002) Dermer C. D., Schlickeiser R., 2002, ApJ, 575, 667
  * Dermer et al. (2009) Dermer C. D., Finke J. D., Krug H., Böttcher M., 2009, ApJ, 692, 32
  * Dickey & Lockman (1990) Dickey J. M., Lockman F. J., 1990, ARA&A, 28, 215
  * Diltz & Böttcher (2016) Diltz C., Böttcher M., 2016, ApJ, 826, 54
  * Donea & Protheroe (2003) Donea A.-C., Protheroe R. J., 2003, Astroparticle Physics, 18, 377
  * Dutka et al. (2013) Dutka M. S., et al., 2013, ApJ, 779, 174
  * Edelson & Krolik (1988) Edelson R. A., Krolik J. H., 1988, ApJ, 333, 646
  * Fan et al. (2006) Fan J. H., et al., 2006, PASJ, 58, 797
  * Foffano et al. (2019) Foffano L., Prandini E., Franceschini A., Paiano S., 2019, MNRAS, 486, 1741
  * Fossati et al. (1998) Fossati G., Maraschi L., Celotti A., Comastri A., Ghisellini G., 1998, MNRAS, 299, 433
  * Garrappa & Valverd (2023) Garrappa S., Valverd J., 2023, The Astronomer’s Telegram, 15859, 1
  * Gehrels et al. (2004) Gehrels N., et al., 2004, ApJ, 611, 1005
  * Ghisellini (2012) Ghisellini G., 2012, MNRAS, 424, L26
  * Ghisellini & Madau (1996) Ghisellini G., Madau P., 1996, MNRAS, 280, 67
  * Ghisellini & Maraschi (1989) Ghisellini G., Maraschi L., 1989, ApJ, 340, 181
  * Ghisellini & Tavecchio (2008) Ghisellini G., Tavecchio F., 2008, MNRAS, 387, 1669
  * Ghisellini et al. (2009) Ghisellini G., Tavecchio F., Ghirlanda G., 2009, MNRAS, 399, 2041
  * Ghisellini et al. (2010) Ghisellini G., Tavecchio F., Foschini L., Ghirlanda G., Maraschi L., Celotti A., 2010, MNRAS, 402, 497
  * Ghisellini et al. (2011) Ghisellini G., Tavecchio F., Foschini L., Ghirland a G., 2011, MNRAS, 414, 2674
  * Goyal (2021) Goyal A., 2021, ApJ, 909, 39
  * Gu et al. (2006) Gu M. F., Lee C. U., Pak S., Yim H. S., Fletcher A. B., 2006, A&A, 450, 39
  * Harrison et al. (2013) Harrison F. A., et al., 2013, ApJ, 770, 103
  * Hewett & Wild (2010) Hewett P. C., Wild V., 2010, MNRAS, 405, 2302
  * Impey & Neugebauer (1988) Impey C. D., Neugebauer G., 1988, AJ, 95, 307
  * Kalberla et al. (2005) Kalberla P. M. W., Burton W. B., Hartmann D., Arnal E. M., Bajaja E., Morras R., Pöppel W. G. L., 2005, A&A, 440, 775
  * Kochanek et al. (2017) Kochanek C. S., et al., 2017, PASP, 129, 104502
  * Konigl (1981) Konigl A., 1981, ApJ, 243, 700
  * Liao et al. (2014) Liao N. H., Bai J. M., Liu H. T., Weng S. S., Chen L., Li F., 2014, ApJ, 783, 83
  * Liodakis et al. (2018) Liodakis I., Romani R. W., Filippenko A. V., Kiehlmann S., Max-Moerbeck W., Readhead A. C. S., Zheng W., 2018, MNRAS, 480, 5517
  * Liodakis et al. (2019) Liodakis I., Romani R. W., Filippenko A. V., Kocevski D., Zheng W., 2019, ApJ, 880, 32
  * MacDonald et al. (2015) MacDonald N. R., Marscher A. P., Jorstad S. G., Joshi M., 2015, ApJ, 804, 111
  * Mannheim (1993) Mannheim K., 1993, A&A, 269, 67
  * Mao et al. (2016) Mao P., Urry C. M., Massaro F., Paggi A., Cauteruccio J., Künzel S. R., 2016, ApJS, 224, 26
  * Marscher & Gear (1985) Marscher A. P., Gear W. K., 1985, ApJ, 298, 114
  * Mattox et al. (1996) Mattox J. R., et al., 1996, ApJ, 461, 396
  * Max-Moerbeck et al. (2014a) Max-Moerbeck W., et al., 2014a, MNRAS, 445, 428
  * Max-Moerbeck et al. (2014b) Max-Moerbeck W., Richards J. L., Hovatta T., Pavlidou V., Pearson T. J., Readhead A. C. S., 2014b, MNRAS, 445, 437
  * 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
  * Mukherjee & VERITAS Collaboration (2017) Mukherjee R., VERITAS Collaboration 2017, The Astronomer’s Telegram, 11075, 1
  * Negi et al. (2022) Negi V., Joshi R., Chand K., Chand H., Wiita P., Ho L. C., Singh R. S., 2022, MNRAS, 510, 1791
  * Nolan et al. (2012) Nolan P. L., et al., 2012, ApJS, 199, 31
  * Padovani & Giommi (1995) Padovani P., Giommi P., 1995, ApJ, 444, 567
  * Paliya et al. (2015) Paliya V. S., Sahayanathan S., Stalin C. S., 2015, ApJ, 803, 15
  * Paliya et al. (2016) Paliya V. S., Diltz C., Böttcher M., Stalin C. S., Buckley D., 2016, ApJ, 817, 61
  * Paliya et al. (2017) Paliya V. S., Marcotulli L., Ajello M., Joshi M., Sahayanathan S., Rao A. R., Hartmann D., 2017, ApJ, 851, 33
  * Patel & Chitnis (2020) Patel S. R., Chitnis V. R., 2020, MNRAS, 492, 72
  * Patel et al. (2018) Patel S. R., Chitnis V. R., Shukla A., Rao A. R., Nagare B. J., 2018, ApJ, 866, 102
  * Poutanen & Stern (2010) Poutanen J., Stern B., 2010, ApJ, 717, L118
  * Prince (2019) Prince R., 2019, ApJ, 871, 101
  * Prince (2023) Prince R., 2023, The Astronomer’s Telegram, 15854, 1
  * Pushkarev et al. (2009) Pushkarev A. B., Kovalev Y. Y., Lister M. L., Savolainen T., 2009, A&A, 507, L33
  * Rajput & Pandey (2021) Rajput B., Pandey A., 2021, Galaxies, 9, 118
  * Rajput et al. (2019) Rajput B., Stalin C. S., Sahayanathan S., Rakshit S., Mandal A. K., 2019, MNRAS, 486, 1781
  * Rajput et al. (2020) Rajput B., Stalin C. S., Sahayanathan S., 2020, MNRAS, 498, 5128
  * Safna et al. (2020) Safna P. Z., Stalin C. S., Rakshit S., Mathew B., 2020, MNRAS, 498, 3578
  * Sahakyan (2020) Sahakyan N., 2020, A&A, 635, A25
  * Sbarrato et al. (2012) Sbarrato T., Ghisellini G., Maraschi L., Colpi M., 2012, MNRAS, 421, 1764
  * Scarpa & Falomo (1997) Scarpa R., Falomo R., 1997, A&A, 325, 109
  * Shah et al. (2017) Shah Z., Sahayanathan S., Mankuzhiyil N., Kushwaha P., Misra R., Iqbal N., 2017, MNRAS, 470, 3283
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shappee et al. (2014) Shappee B. J., et al., 2014, ApJ, 788, 48
  * Shaw et al. (2012) Shaw M. S., et al., 2012, ApJ, 748, 49
  * Sikora et al. (2008) Sikora M., Moderski R., Madejski G. M., 2008, ApJ, 675, 71
  * Stocke et al. (1991) Stocke J. T., Morris S. L., Gioia I. M., Maccacaro T., Schild R., Wolter A., Fleming T. A., Henry J. P., 1991, ApJS, 76, 813
  * Thompson et al. (1995) Thompson D. J., et al., 1995, ApJS, 101, 259
  * Tramacere (2020) Tramacere A., 2020, JetSeT: Numerical modeling and SED fitting tool for relativistic jets, Astrophysics Source Code Library, record ascl:2009.001 (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
  * Tripathi et al. (2023) Tripathi T., Devanand P. U., Gupta A. C., Kishore S., Dogra K., Krishna Mohana A., Dhiman V., 2023, The Astronomer’s Telegram, 15875, 1
  * Urry & Mushotzky (1982) Urry C. M., Mushotzky R. F., 1982, ApJ, 253, 38
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, PASP, 107, 803
  * Urry et al. (2000) Urry C. M., Scarpa R., O’Dowd M., Falomo R., Pesce J. E., Treves A., 2000, ApJ, 532, 816
  * Vermeulen et al. (1995) Vermeulen R. C., Ogle P. M., Tran H. D., Browne I. W. A., Cohen M. H., Readhead A. C. S., Taylor G. B., Goodrich R. W., 1995, ApJ, 452, L5
  * Wu et al. (2011) Wu J., Zhou X., Ma J., Jiang Z., 2011, MNRAS, 418, 1640
  * Zhang et al. (2015) Zhang B.-K., Zhou X.-S., Zhao X.-Y., Dai B.-Z., 2015, Research in Astronomy and Astrophysics, 15, 1784
  * de Jaeger et al. (2023) de Jaeger T., et al., 2023, MNRAS, 519, 6349