# Multi-wavelength temporal variability of the blazar PKS 1510$-$089

Q. Yuan Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 150
Science 1-Street, Urumqi 830011, China  University of Chinese Academy of
Sciences, 19A Yuquan Road, Beijing 100049, China  Pankaj Kushwaha DST-INSPIRE
Faculty  Department of Physical Sciences, Indian Institute of Science
Education and Research Mohali, Knowledge City, Sector 81, SAS Nagar, Punjab
140306, India  Aryabhatta Research Institute of Observational Sciences
(ARIES), Manora Peak, Nainital 263 001, India  Alok C. Gupta Key Laboratory
for Research in Galaxies and Cosmology, Shanghai Astronomical Observatory,
Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China
Aryabhatta Research Institute of Observational Sciences (ARIES), Manora Peak,
Nainital 263 001, India  Ashutosh Tripathi Department of Physics, Southern
Methodist University, 3215 Daniel Ave., Dallas, TX 75205, USA  Center for
Field Theory and Particle Physics and Department of Physics, Fudan University,
2005 Songhu Road, Shanghai 200438, China  Paul J. Wiita Department of Physics,
The College of New Jersey, 2000 Pennington Rd., Ewing, NJ 08628-0718, USA  M.
Zhang Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 150
Science 1-Street, Urumqi 830011, China  Key Laboratory for Radio Astronomy,
Chinese Academy of Sciences, 2 West Beijing Road, Nanjing 210008, China
Xinjiang Key Laboratory of Radio Astrophysics, 150 Science 1-Street, Urumqi
830011, China  X. Liu  Xinjiang Astronomical Observatory, Chinese Academy of
Sciences, 150 Science 1-Street, Urumqi 830011, China  Key Laboratory for Radio
Astronomy, Chinese Academy of Sciences, 2 West Beijing Road, Nanjing 210008,
China  Xinjiang Key Laboratory of Radio Astrophysics, 150 Science 1-Street,
Urumqi 830011, China  Anne Lähteenmäki Aalto University Metsähovi Radio
Observatory, Metsähovintie 114, 02540 Kylmälä, Finland  Aalto University
Department of Electronics and Nanoengineering, P.O. BOX 15500, FI-00076 Aalto,
Finland  Merja Tornikoski Aalto University Metsähovi Radio Observatory,
Metsähovintie 114, 02540 Kylmälä, Finland  Joni Tammi Aalto University
Metsähovi Radio Observatory, Metsähovintie 114, 02540 Kylmälä, Finland
Venkatessh Ramakrishnan Aalto University Metsähovi Radio Observatory,
Metsähovintie 114, 02540 Kylmälä, Finland  Finnish Centre for Astronomy with
ESO (FINCA), University of Turku, Vesilinnantie 5, 20014 University of Turku,
Finland  L. Cui  Xinjiang Astronomical Observatory, Chinese Academy of
Sciences, 150 Science 1-Street, Urumqi 830011, China  Key Laboratory for Radio
Astronomy, Chinese Academy of Sciences, 2 West Beijing Road, Nanjing 210008,
China  Xinjiang Key Laboratory of Radio Astrophysics, 150 Science 1-Street,
Urumqi 830011, China  X. Wang  Xinjiang Astronomical Observatory, Chinese
Academy of Sciences, 150 Science 1-Street, Urumqi 830011, China  University of
Chinese Academy of Sciences, 19A Yuquan Road, Beijing 100049, China  M.F. Gu
Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical
Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030,
China  Cosimo Bambi Center for Field Theory and Particle Physics and
Department of Physics, Fudan University, 2005 Songhu Road, Shanghai 200438,
China  A. E. Volvach Radio Astronomy and Geodynamics Department CrAO, Taras
Shevchenko National University of Kyiv, Ukraine

###### Abstract

We perform correlation and periodicity search analyses on long-term multi-band
light curves of the FSRQ 1510$-$089 observed by the space-based Fermi–Large
Area Telescope in $\gamma$-rays, the SMARTS and Steward Observatory telescopes
in optical and near-infrared (NIR) and the 13.7 m radio telescope in Metsähovi
Radio Observatory between 2008 and 2018. The z-transform discrete correlation
function method is applied to study the correlation and possible time lags
among these multi-band light curves. Among all pairs of wavelengths, the
$\gamma$-ray vs. optical/NIR and optical vs. NIR correlations show zero time
lags; however, both the $\gamma$-ray and optical/NIR emissions precede the
radio radiation. The Generalized Lomb-Scargle periodogram, Weighted Wavelet
Z-transform, and REDFIT techniques are employed to investigate the unresolved-
core-emission dominated 37 GHz light curve and yield evidence for a quasi-
period around 1540 days, although given the length of the whole data set it
cannot be claimed to be significant. We also investigate the optical/NIR color
variability and find that this source shows a simple redder-when-brighter
behavior over time, even in the low flux state.

Active galactic nuclei (16); Blazars (164); Flat-spectrum radio quasars
(2163); Observational astronomy (1145)

††software: Python

##  1 Introduction

Blazars are radio-loud (RL) active galactic nuclei (AGNs) whose relativistic
jets are seen at a small angle to the line of sight (Urry & Padovani, 1995).
Traditionally, blazars are subclassified as flat spectrum radio quasars
(FSRQs) and BL Lac objects (BL Lacs). There are strong emission lines in the
spectrum of the former while these lines are weak or absent in the latter. The
fluxes from blazars are highly variable in the entire accessible
electromagnetic spectrum, with time scales ranging from minutes to years
(e.g., Urry et al., 1993; Petry et al., 2000; Katarzyński et al., 2001;
Krawczynski et al., 2001; Aharonian et al., 2002, 2005a; Błażejowski et al.,
2005; Aharonian et al., 2005b; Rebillot et al., 2006; Aharonian et al., 2007;
Fossati et al., 2008; Horan et al., 2009; Aharonian et al., 2009; Gupta et
al., 2017, 2019, 2022, and references therein).

The Spectral Energy Distributions (SEDs) of blazars are characterized by a
broad double-peaked structure. The low-energy hump spans radio to ultraviolet
(UV) or X-ray bands with a peak between near-infrared (NIR) and UV/X-ray. This
low-energy hump is synchrotron radiation from the relativistic non-thermal
electrons in the jet. The high-energy hump, on the other hand, extends from
X-rays to GeV/TeV $\gamma$-rays, peaking between hard X-rays and
$\gamma$-rays. The origin of the second hump is still unclear; both lepton and
hadron based emission scenarios have been proposed (e.g., Böttcher et al.,
2013). In leptonic scenarios, the mechanism of the high energy emission is the
inverse-Compton scattering of low-energy seed photons by the same relativistic
electrons that produces the synchrotron radiation. Those low-energy seed
photons may be the synchrotron photons produced in the jet (synchrotron self-
Compton, SSC) or be external photons originating in the local environment
(external Compton, EC) which includes the accretion disk, the broad line
region (BLR) and the dusty torus (DT) or conceivably even the cosmic microwave
background (e.g., Tavecchio et al., 2000; Kushwaha et al., 2013, and
references therein). In the hadronic scenario, the second hump is due to the
proton synchrotron or the proton-induced particle cascades (e.g., Mannheim &
Biermann, 1992; Mücke & Protheroe, 2001).

Though broadly stochastic in nature, flux variability in the various wavebands
differs from source to source. Even for the same source, there can be
significant differences in the correlation of light curves between multiple
bands during different observation periods, so studying multi-band spectral
and temporal variations offer a non-direct way to peek into unresolved
processes and spatial scales hidden in compact systems. In general, the
emission variability of blazars may caused by a combination of factors
internal to the jet, such as shocks, turbulence, or magnetic reconnection, and
external ones, such as changes in the jet direction with respect to the
observer (e.g. Böttcher, 2019). Blazar variability patterns pose challenges to
some existing blazar radiation models and can help discriminate between them.

The location of the $\gamma$-ray emission site(s) in blazars has been the
subject of substantial debate (e.g. Kramarenko et al., 2022; Acharyya et al.,
2021; Arsioli & Chang, 2018; Costamante et al., 2018; Hodgson et al., 2017;
Kushwaha et al., 2013, and references therein). It is well known that the
optical emission of blazars is usually dominated by relativistic jets, but
there are exceptions. For example, Fernandes et al. (2020) found that the
optical emission in the FRSQ 3C 273 was dominated by the accretion disk rather
than the jet over the entire time span of their study. A detailed multi-
wavelength cross-correlation study can analyze the radiation mechanisms in
different bands and provide constraints on the location of the $\gamma$-ray
emission region(s).

PKS 1510$-$089 is a bright FSRQ located at redshift, $z=0.36$ (Thompson et
al., 1990; Tanner et al., 1996). It is one of the nine FSRQs detected in the
very-high-energy (VHE $\geq$ 100 Gev) range by the Imaging Atmospheric
Cherenkov Telescopes (IACTs)111http://tevcat.uchicago.edu and also one of the
persistent MeV-GeV emitters (e.g. Kushwaha et al., 2017b). The High Energy
Spectroscopic System (H.E.S.S.) collaboration discovered VHE $\gamma-$ray
emission $>$ 0.1 TeV from PKS 1510$-$089 in its 15.8 hours of observation (at
9.2$\sigma$)(H. E. S. S. Collaboration et al., 2013). The source emitted
persistent VHE $\gamma-$ray emission in MAGIC observations during its low flux
state for an extended period of time during 2012 – 2017 (MAGIC Collaboration
et al., 2018). PKS 1510$-$089 has witnessed several simultaneous multi-
wavelength observing campaigns that obtained flux and spectral variabilities,
performed SED modeling, and discussed the diverse emission mechanisms at
different epochs of observations (e.g., Ahnen et al., 2017; Aleksić et al.,
2014; Orienti et al., 2013a; Chen et al., 2012; Abdo et al., 2010; Marscher et
al., 2010; Kataoka et al., 2008a, and references therein). Recently, transient
$\gamma-$ray quasi-periodic oscillations (QPOs) with periods of 3.6 days and
92 days were reported in its Fermi-LAT observations covering $\sim$ 12 years
(Roy et al., 2022a).

PKS 1510$-$089 has had a substantial amount of well-synchronized sampled data
taken in optical and NIR bands. It is also a target source for numerous blazar
radio monitoring programs, such as the Fermi-GST AGN Multi-frequency
Monitoring Alliance (F-GAMMA) program (Fuhrmann et al., 2016; Angelakis et
al., 2019), the ongoing Owens Valley Radio Observatory (OVRO) 40m monitoring
program (Richards et al., 2011), and the Metsähovi Radio Observatory blazar
monitoring program. In addition to these monitoring campaigns, there are quite
a few higher-resolution very-long-baseline interferometry (VLBI) observations:
the VLBA (Very Long Baseline Array) 2 cm Survey (Kellermann et al., 1998) and
its successor – MOJAVE (Monitoring Of Jets in Active galactic nuclei with VLBA
Experiments) program (Lister et al., 2018); the VLBA-BU-BLAZAR program and its
successor – BEAM-ME (Blazars Entering the Astrophysical Multi-Messenger Era)
program (Jorstad & Marscher, 2016; Jorstad et al., 2017; Weaver et al., 2022).
The correlations between different bands have already been studied for PKS
1510$-$089\. Abdo et al. (2010) reported that the $\gamma$-ray light curve had
a complicated correlation with other wavelengths during the period between
September 2008 and June 2009. They saw no correlation with the X-ray band, a
weak correlation with the UV band, and a significant correlation with the
optical band, with the $\gamma$-ray light curve preceding the optical one by
about 13 days. In observations taken during 2009–2013,Beaklini et al. (2017)
detected a correlation between the radio and $\gamma$-ray flares where the
radio flares are delayed by approximately 54 days relative to the $\gamma$-ray
flares. In another study based on observations taken during July 2012 to
October 2014, the optical and $\gamma$-ray bands showed correlations with zero
time delay (Ramakrishnan et al., 2016).

Broadband SED studies using multi-wavelength data at different source activity
states e.g. August 2006 (Kataoka et al., 2008b), August 2008 to May 2012
(Brown, 2013), and very high energy $\gamma-$ray flares in March 2009
(Barnacka et al., 2014) indicate that the radio and optical/NIR emission is
synchrotron radiation from the jet, while the high energy emission (X-ray and
$\gamma$-ray) arises from the EC scattering of seed photos from the BLR and
DT.

This paper focuses on multi-band emission connections over lengthy timescales,
as we try to better understand blazar radiation mechanisms and jet kinematic
behaviors. Details of the multi-band data we include and their reduction
procedures for a ten-year-long observation period are in Section 2, with a
study of spectral variations given in Section 2.5. In section 3, we present
the correlation analysis method and report results among different bands. The
subsequent section deals with the methods used to search for QPOs and the
results. We discuss these results in Section 5 and summarize our findings in
Section 6.

##  2 MULTI-WAVELENGTH ARCHIVAL DATA AND REDUCTION

###  2.1 Fermi $\gamma$-ray data

The Fermi mission is a space-based $\gamma$-ray observatory launched in 2008
(Abdo et al., 2010). It is sensitive to gamma-ray photons with energies from
$\gtrsim 20$ MeV to $\sim$ TeV. It continuously scans the sky for
$\gamma$-rays and covers the entire sky region within 90 minutes.

We analyzed the $\gamma$-ray data from the Large Area Telescope (LAT)
following the standard analysis procedures employing a “binned likelihood”
analysis with the Fermipy (v1.0.1) software. We downloaded the PASS8 (P8R3)
instrument response function to process data from Fermi-LAT. First, we applied
the intended cuts as recommended by selecting only the events marked as
“SOURCE” class (evclass=128 and evtype=3) with energies above 100 MeV from a
region of interest of $15^{\circ}$ centered on the source location
($\alpha_{2000.0}=$ 15h 12m 52.2s, $\delta_{2000.0}=-09^{\circ}\ 06^{\prime}\
21.6^{\prime\prime}$). A zenith angle cut of $90^{\circ}$ was also applied as
recommended by the LAT team. The good time intervals (GTIs) were generated
using the standard criteria ”$\rm(DATA\\_QUAL\>0)\&\&(LAT\\_CONFIG==1)$”. For
spectral modeling, an XML file with sources from the fourth LAT point source
catalog (4FGL; Abdollahi et al., 2020) was used. Apart from the $\gamma$-ray
emitting point sources, the XML model file also included the galactic and
extragalactic isotropic contributions through the template file
“1gll_psc_v21.fits” and “iso_P8R3_SOURCE_V2_v1.txt”, provided by the LAT team.
The source was modeled with a log-parabola spectrum, as in the catalog. A
test-statistic (TS) value of $\geq 9$ was used to select the data for our
analysis. The light curve (photon flux vs time) extracted following this
procedure is shown in the top panel of Fig. 1.

Figure 1: Multiband light curves of PKS 1510$-$089 from $\gamma$-ray to radio.
Vertical lines divide the segments considered for temporal analysis in our
work, the purple-shaded regions represent the $\gamma$-ray flaring states. The
green horizontal dotted line indicates the constant flux value that divides
the high and low $\gamma$-ray states.

###  2.2 Optical and NIR data

The optical R-band data and NIR J-band data for PKS 1510$-$089 are taken from
the public archives of the Small and Moderate Aperture Research Telescope
System (SMARTS)222http://www.astro.yale.edu/smarts/glast/home.php and the
Steward Observatory telescopes. SMARTS consists of four meter-class telescopes
at the Cerro Tololo Inter-American Observatory (CTIO) in Chile: the 0.9 m, 1.0
m, 1.3 m, and 1.5 m, which observed the Fermi-LAT monitored blazars at both
optical and NIR wavelengths. Details about the SMARTS telescopes, detectors,
observations, and data analysis are given in Bonning et al. (2012) and Buxton
et al. (2012).

The Steward Observatory of the University of Arizona uses the 2.3 m Bok and
1.54 m Kuiper telescopes to carry out optical photometric and polarimetric
observations of a large number of blazars using the spectropolarimeter (SPOL).
Details about these telescopes, instruments, observations, and data analysis
methods are given in Smith et al. (2009). These R-band observations are
combined with those of SMARTS in the second panel of Fig. 1 and the SMARTS
J-band data are in the third panel. The Galactic dust reddening and extinction
for a line of sight must be accounted for. The total Galactic visual
extinction is estimated and corrected using the extinction map given by
Schlegel et al. (1998). Newer estimates of Galactic dust extinction from
Schlafly & Finkbeiner (2011) are now provided alongside those of Schlegel et
al. (1998). The magnitudes are corrected for galactic extinction using an
online tool333https://ned.ipac.caltech.edu/extinction_calculator based on
Schlafly & Finkbeiner (2011).

###  2.3 Radio Data

The 37 GHz observations of the blazar PKS 1510$-$089 were made with the 13.7 m
diameter Aalto University Metsähovi radio telescope, which is a radome
enclosed Cassegrain type antenna situated in Finland (60∘ 13${}^{{}^{\prime}}$
04${}^{{}^{\prime\prime}}$ N, 24∘ 23${}^{{}^{\prime}}$
35${}^{{}^{\prime\prime}}$ E). The receivers have HEMPT (high electron
mobility pseudomorphic transistor) front ends operating at room temperature.
The observations are Dicke switched ON–ON observations, alternating the source
and the sky in each feed horn. A typical integration time to obtain one flux
density data point is between 1200 and 1800 s. The detection limit of the
telescope at 37 GHz is on the order of 0.2 Jy under optimal conditions, but is
heavily weather dependent. Data points with a signal-to-noise ratio $<$ 4 are
handled as non-detections. The flux density scale is set by observations of
the H II region DR 21. Sources NGC 7027, 3C 274 and 3C 84 are used as
secondary calibrators. A detailed description of the data reduction and
analysis is given in Teraesranta et al. (1998). The error estimate in the flux
density includes the contribution from the measurement $RMS$ and the
uncertainty of the absolute calibration. Radio data at 37 GHz are plotted in
the bottom panel of Fig. 1.

###  2.4 Multi-wavelength light curves

Fig.1 shows the multi-wavelength light curves of PKS 1510$-$089 between 2008
and 2018. The aggregated data show frequent changes in the fluxes in different
electromagnetic (EM) bands, with the high (active or flaring) states of the
$\gamma$-ray fluxes shaded in purple (Fig. 1); that shading is extended to the
other bands for convenient visualization. To examine the apparent correlations
between variations at different wavelengths in detail, we divide the complete
data into four individual segments, denoted by vertical dashed lines, each of
which has at least one high $\gamma$-ray state and concludes in a gap in the
optical/NIR corresponding to the end of an observing season. These segments
are chosen to be long enough to allow for detection of possible changes in
correlations between different bands, as seen in our earlier work (e.g., Gaur
et al., 2014; Gupta et al., 2017; Kushwaha et al., 2017a, 2018a, 2018b).

Table 1: Start and Stop dates of the four individual segments SEGMENT | Start Date | Stop Date | Start MJD | Stop MJD
---|---|---|---|---
SEGMENT 1 | 2008 Jan 28 | 2009 Oct 11 | 54493 | 55115
SEGMENT 2 | 2009 Oct 18 | 2012 Oct 02 | 55122 | 56202
SEGMENT 3 | 2012 Oct 02 | 2014 Dec 01 | 56202 | 56992
SEGMENT 4 | 2014 Dec 01 | 2017 Sep 18 | 56992 | 58014
Table 2: Each sub-segment is defined as the start and the end times of the high state on $\gamma$-ray emission |  | Start Date | Stop Date | Start MJD | Stop MJD
---|---|---|---|---|---
| SEGMENT 1-1 | 2008 Sep 18 | 2008 Sep 28 | 54727 | 54737
SEGMENT 1 | SEGMENT 1-2 | 2009 Jan 16 | 2009 Feb 16 | 54847 | 54857
| SEGMENT 1-3 | 2009 Mar 17 | 2009 Apr 16 | 54907 | 54937
| SEGMENT 1-4 | 2009 Apr 26 | 2009 Jul 15 | 54947 | 55027
| SEGMENT 2-1 | 2009 Dec 12 | 2010 Jan 21 | 55177 | 55217
SEGMENT 2 | SEGMENT 2-2 | 2011 Jul 5 | 2011 Aug 24 | 55747 | 55797
| SEGMENT 2-3 | 2011 Oct 13 | 2011 Nov 22 | 55847 | 55887
| SEGMENT 2-4 | 2012 Jan 31 | 2012 Mar 31 | 55957 | 56017
SEGMENT 3 | SEGMENT 3-1 | 2013 Sep 2 | 2013 Nov 11 | 56537 | 56607
SEGMENT 4 | SEGMENT 4-1 | 2015 Mar 16 | 2015 Aug 23 | 57097 | 57257
| SEGMENT 4-2 | 2016 Sep 6 | 2016 Sep 26 | 57637 | 57657

The start and stop dates of the individual segments are given in Table 1. The
flaring states’ start and end times within each segment define sub-segments
that are listed in Table 2. There are a variety of ways to decide which data
points belong to the high and low (quiescent) states (Prince et al., 2019;
Meyer et al., 2019) (e.g., fractional rms variability, Bayesian blocks, etc).
Here we follow the work of Kushwaha et al. (2016) who showed that the
$\gamma$-ray flux distribution of PKS 1510$-$089 showed two distinct log-
normal profiles, one for the high flux level and one for the low one. The
intersection of these two log-normal fits, at $\rm 10^{-5.853}=1.401\times
10^{-6}\ ph\ cm^{-2}\ s^{-1}$, is taken as the separation between flux levels
and is shown as a green horizontal dotted line running across Panel(a) of Fig.
1. A quick look also reveals that flux changes in different bands are
apparently concurrent (except for the radio) and a more careful inspection
also shows that the degree of change differs between the EM bands. We must
note that data in different EM bands are extracted using different time-bin
durations and thus, are not strictly simultaneous. For example, at
$\gamma$-ray energies, a 10-day duration is used to ensure sufficient counts
for a data point, while such good measurements take just a few minutes at
optical and NIR bands.

####  2.4.1 Segment 1 (MJD 54493 – MJD 55115)

The first segment includes four flaring states. Visual inspection of the
multi-wavelength light curves shows a close correlation between the
$\gamma$-ray, optical, NIR, and radio emissions. The optical and NIR bands
also brighten during the high activity of the last three $\gamma$-ray high
states. The observed increase in radio flux seems to start quasi-
simultaneously with the third $\gamma$-ray high state and reached its peak
during the period of the fourth $\gamma$-ray high state.

####  2.4.2 Segment 2 (MJD 55122 – MJD 56202)

The source showed intense activity in this period, with four major
$\gamma$-ray flaring states. The normalized variability amplitudes are not
consistent between the bands. In addition, there are seasonal gaps in optical
and NIR light curves during the third of the four $\gamma$-ray flaring states
due to the target’s proximity to the Sun. Therefore, it is difficult to tell
by inspection whether there are any correlations between different bands.

####  2.4.3 Segment 3 (MJD 56202 – MJD 56992)

In this segment, the fluctuating trends in the optical and NIR light curves
visually coincide with the single $\gamma$-ray flaring state but the radio
flux density is not as clearly correlated.

####  2.4.4 Segment 4 (MJD 56992 –MJD 58104)

During the first high state of segment 4, four rapid sub-flares were already
identified as A (MJD 57100 to MJD 57150), B (MJD 57150 to MJD 57180), C (MJD
57208 to MJD 57235), and D (MJD 57235 to MJD 57260) in Prince et al. (2019).
The four flares described in their paper appear as four peaks in our
$\gamma$-ray light curve constructed using the 10-day bins. During the high
state period, the optical and NIR emissions also exhibit rapid and complex
variations. Visual inspection shows that the first three $\gamma$-ray
radiation peaks appear almost quasi-simultaneously in the optical band, while
the first two peaks also appear almost quasi-simultaneously in the NIR band.
Meanwhile, it can be clearly seen that the strong $\gamma$-ray emission occurs
during the rising phase of the radio outburst. During the second high state of
segment 4, there are three near-simultaneous low-frequency counterparts to the
$\gamma$-ray flaring state. The prominent radio flare shows longer rise and
decay time scales than $\gamma$-ray and optical bands.

Figure 2: Color index (R $-$ J) vs. J magnitude in each segment; the color
bars indicate the progression of time in MJD. Figure 3: Color index (R $-$ J)
vs. J magnitude, grouped by into low flux states (R magnitude $\textgreater$
16, blue circles) and high flux states (R magnitude $\textless$ 16, red
squares).

###  2.5 Spectral Variations

We calculate the (R$-$J) color index of the datasets to examine color
variability. Some studies indicate that the observed optical/NIR color changes
depending on the current combination of the relativistic jet’s synchrotron
emission and the accretion disk’s thermal emission (Isler et al., 2017; Sarkar
et al., 2019). Compared to the complex color trends found in studies of some
other blazars, PKS 1510$-$089 shows a simple color behavior. It is clear that
the color trends of all segments exhibit a redder-when-brighter (RWB) behavior
(Fig. 2), i.e., the color index increases with increasing brightness.

We arbitrarily define a low state when the R-band magnitude is greater than
16, and vice versa as the high state. During the optical low state observation
periods, the R-band and J-band light curves do not vary rapidly or
drastically. When the optical emission is in its low state, the accretion disk
emission can dominate and the accretion disk luminosity fluctuates over a wide
range of time scales (Lira et al., 2011). Therefore, we investigate the
correlation between the (R$-$J) color index and the J magnitude from another
perspective using the above state division, which is displayed in Fig. 3.

Figure 4: ZDCFs between different wavebands over the entire time of the
observations; a ZDCF peak at positive time lag means the LC1 precedes the LC2
where panels are labeled as ‘LC1’ vs ‘LC2’. The color contours denote the
distribution of random cross-correlations obtained by Monte Carlo simulations;
from dark to light these are 1$\sigma$, 2$\sigma$, and 3$\sigma$.

##  3 VARIABILITY ANALYSIS AND RESULTS

###  3.1 Method

Figure 5: ZDCFs between different wavebands in the four temporal segments; a
ZDCF peak at positive time lag means the LC1 precedes the LC2 where panels are
labeled as ’LC1’ vs ’LC2’. The color contours denote the distribution of
random cross-correlations obtained by Monte Carlo simulations; here the dark
and light are 2$\sigma$ and 3$\sigma$, respectively.

Traditional timing analysis requires that a time series is sampled evenly in
the time domain. In astronomical observations, due to the observation
schedule, weather and technical issues, it is more common that the sampling of
the light curve of each campaign is uneven. We apply the z-transformed
discrete correlation function (ZDCF) method (Alexander, 1997, 2013) to search
for the correlation and time lag between the light curves of each pair of
bands. This method uses equal population binning and Fisher’s z-transform
(Kendall & Stuart, 1969, 1973, and references therein) to correct several
biases of the discrete correlation function (DCF) method (Edelson & Krolik,
1988). This method also ensures that the statistical significance of each bin
is high enough by changing the bin width, a form of adaptive binning. The
specific binning algorithm we used is described in detail in Kramarenko et al.
(2022). After binning, the Pearson’s correlation coefficient ($r$) is
calculated and transformed to the z-space:

$z=\frac{1}{2}\log\left(\frac{1+r}{1-r}\right),~{}\zeta=\frac{1}{2}\log\left(\frac{1+\rho}{1-\rho}\right),~{}r=\tanh
z,$ (1)

where $\rho$ is the unknown population correlation coefficient of the bin. The
ZDCF method uses the Ansatz that $\rho=r$ for the transformation.

The mean of $z$ ($\bar{z}$) and the variance of $z$ ($s_{z}$) are estimated by

$\bar{z}=\zeta+\frac{\rho}{2(n-1)}\times\left[1+\frac{5+\rho^{2}}{4(n-1)}+\frac{11+2\rho^{2}+3\rho^{4}}{8(n-1)^{2}}+\cdots\right],$
(2)

and

$s_{z}^{2}=\frac{1}{n-1}\left[1+\frac{4-\rho^{2}}{2(n-1)}+\frac{22-6\rho^{2}-3\rho^{4}}{6(n-1)^{2}}+\cdots\right].$
(3)

The error on the ZDCF values is determined from 300 Monte Carlo simulations
where light curves are generated by randomly adding error to each data point
according to its observed error, followed by finding their ZDCF estimates. The
mean and variance are then estimated in the z-space, which is then transformed
back to the real space to estimate the error (see Eqn. (8) of Alexander,
2013).

To estimate the significance of the correlation coefficients, we follow the
Monte Carlo method described by Max-Moerbeck et al. (2014a) based on the
simulated light curves having variability and statistical properties similar
to the observed light curve. In this work, we simulated a total of 1000 light
curves with the same power spectral density (PSD) and probability density
function (PDF) as the observed light curves in each band using the algorithm
presented by Emmanoulopoulos et al. (2013), as realized by the
DELightcurveSimulation code (Connolly,
2015)444https://github.com/samconnolly/DELightcurveSimulation. We estimate the
underlying PSD by fitting a smoothly bending power-law model plus a constant,
c:

$\mathcal{P}(f;\bm{\gamma},c)=\frac{Af^{-\alpha_{\text{low
}}}}{1+\left(f/f_{\text{bend }}\right)^{\alpha_{\text{high
}}-\alpha_{\text{low }}}}+c,$ (4)

where $\bm{\gamma}=\left(A,f_{\text{bend }},\alpha_{\text{low
}},\alpha_{\text{high }}\right)$ represents the model parameters, which are
normalization, bend frequency, low-frequency slope, and high-frequency slope,
respectively; $c$ is an additional constant Poisson noise. The optimal
parameters of this PSD model were obtained by maximum likelihood using the
Basin-Hopping algorithm and the Nelder-Mead minimization algorithm provided by
the Python package SciPy. As for the PDF model, we use a mixed distribution
model consisting of a gamma distribution and a log-normal distribution:

$f_{\text{mix
}}(x)=w_{\Gamma}\frac{\theta^{-\kappa}\mathrm{e}^{-x/\theta}x^{\kappa-1}}{\Gamma(\kappa)}+w_{\ln\mathcal{N}}\frac{\mathrm{e}^{-(\ln
x-\mu)^{2}/\left(2\sigma^{2}\right)}}{\sqrt{2\pi}x\sigma},$ (5)

where $\kappa$ and $\theta$ are the shape and the scale parameters of the
gamma distribution; $\mu$ and $\sigma^{2}$ represent the mean and the variance
of the log-normal distribution; $w_{\Gamma}$ and $w_{\ln\mathcal{N}}$ are the
weights of the two distributions, respectively, where their sum is 1. The
optimal PDF model parameters are obtained through a maximum likelihood
analysis using the BFGS (Broyden-Fletcher-Goldfarb-Shannon) minimization
algorithm in the Python package SciPy. Table 3 lists the best-fit parameters
for both models across all bands.

Table 3: The best-fit parameters for both models across all bands | PSD | PDF
---|---|---
| A(Hz-1) |  $f_{\text{bend }}$(Hz) | $\alpha_{\text{low }}$ | $\alpha_{\text{high }}$ | c(Hz-1) | $\kappa$ | $\theta$ | $\mu$ | $\sigma$ | $w_{\Gamma}$ | $w_{\ln\mathcal{N}}$
$\gamma$-ray | 22.8056 | 0.0070 | 0.1202 | 1.7735 | 0.1288 | 5.9916 | 6.01437 | 0.7580 | 4.5668 | 0.0994 | 0.9006
R | 0.0814 | 0.1789 | 0.8749 | 5.8442 | 0.0542 | 3.8340 | 0.7921 | 0.2269 | 1.3622 | 0.4850 | 0.515
J | 0.1263 | 0.1048 | 0.9620 | 2.6107 | 0.0293 | 3.0122 | 1.5823 | 0.2814 | 2.1324 | 0.6625 | 0.3375
37 GHz | 3.4348 | 0.0015 | 0.3314 | 1.8827 | 0.0070 | 5.3792 | 0.6663 | 0.1914 | 2.741 | 1.1674 | $-0.1674$

We calculated the correlation coefficients between the simulated light curves
in the exact same way as done for the observed data. After that, significance
levels for observed data correlation coefficients are derived from the
distribution of simulated correlation coefficients for each time lag. It
should be noted that significance estimated this way assumes completely random
variations in the different bands which is not true; e.g., NIR J and optical R
bands both are dominated by the same synchrotron component. Similarly,
$\gamma$-ray fluxes are related to the optical ones if the same particle
distribution is responsible for both the synchrotron and inverse Compton
components. The synchrotron radio emission probably has contributions from a
much bigger region and also encodes opacity effects. Thus, the quoted
significances are overestimates in general.

###  3.2 Results

The ZDCF results for the entire span of these observations and for the four
individual segments are shown in Figs. 4 and 5, respectively. The time lags at
which the ZDCFs peak are reported in Table 4. Any ZDCF peak within the
sampling time resolution is not considered to indicate an actual time lag
between the two bands’ emissions. We only consider the highest ZDCF peak
closest to 0 lag and discard any peaks near the edges of the temporal span
considered (e.g., $\pm 300$ d for segment1).

####  3.2.1 Interband correlations over the entire time

We first consider interband correlations over the entire span of these
observations. We find correlations close to zero lag with significance larger
than 2$\sigma$ between all frequency bands except for $\gamma$-ray and 37 GHz,
with four of them exceeding 3$\sigma$: $\gamma$-ray and R-band, $\gamma$-ray
and J-band, R-band and J-band, and R-band and 37 GHz (Fig. 4). The result that
these ZDCF peaks near zero lag rarely exceed 0.5 can be understood as the
cumulative effect of incorporating the gaps in the R and J band data between
each segment. We know that even for the same blazar, there can be significant
differences in multi-band correlations during different observational
intervals. Therefore, we study the correlation of light curves between
different bands over the four segments described above, thereby hoping to
minimize the impact of missing data on the intrinsic correlations.

Table 4: Time lags, in days, of the ZDCF peaks between different bands Light curves | ENTIRE TIME | SEGMENT 1 | SEGMENT 2 | SEGMENT 3 | SEGMENT 4
---|---|---|---|---|---
$\gamma$-ray versus R | $13.27_{-4.40}^{+22.60}$ | $-0.81_{-6.85}^{+2.57}$ | $-4.39_{-13.72}^{+2.28}$ | $-4.58_{-6.52}^{+2.96}$ | $-3.82_{-2.52}^{+4.47}$
$\gamma$-ray versus J | $2.81_{-3.17}^{+16.52}$ | $6.57_{-4.23}^{+2.76}$ | $9.85_{-6.94}^{+2.24}$ | $-5.42_{-5.01}^{+3.79}$ | $8.13_{-2.52}^{+4.47}$
$\gamma$-ray versus 37 GHz | $--$ | $51.65_{-3.65}^{+5.24}$ | $39.03_{-5.79}^{+3.77}$ | $--$ | $--$
R versus J | $-25.97_{-8.08}^{+55.87}$ | $0.81_{-11.81}^{+2.17}$ | $0.16_{-0.16}^{+11.78}$ | $2.05_{-14.01}^{+3.75}$ | $-0.42_{-1.53}^{+10.52}$
R versus 37 GHz | $53.5_{-28.02}^{+5.96}$ | $34.01_{-3.48}^{+15.63}$ | $10.89_{-18.07}^{+3.00}$ | $42.20_{-2.58}^{+8.25}$ | $--$
J versus 37 GHz | $39.64_{-84.76}^{+16.01}$ | $41.25_{-2.43}^{+7.56}$ | $19.18_{-32.35}^{+6.55}$ | $38.84_{-3.12}^{+4.91}$ | $--$

####  3.2.2 Interband correlations for each segment

We now consider the results for each segment in examining the correlations
between the bands. In most cases, the significance of the highest peak ZDCF
peak closest to zero lag is greater than 2$\sigma$ and has a value exceeding
0.6.

_$\gamma$ -ray versus optical/NIR:_
There are significant correlations between the $\gamma$-ray photon flux and
the optical R-band and NIR J-band flux densities for all segments. Recall that
we use a time bin of 10 d to reconstruct the $\gamma$-ray light curve. The
results given in the first two rows of Fig. 5 and Table 4 show that all these
ZDCF peaks are within that 10 day interval and hence we measure essentially no
time lag between the $\gamma$-ray and optical/NIR emissions in every segment.

_$\gamma$ -ray versus radio:_
The $\gamma$-ray versus radio ZDCF peaks at non-zero time lags in the first
two segments, with the $\gamma$-ray emission leading the radio emission by
about 50 days in segment 1 and about 40 days in segment 2. In segments 3 and
4, no statistically robust correlations between the $\gamma$-ray and 37 GHz
bands is seen.

_Optical versus NIR:_
For all segments, the optical versus NIR ZDCF peaks are very close zero time
lags, indicating that the R-band and J-band emissions are simultaneous or at
least that any delay between these two light curves occurs on the time-scale
smaller than the cadence of observations.

_Optical/NIR versus radio:_
In segment 1 and segment 3, both the optical and NIR versus radio ZDCF peaks
correspond to about 40 days, with the optical/NIR emission leading the radio
emission. During the MJD 55800 $\sim$ MJD 55900 portion of segment 2 the radio
emission rises rapidly, and optical and NIR data for this period happen to be
missing because of the source’s proximity to the Sun. Therefore it is not
surprising that the ZDCF peak for segment 2 shows only a weak correlation
coefficient between these bands and the errors of the time lags corresponding
to those peaks are very large. In segment 4, as mentioned above, the radio
flares show longer rise and decay timescales than do the changes in the
$\gamma$-ray and optical/NIR bands. Therefore for this fourth temporal
segment, the optical and NIR correlations and broad nominal positive lags we
obtain to the radio emission by applying mathematical methods such as ZDCF are
unconvincing.

##  4 Quasi-periodic Oscillation Analysis

Over the last 15 years or so a modest number of QPO detections have been
reported in blazars in different EM bands. (e.g., Gupta, 2014, 2017; Sarkar et
al., 2019; Roy et al., 2022a, b, and references therein). PKS 1510$-$089 is
among a few blazars which has had published claims of detection of QPOs in
different EM bands in different temporal spans. In 22 and 37 GHz radio bands
data taken during 1995 – 2005, gave indications of QPOs with two periods, P1 =
0.92 yr and P2 = 1.84 yrs, where P1 is a harmonic of P2 (Xie et al., 2008). In
UMRAO (University of Michigan Radio Astronomical Observatory) data taken in 8
and 14.5 GHz during 1974 – 2011, indications of multiple QPO periods in the
period range 430 to 1080 days were found (Li et al., 2021). In a poorly
sampled optical data taken during 1999 – 2001, QPO with period of $\sim$ 0.92
yr was suggested (Xie et al., 2002), though the data train was not long enough
to make this a convincing case. Other poorly sampled NIR K band and optical R
band light curves of the source taken during the longer span of 2006 – 2014
apparently displayed multiple periods ranging from 203 to 490 days
(Sandrinelli et al., 2016). In the 2006 – 2014 Fermi-LAT light curve of this
blazar, a QPO with a period of $\sim$ 115 days was claimed (Sandrinelli et
al., 2016; Castignani et al., 2017). Recently, transient $\gamma-$ray QPOs
with periods of 3.6 days and 92 days in Fermi-LAT data of the source were
claimed in two different portions of the light curve (Roy et al., 2022a).

In light of these previous QPO signals from the blazar PKS 1510$-$089 we now
consider for further analysis the independent new data we have presented here:
the Metsähovi 37 GHz radio light curve between 2008 – 2018 (plotted in the
bottom panel of Fig. 1). Modulations indicating a possible quasi-periodic
component might be identified by visually inspecting the flux variations with
time. We used several methods, the Generalized Lomb-Scargle periodogram
(GLSP), REDFIT, and the Weighted Wavelet Z-transform (WWZ) to assess the
possibility of quasi-periodicity in the frequency domain (and also, for WWZ,
the time domain). All of these methods can be used for unevenly sampled data
which is normally the case for radio light curves observed over the range of
few years to decades.

Figure 6: GLSP of the radio light curve of PKS 1510$-$089 for the duration
2008–2018. The black histogram denotes the normalized GLSP. The solid red,
blue and magenta curves denote 3$\sigma$, 4$\sigma$ and 5$\sigma$ confidence
intervals, respectively. Figure 7: Wavelet result for the duration 2008 – 2018
light curve. Left panel: WWZ plot with red denoting the maximum wavelet power
which decreases towards violet and black. Right panel: Time averaged WWZ
(solid black curve); the blue, green and red dotted curves represent
3$\sigma$, 4$\sigma$ and 5$\sigma$ confidence intervals, respectively. Figure
8: REDFIT result for the 2008 – 2018 light curve. The black solid curve and
the red dashed curve represents the spectrum of input data and bias corrected
spectrum respectively. The blue dashed-dotted curve denotes the AR(1)
theoretical spectrum and the solid magenta curve is the 95$\%$ confidence
curve.

In this work, we employed the GLSP routine of the
PYASTRONOMY555https://github.com/sczesla/PyAstronomy package which is based on
Zechmeister & Kürster (2009). Fig. 6 shows the GLSP plotted against frequency
(black curve). The red, blue and magenta curves respectively represent
3$\sigma$, 4$\sigma$ and 5$\sigma$ significance levels. The peak of the
periodogram is found at the frequency of $1566^{+439}_{-184}$ days at about
4$\sigma$ nominal significance. But we note that with a data train that does
not exceed 4000 d, this possible period amounts to less than 3 cycles and so
can only be considered to hint at a QPO.

The significance of any peak that could be a QPO is estimated by simulating
light curves having properties similar to those of the original light curve.
To quantify the significance of a peak it must be assumed that the light curve
of blazars can be well represented by the stochastic processes occurring in
the accretion disk and the associated jets. For blazars, the PSD over a range
of frequencies also can usually be well fit by a power law with negative
spectral index $\alpha$ (Gupta, 2017; Bhatta & Dhital, 2019, and references
therein). So we use a simple power law as the underlying model for the power
spectrum Any QPO is likely to arise due to coherent disk or jet processes and
is different from the stochastic process and thus would appear as peaks in the
power spectrum (Ackermann et al., 2012), the significance of which can be
assessed. The degree to which any QPO peaks stand above the power law is used
to assess the significance of the periodicities (see for details, Tripathi et
al., 2021). Wavelet functions are commonly used to examine any evolution in
frequency and amplitude of a QPO signal (Torrence & Compo, 1998). As the data
employed in this work are unevenly sampled, we use an improved version of the
wavelet approach, namely the WWZ.

We calculated the WWZ666https://www.aavso.org/software-directory and time-
averaged WWZ for the simulated light curves as well as for the actual data,
following Tripathi et al. (2021) to calculate the significance levels. If the
WWZ amplitude is marginalized over the whole length of the observations, one
gets the WWZ as a function of frequency which is essentially the PSD of the
WWZ, which should follow the power law and will be distributed as $\chi^{2}$
with 2 degrees of freedom. Again we use a simple power law as the underlying
model to fit the time-averaged WWZ. The left panel of Fig. 7 shows the wavelet
density plot. Interestingly, there is only a single peak at $1538^{+24}_{-23}$
days and it is persistent throughout the data stream. The right panel of the
Fig. 7 plots the time-averaged WWZ against the period. The significance of the
highly persistent signal found in the wavelet density plot nominally exceeds
5$\sigma$.

We also employed the REDFIT method777The code is found at
https://www.marum.de/Prof.-Dr.-michael-schulz/Michael-Schulz-Software.html
which essentially compares the data to a first-order auto-regressive process
(AR1) (Schulz & Mudelsee, 2002; Hong et al., 2018; Gupta et al., 2018, and
references therein) and calculates significance based on $\chi^{2}$
distributions. In AR analyzes, the flux at a given point of time is compared
to the data at the past times using a regressive relation and can include any
number of such past values. In the simplest AR1 method, only one value
previous to the data is used to compute the theoretical light curve and
subsequently, the AR1 spectrum. For time $t_{i}$ (i=1,2,…$N$), the AR1 process
$r$ can be written as

$\displaystyle r(t_{i})=\rho_{i}r(t_{i-1})+\epsilon(t_{i})$ (6)
$\displaystyle\rho_{i}=exp(-(t_{i}-t_{i-1})/\tau),$ (7)

where $\epsilon$ is the Gaussian noise such that the AR1 process has unit
variance and zero mean, $\rho$ denotes the autocorrelation coefficient, and
$\tau$ represents the characteristic time scale for the AR1 process. The power
spectrum $S_{rr}(f_{i})$, corresponding to the AR1 process described in Eq. 4
for frequency $f_{i}$ up to the Nyquist frequency, $f_{nyq}$, is given as

$S_{rr}(f_{i})=S_{0}\frac{1-\rho^{2}}{1-2\rho~{}cos(\pi
f_{i}/f_{Nyq})-\rho^{2}}$ (8)

The estimated values of $\rho$ and $\tau$ for this light curve are found to be
0.92 and 206 days, respectively. Fig. 8 plots the normalized REDFIT results
against frequency. The solid black and red dashed curve respectively represent
the spectrum of the input data and the bias corrected spectrum. The blue
dashed-dotted curve shows the theoretical AR(1) spectrum. The magenta curve
represents a 2$\sigma$ confidence curve. The peak around 360 days has
significance exceeding 2$\sigma$ but it is suspiciously close to 1 year. The
otherwise strongest signal is at $1465^{+56}_{-53}$ days and while it is less
than 2$\sigma$ it is consistent with the results from GLSP and WWZ methods.

##  5 DISCUSSION

###  5.1 Unveiling the physics behind the spectral variations

Recently, Otero-Santos et al. (2022) used a powerful statistical tool – non-
negative matrix factorization (NMF, see Paatero & Tapper, 1994; Ivezić et al.,
2014) to reproduce the spectral variability for some blazars over ten years
using two to four components. The photometric and spectropolarimetric data
used in their study (taken from Steward Observatory) are partly homologous to
the optical R-band data in this paper, with temporal coverages slightly (by
less than a year) longer than ours. This study is critical to understanding
the origin of blazar variability. In PKS 1510$-$089, the spectral variability
is explained as the sum of a bright BLR component, a power-law (PL) component
accounting for the non-thermal synchrotron radiation of the relativistic jet,
and the blue accretion disk contribution that shows slight variation and lower
brightness than the BLR and the PL for the jet. Otero-Santos et al. (2022)
concluded that the contribution of the jet component is related to color and
showed an RWB trend consistent with our results. For PKS 1510$-$089, which is
an LSP FSRQ (Ajello et al., 2020) and which has an accretion disk temperature
estimated as 40,000 K (Abdo et al., 2010), the synchrotron peak and the
thermal accretion disk spectrum peak are located at the NIR and the far UV,
respectively. We note that the thermal disk spectrum is bluer, especially in
this source, where the accretion disk temperature is so high. However, non-
thermal photons have a redder spectral shape. So observed colors are
influenced by changes in accretion disk emission, jet emission, or both. It is
reasonable to expect that we would observe a bluer-when-brighter (BWB) trend
in the low flux regime if thermal radiation from optical/NIR wavelengths is
much stronger than or even swamps the jet radiation. However, our
investigation finds the opposite RWB trend. We can infer that the Doppler-
boosted jet emission dominates the entire radiation process, even the optical
band. Moreover, in the optical/NIR high flux state, RWB behavior is also
displayed throughout the observation period. So we also can infer that the
high flux states of optical/NIR emission come from enhanced jet activity, and
increased non-thermal emission results in a reddening of the color.

###  5.2 Multi-band cross-correlation analysis

Blazars are the most powerful, persistent extra-galactic broadband sources.
The region responsible for the broadband emission is highly compact – beyond
the resolution limit of modern facilities even for the nearest sources, e.g.,
M87 and Sgr A* (Event Horizon Telescope Collaboration et al., 2019a, b; Event
Horizon Telescope Collaboration et al., 2022). However, the jet region is
optically thin at NIR and higher energies and thus, temporal flux variability
allows us to access the spatial extension buried under the resolution limit
and is currently the only way to infer these spatial scales. Further, though
rapid and strong flux variability is one of the defining characteristics of
blazars, these high states are rarely accompanied by the strong spectral
changes that would indicate a departure from its SED class.

####  5.2.1 $\gamma$-ray versus optical/NIR

The notable correlation between $\gamma$-ray and optical/NIR emissions
indicate simultaneous variations (within the sampling intervals) and thus
strongly supports the leptonic model. Ramakrishnan et al. (2016) studied 15
blazars and reported two sources with a strong correlation in optical and
$\gamma$-ray emissions with zero time delay, one of which is PKS 1510$-$089\.
We validated these previous results using high-quality data over a more
extended period. In the leptonic scenario, the high-energy emission is
produced by relativistic electrons scattering seed photons from low-energy to
high-energy through the inverse Compton (IC) process. Studies employing SED
modeling of PKS 1510$-$089 shows that synchrotron radiation from relativistic
electrons can adequately describe the mm to optical emission, while high-
energy emissions (X-ray to $\gamma$-ray) are not. These previous studies have
shown that X-ray and $\gamma$-ray emissions are also an EC process of the BLR
and DT photons (Kataoka et al., 2008b; Abdo et al., 2010; Brown, 2013;
Barnacka et al., 2014; Prince et al., 2019). Moreover, Prince et al. (2019)
inferred that the $\gamma$-ray and optical emissions are produced in the same
region from the flux doubling times of the $\gamma$-ray and optical and UV
light curves.

Figure 9: Archived 15 GHz VLBI data light curves for PKS 1510$-$089\. Red
circles represent the fitted core features flux density from Angelakis et al.
(2019), and blue circles represent the VLBI total flux density from the MOJAVE
web page https://www.cv.nrao.edu/MOJAVE/sourcepages/1510-089.shtml (Lister et
al., 2018). The values that make up the upper border of the green shade
represent the ratio of the core flux density to the total flux density for
each epoch.

####  5.2.2 Optical versus NIR

As stated above, the synchrotron radiation of the relativistic jet is the
primary source of its optical and NIR radiation (e.g. Blandford & Rees, 1978).
However, additional contributions are expected from thermal radiation from the
accretion disk and torus. A statistical study suggests that the host galaxy
and BLR may also contribute significantly to the optical band (Otero-Santos et
al., 2022). Quantifying all possible optical and NIR emission contributions
has been a challenge. The combination of the color trend and correlation
analysis of optical and NIR light curves is an effective way to understand the
radiation mechanism of this source. We determined that the jet emission
dominates the optical and NIR color variability behavior throughout the
observation period.

The correlation between the optical B-band and NIR J-band observation data of
the SMARTS monitoring program from 2008 to 2010 was studied at an early stage
(Bonning et al., 2012). We use the R-band data released by this program and
supplementary data from the Steward Observatory to revisit the correlation
between the optical and NIR emissions with longer observation times and denser
sampling intervals. We confirm the results of Bonning et al. (2012) that there
is a very good correlation between the light curves of optical and NIR bands.
Significant correlations with zero time lag are seen between the optical
R-band and NIR J-band light curves, indicating that the bulk of these
emissions are simultaneous or that any time lag occurs on a timescale smaller
than the observation cadence. It can further be inferred that the main
emission region of the optical and NIR photons is co-spatial and originates in
the jet.

####  5.2.3 $\gamma$-ray versus radio

Several studies have reported the time delay between the $\gamma$-ray emission
and the radio emission measured by both a single radio dish and VLBI. For
single dish data, $\gamma$-ray flares occur earlier than radio flares (Max-
Moerbeck et al., 2014b; Fuhrmann et al., 2014). And for VLBI data, Pushkarev
et al. (2010) used a large sample of 183 bright Fermi-detected sources and
found that the correlations were strong between the $\gamma$-ray emission and
the VLBA core component emission with $\gamma$-ray emission preceding radio
emission. Recently, they again performed the correlation analysis for the
larger samples that have been accumulated and found the time delay is in good
agreement with their previous results (Kramarenko et al., 2022). The light
curve obtained from single-dish observations or VLBA observations yielded
similar results when correlated with the $\gamma$-ray light curve. The above
results are supported for two reasons. Firstly, the radio core dominates the
radio emission, which makes the core flux density and total flux density light
curves’ characteristics similar in VLBA observations. Secondly, the parsec-
scale emission region dominates the radio flux density variation, which
reflects that the characteristics of total light curves obtained by the
single-dish observation and VLBA observation are similar. Similar
characteristics refer to the flares occurring almost quasi-simultaneously,
with roughly the same upward and downward trends, although the above
characteristics vary from source to source.

We can compare the magnitudes of core and total flux densities of PKS
1510$-$089 at multiple epochs when 15 GHz VLBI observations were made by
collecting data released by the MOJAVE program (see Fig. 9). Here we also
calculate the values ($S_{\mathrm{core}}/S_{\mathrm{VLBA}}$) that characterize
the core dominance degree of each epoch (also see Fig. 9), where
$S_{\mathrm{core}}$ and $S_{\mathrm{VLBA}}$ are the core and total flux
density, respectively. It can be clearly seen that this source is core
dominated. In addition, for PKS 1510$-$089, it has been demonstrated that the
parsec-scale emission region dominates the radio flux density variability
(Orienti et al., 2013b). So, the above $S_{\mathrm{core}}$ and
$S_{\mathrm{VLBA}}$ variations are not much different from the radio emission
measured by a single antenna. We therefore can reasonably predict that if
there were even more densely sampled VLBI data available we would obtain
$\gamma$-ray and radio radiation correlations that are consistent with single-
dish data analyses. We should remember that intensive monitoring is essential
to understand the correlation between the cross-band activity for the long-
term trends, so our choice of the well-sampled 37 GHz single-dish data to
perform the correlation analysis is sensible.

In the classical and commonly used shock-in-jet model (e.g., Marscher & Gear,
1985; Türler et al., 2000; Fromm et al., 2011), the shock propagates down a
conical jet and accelerates the relativistic particles at the shock front.
These particles propagate behind the shock front and lose their energy through
different energy loss mechanisms, such as adiabatic expansion and synchrotron
radiation. The most widely-accepted view is that the moving shock is an
intrinsic factor responsible for the observed radio flux density variation.

As mentioned above, the current investigations found that $\gamma$-ray
emission precedes the radio emission. Several studies have attempted to
interpret this result, suggesting that $\gamma$-ray luminosity variations may
be connected to the same shocked radio features and connect these two
frequency emissions through the shock (Fuhrmann et al., 2014; Ramakrishnan et
al., 2015). Suppose radio and $\gamma$-ray emissions are triggered by shocks
propagating in the relativistic jet. In this scenario, the time lag is related
to the distance between the radio core and the location where the $\gamma$-ray
emission is produced. This can also be understood as the unresolved core
centroid moving towards the position of the black hole as observing frequency
increases (Blandford & Königl, 1979; Konigl, 1981), i.e., the so-called ‘core
shift’ effect. The $\gamma$-ray photons escape from the jet immediately when
the shock arrives, and it takes a while for the shock (perturbation) to
propagate farther along the jet until it reaches the $\tau=1$ surface from
which radio photons at the specified radio frequency can escape.

In our study, we also found that the $\gamma$-ray emission leading the radio
emission in segment 1 and segment 2. Moreover, as shown above, the flux
density variability of PKS 1510$-$089 is dominated by the flux density of the
core fitted in the VLBI images. Hence, a possible explanation for these
results may be the above scenario involving a shock. We estimate the
separation of the $\gamma$-ray and the 37 GHz emission region following the
method and assumption of Pushkarev et al. (2010),

$\Delta r=r_{\gamma}-r_{\mathrm{Radio}}=\frac{\delta\Gamma\beta c\Delta
t_{\mathrm{Radio}-\gamma}^{\mathrm{obs}}}{1+z}=\frac{\beta_{\mathrm{app}}c\Delta
t_{\mathrm{Radio}-\gamma}^{\mathrm{obs}}}{(1+z)\sin\theta}\\\ $ (9)

where $\beta_{\mathrm{app}}$ is the apparent speed, $z$ is the redshift,
$\theta$ is the viewing angle and the $\Delta
t_{\mathrm{Radio}-\gamma}^{\mathrm{obs}}$ is the time lag in the observer’s
frame. Using $\beta_{\mathrm{app}}$ $=$ 28 (Lister et al., 2021), $\theta$ $=$
2.5∘ (Homan et al., 2021) and $\Delta
t_{\mathrm{Radio}-\gamma}^{\mathrm{obs}}$ $=$ 45 days which is the average of
the time lags of the two segments with significant correlation, the distance
between the two emission regions is estimated as 17.82 pc, and the
corresponding projected distance is 0.78 pc.

From the core shift measure ($\Omega$) (defined in Lobanov, 1998) for PKS
1510$-$089 by Pushkarev et al. (2012), we can estimate the distance between
the radio core and the true jet apex as 8.36 pc. So, the value of the distance
between the $\gamma$-ray emission region and the jet apex is a negative value.
Above results are in agreement with Kramarenko et al. (2022) findings which
calculated the radio lagging $\gamma$-ray emission $111_{-30}^{+16}$ days in
the observer’s frame and the $\gamma$-ray emission generation region is
$-12.00_{-8.65}^{+13.71}$ pc from the central engine. The jet is widely
thought to be launched in the vicinity of the supermassive black hole within
100 $R_{s}$ (Schwarzschild radius) (Meier et al., 2001). Referring to previous
studies, the black hole mass of PKS 1510$-$089 is estimated to be on the order
of $10^{7}$ to $10^{8}$ solar masses (Oshlack et al., 2002; Xie et al., 2005;
Liu et al., 2006; Abdo et al., 2010; Park & Trippe, 2017; Castignani et al.,
2017; Rakshit, 2020). Therefore, the jet apex is $10^{-5}$ to $10^{-6}$ pc
away from the central engine or supermassive black hole (SMBH), which is
negligible. The negative value of the distance from the location of the
$\gamma$-ray emission production to the central engine might further suggest
that the conjecture that the perturbation propagates in a straight jet
geometry is too simple; Kramarenko et al. (2022) offer several alternative
interpretations.

####  5.2.4 Optical/NIR versus radio

For this source, as discussed above, the optical and NIR emissions are
dominated by non-thermal synchrotron emission from the relativistic jet and
their emission region is almost co-spatial. We see that the optical and radio
ZDCF profiles are very similar to the NIR and radio ZDCF profiles, while the
lack of some J-band data near MJD 57500 compared to the R-band data results in
some visible difference between their profiles in segment 4. Meanwhile, due to
the seasonal gap in the data of this particular segment, we cannot arrive at a
firm conclusion about the correlation results in this segment. However, the
remaining three segments show that radio is lagging the optical/NIR radiation.
These results are in agreement with the above model; in other words, the
scenario in which radio emission lags behind optical/NIR emission is analogous
to the scenario in which $\gamma$-ray emission lagging radio emission.
Although no firm conclusions have been drawn to date on the correlations
between optical/NIR and radio emissions for blazars in general, correlation
analyses beween optical/NIR and radio bands in large samples or individual
sources have revealed radio emission lagging the optical/NIR in a substantial
fraction of cases (e.g. Zhang et al., 2017; Sarkar et al., 2019).

###  5.3 QPO analysis

We found a hint of a possible QPO of about 1580 days in the 37 GHz radio
emission, but it might only be confirmed through even longer future studies.
Since the radio is synchrotron emission from the jet, jet-precession seems a
potential interpretation of such a lengthy putative QPO. However, almost the
entire multi-band emission is from the jet and thus, similar temporal profiles
are expected across the entire EM bands unless there is a strong spectral
change, which does not seem to be the case here. Though the radio emission
correlates very well with emission in other bands, the temporal profile is
very different, with much smoother and broader peaks. Also, a large fraction
of the radio photons often come from the unresolved core, as in the present
case (see Fig. 9).

##  6 Summary

In this study, we performed a correlation analysis of the multi-band light
curves of the FSRQ PKS 1510$-$089, for $\gamma$-ray, optical, NIR, and 37 GHz
radio bands data sets spanning around ten years (February 2008 to September
2017). The $\gamma$-ray data come from the Large Area Telescope onboard the
space-based Fermi Gamma-ray Space Telescope, the optical R-band data and NIR
J-band data are from the public archive of the SMARTS and the Steward
Observatory telescopes, and the radio observations at 37 GHz were made at the
14 m radio telescope at the Metsähovi Radio Observatory of Aalto University.
We also employed several methods to search for the QPO phenomenon in the radio
light curve between 2008 – 2018. Through a series of analyses, we obtain the
following results:

$\bullet$ This source shows RWB behavior whether compared during individual
segments or grouped by brightness states, indicating that jet emissions
dominate accretion disk emissions throughout the observation period, even in
low flux states.
$\bullet$ All segments show strong $\gamma$-ray and optical/NIR correlations,
and the correlation between $\gamma$-ray and radio is shown in the first two
segments with a non-zero time lag. The existence of these correlations lends
support to the lepton model of $\gamma$-ray emission. For PKS 1510$-$089, the
nominal distance from the location of the $\gamma$-ray emission production to
the central engine is a negative value using a simple model of shock
propagating down a conical jet. This could indicate that the morphology of the
conical jet may have been altered by jet dynamics, making the model unable to
support such a complex situation.
$\bullet$ The correlation analysis reveals a clear, and unsurprising,
correlation between the optical R-band and NIR J-band emissions, implying that
these flares are simultaneous, or at least that the time lag between these
light curves is shorter than the cadence of the observations. Additionally,
the optical and NIR flares are caused mainly by synchrotron radiation of the
jet, and their emission regions are almost co-spatial.
$\bullet$ We confirm that variations in the optical/NIR band lead to
variations in the 37 GHz radio band, and these results again suggest that the
origin of the radiation that dominates these three frequencies is the same.
This correlation indirectly indicates that other components, such as thermal
radiation, make lesser contributions to the total optical and NIR fluxes. This
conclusion of jet dominance is further supported by the color index changes in
the optical and NIR bands.
$\bullet$ A nominal QPO signal is found in the radio light curve, but the
period of $\sim 1580$ d is too long compared to the length of the observations
to be trusted.

## Acknowledgments

This publication makes use of data obtained at Metsähovi Radio Observatory,
operated by Aalto University in Finland. The various diligent observers of
Aalto University in Finland are thankfully acknowledged. This paper has made
use of up-to-date SMARTS optical/near-infrared light curves that are available
at www.astro.yale.edu/smarts/glast/home.php. SMARTS observations of Large Area
Telescope-monitored blazars are supported by Yale University and Fermi GI
grant NNX 12AP15G, and the SMARTS 1.3-m observing queue received support from
NSF grant AST-0707627. Data from the Steward Observatory spectropolarimetric
monitoring project were used. This programme is supported by Fermi GI grants
NNX08AW56G, NNX09AU10G, NNX12AO93G and NNX15AU81G. This work has made use of
publicly available Fermi-LAT data obtained from FSSC’s website data server and
provided by NASA Goddard Space Flight Center. This research has made use of
data from the MOJAVE database that is maintained by the MOJAVE team (Lister et
al., 2018).

This work is supported by the National Key R&D Intergovernmental Cooperation
Program of China (2022YFE0133700), the Regional Collaborative Innovation
Project of Xinjiang Uyghur Autonomous Region (2022E01013), the National
Natural Science Foundation of China (12173078, 11773062), and the Chinese
Academy of Sciences (CAS) ‘Light of West China’ Program (2017-XBQNXZ-A-008,
2021-XBQNXZ-005). PK acknowledges support from the Department of Science and
Technology (DST), government of India, through the DST-INSPIRE Faculty grant
(DST/INSPIRE/04/2020/002586). ACG is partially supported by Chinese Academy of
Sciences (CAS) President’s International Fellowship Initiative (PIFI) (grant
no. 2016VMB073). The work of AT was supported by the Innovation Program of the
Shanghai Municipal Education Commission, Grant No. 2019-01-07-00-07-E00035,
and the National Natural Science Foundation of China (NSFC), grant no.
11973019. MFG acknowledges support from the National Natural Science
Foundation of China (grant no. 11873073), Shanghai Pilot Program for Basic
Research Chinese Academy of Science, Shanghai Branch (JCYJ-SHFY2021-013), and
the science research grants from the China Manned Space Project (No.
CMSCSST-2021-A06).

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, ApJ, 721, 1425, doi: 10.1088/0004-637X/721/2/1425
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Acharyya et al. (2021) Acharyya, A., Chadwick, P. M., & Brown, A. M. 2021, MNRAS, 500, 5297, doi: 10.1093/mnras/staa3483
  * Ackermann et al. (2012) Ackermann, M., Ajello, M., Ballet, J., et al. 2012, ApJ, 751, 159, doi: 10.1088/0004-637X/751/2/159
  * Aharonian et al. (2002) Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2002, A&A, 393, 89, doi: 10.1051/0004-6361:20021005
  * Aharonian et al. (2005a) Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2005a, A&A, 430, 865, doi: 10.1051/0004-6361:20041853
  * Aharonian et al. (2005b) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2005b, A&A, 442, 895, doi: 10.1051/0004-6361:20053353
  * Aharonian et al. (2007) —. 2007, ApJ, 664, L71, doi: 10.1086/520635
  * Aharonian et al. (2009) Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 502, 749, doi: 10.1051/0004-6361/200912128
  * Ahnen et al. (2017) Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2017, A&A, 603, A29, doi: 10.1051/0004-6361/201629960
  * Ajello et al. (2020) Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105, doi: 10.3847/1538-4357/ab791e
  * Aleksić et al. (2014) Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014, A&A, 569, A46, doi: 10.1051/0004-6361/201423484
  * Alexander (1997) Alexander, T. 1997, in Astrophysics and Space Science Library, Vol. 218, Astronomical Time Series, ed. D. Maoz, A. Sternberg, & E. M. Leibowitz, 163, doi: 10.1007/978-94-015-8941-3_14
  * Alexander (2013) Alexander, T. 2013, arXiv e-prints, arXiv:1302.1508.  https://arxiv.org/abs/1302.1508
  * Angelakis et al. (2019) Angelakis, E., Fuhrmann, L., Myserlis, I., et al. 2019, A&A, 626, A60, doi: 10.1051/0004-6361/201834363
  * Arsioli & Chang (2018) Arsioli, B., & Chang, Y. L. 2018, A&A, 616, A63, doi: 10.1051/0004-6361/201833005
  * Barnacka et al. (2014) Barnacka, A., Moderski, R., Behera, B., Brun, P., & Wagner, S. 2014, A&A, 567, A113, doi: 10.1051/0004-6361/201322205
  * Beaklini et al. (2017) Beaklini, P. P. B., Dominici, T. P., & Abraham, Z. 2017, A&A, 606, A87, doi: 10.1051/0004-6361/201731118
  * Bhatta & Dhital (2019) Bhatta, G., & Dhital, N. 2019, doi: 10.3847/1538-4357/ab7455
  * Blandford & Königl (1979) Blandford, R. D., & Königl, A. 1979, ApJ, 232, 34, doi: 10.1086/157262
  * Blandford & Rees (1978) Blandford, R. D., & Rees, M. J. 1978, in BL Lac Objects, ed. A. M. Wolfe, 328–341
  * Błażejowski et al. (2005) Błażejowski, M., Blaylock, G., Bond, I. H., et al. 2005, ApJ, 630, 130, doi: 10.1086/431925
  * Bonning et al. (2012) Bonning, E., Urry, C. M., Bailyn, C., et al. 2012, ApJ, 756, 13, doi: 10.1088/0004-637X/756/1/13
  * Böttcher (2019) Böttcher, M. 2019, Galaxies, 7, 20, doi: 10.3390/galaxies7010020
  * Böttcher et al. (2013) Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, ApJ, 768, 54, doi: 10.1088/0004-637X/768/1/54
  * Brown (2013) Brown, A. M. 2013, MNRAS, 431, 824, doi: 10.1093/mnras/stt218
  * Buxton et al. (2012) Buxton, M. M., Bailyn, C. D., Capelo, H. L., et al. 2012, AJ, 143, 130, doi: 10.1088/0004-6256/143/6/130
  * Castignani et al. (2017) Castignani, G., Pian, E., Belloni, T. M., et al. 2017, A&A, 601, A30, doi: 10.1051/0004-6361/201629775
  * Chen et al. (2012) Chen, X., Fossati, G., Böttcher, M., & Liang, E. 2012, MNRAS, 424, 789, doi: 10.1111/j.1365-2966.2012.21283.x
  * Connolly (2015) Connolly, S. D. 2015, A Python Code for the Emmanoulopoulos et al. [arXiv:1305.0304] Light Curve Simulation Algorithm, arXiv, doi: 10.48550/ARXIV.1503.06676
  * Costamante et al. (2018) Costamante, L., Cutini, S., Tosti, G., Antolini, E., & Tramacere, A. 2018, MNRAS, 477, 4749, doi: 10.1093/mnras/sty887
  * Edelson & Krolik (1988) Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646, doi: 10.1086/166773
  * Emmanoulopoulos et al. (2013) Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013, MNRAS, 433, 907, doi: 10.1093/mnras/stt764
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019a, ApJ, 875, L2, doi: 10.3847/2041-8213/ab0c96
  * Event Horizon Telescope Collaboration et al. (2019b) —. 2019b, ApJ, 875, L4, doi: 10.3847/2041-8213/ab0e85
  * Event Horizon Telescope Collaboration et al. (2022) —. 2022, ApJ, 930, L15, doi: 10.3847/2041-8213/ac667410.3847/2041-8213/ac667210.3847/2041-8213/ac6736
  * Fernandes et al. (2020) Fernandes, S., Patiño-Álvarez, V. M., Chavushyan, V., Schlegel, E. M., & Valdés, J. R. 2020, MNRAS, 497, 2066, doi: 10.1093/mnras/staa2013
  * Fossati et al. (2008) Fossati, G., Buckley, J. H., Bond, I. H., et al. 2008, ApJ, 677, 906, doi: 10.1086/527311
  * Fromm et al. (2011) Fromm, C. M., Perucho, M., Ros, E., et al. 2011, A&A, 531, A95, doi: 10.1051/0004-6361/201116857
  * Fuhrmann et al. (2014) Fuhrmann, L., Larsson, S., Chiang, J., et al. 2014, MNRAS, 441, 1899, doi: 10.1093/mnras/stu540
  * Fuhrmann et al. (2016) Fuhrmann, L., Angelakis, E., Zensus, J. A., et al. 2016, A&A, 596, A45, doi: 10.1051/0004-6361/201528034
  * Gaur et al. (2014) Gaur, H., Gupta, A. C., Wiita, P. J., et al. 2014, ApJ, 781, L4, doi: 10.1088/2041-8205/781/1/L4
  * Gupta (2014) Gupta, A. C. 2014, Journal of Astrophysics and Astronomy, 35, 307, doi: 10.1007/s12036-014-9219-7
  * Gupta (2017) Gupta, A. C. 2017, Galaxies, 6, 1, doi: 10.3390/galaxies6010001
  * Gupta et al. (2018) Gupta, A. C., Tripathi, A., Wiita, P. J., et al. 2018, A&A, 616, L6, doi: 10.1051/0004-6361/201833629
  * Gupta et al. (2017) Gupta, A. C., Mangalam, A., Wiita, P. J., et al. 2017, MNRAS, 472, 788, doi: 10.1093/mnras/stx2072
  * Gupta et al. (2019) Gupta, A. C., Gaur, H., Wiita, P. J., et al. 2019, AJ, 157, 95, doi: 10.3847/1538-3881/aafe7d
  * Gupta et al. (2022) Gupta, A. C., Kushwaha, P., Carrasco, L., et al. 2022, ApJS, 260, 39, doi: 10.3847/1538-4365/ac6c2c
  * H. E. S. S. Collaboration et al. (2013) H. E. S. S. Collaboration, Abramowski, A., Acero, F., et al. 2013, A&A, 554, A107, doi: 10.1051/0004-6361/201321135
  * Hodgson et al. (2017) Hodgson, J. A., Krichbaum, T. P., Marscher, A. P., et al. 2017, A&A, 597, A80, doi: 10.1051/0004-6361/201526727
  * Homan et al. (2021) Homan, D. C., Cohen, M. H., Hovatta, T., et al. 2021, ApJ, 923, 67, doi: 10.3847/1538-4357/ac27af
  * Hong et al. (2018) Hong, S., Xiong, D., & Bai, J. 2018, AJ, 155, 31, doi: 10.3847/1538-3881/aa9d89
  * Horan et al. (2009) Horan, D., Acciari, V. A., Bradbury, S. M., et al. 2009, ApJ, 695, 596, doi: 10.1088/0004-637X/695/1/596
  * Isler et al. (2017) Isler, J. C., Urry, C. M., Coppi, P., et al. 2017, ApJ, 844, 107, doi: 10.3847/1538-4357/aa79fc
  * Ivezić et al. (2014) Ivezić, Ž., Connelly, A. J., VanderPlas, J. T., & Gray, A. 2014, Statistics, Data Mining, and Machine Learning in Astronomy (Princeton University Press; 1st edition (January 12, 2014))
  * 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
  * Kataoka et al. (2008a) Kataoka, J., Madejski, G., Sikora, M., et al. 2008a, ApJ, 672, 787, doi: 10.1086/523093
  * Kataoka et al. (2008b) —. 2008b, ApJ, 672, 787, doi: 10.1086/523093
  * Katarzyński et al. (2001) Katarzyński, K., Sol, H., & Kus, A. 2001, A&A, 367, 809, doi: 10.1051/0004-6361:20000538
  * Kellermann et al. (1998) Kellermann, K. I., Vermeulen, R. C., Zensus, J. A., & Cohen, M. H. 1998, AJ, 115, 1295, doi: 10.1086/300308
  * Kendall & Stuart (1973) Kendall, M., & Stuart, A. 1973, Kendall3The Advanced Theory of Statistics1973
  * Kendall & Stuart (1969) Kendall, M. G., & Stuart, A. 1969, The advanced theory of statistics (Vol. 1, Ch. 11), Griffin & Comp., London
  * Konigl (1981) Konigl, A. 1981, ApJ, 243, 700, doi: 10.1086/158638
  * Kramarenko et al. (2022) Kramarenko, I. G., Pushkarev, A. B., Kovalev, Y. Y., et al. 2022, MNRAS, 510, 469, doi: 10.1093/mnras/stab3358
  * Krawczynski et al. (2001) Krawczynski, H., Sambruna, R., Kohnle, A., et al. 2001, ApJ, 559, 187, doi: 10.1086/322364
  * Kushwaha et al. (2016) Kushwaha, P., Chandra, S., Misra, R., et al. 2016, ApJ, 822, L13, doi: 10.3847/2041-8205/822/1/L13
  * Kushwaha et al. (2017a) Kushwaha, P., Gupta, A. C., Misra, R., & Singh, K. P. 2017a, MNRAS, 464, 2046, doi: 10.1093/mnras/stw2440
  * Kushwaha et al. (2013) Kushwaha, P., Sahayanathan, S., & Singh, K. P. 2013, MNRAS, 433, 2380, doi: 10.1093/mnras/stt904
  * Kushwaha et al. (2017b) Kushwaha, P., Sinha, A., Misra, R., Singh, K. P., & de Gouveia Dal Pino, E. M. 2017b, ApJ, 849, 138, doi: 10.3847/1538-4357/aa8ef5
  * Kushwaha et al. (2018a) Kushwaha, P., Gupta, A. C., Wiita, P. J., et al. 2018a, MNRAS, 473, 1145, doi: 10.1093/mnras/stx2394
  * Kushwaha et al. (2018b) —. 2018b, MNRAS, 479, 1672, doi: 10.1093/mnras/sty1499
  * Li et al. (2021) Li, X.-P., Zhao, L., Yan, Y., et al. 2021, Journal of Astrophysics and Astronomy, 42, 92, doi: 10.1007/s12036-021-09773-9
  * Lira et al. (2011) Lira, P., Arévalo, P., Uttley, P., McHardy, I., & Breedt, E. 2011, MNRAS, 415, 1290, doi: 10.1111/j.1365-2966.2011.18774.x
  * Lister et al. (2018) Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12, doi: 10.3847/1538-4365/aa9c44
  * Lister et al. (2021) Lister, M. L., Homan, D. C., Kellermann, K. I., et al. 2021, ApJ, 923, 30, doi: 10.3847/1538-4357/ac230f
  * Liu et al. (2006) Liu, Y., Jiang, D. R., & Gu, M. F. 2006, ApJ, 637, 669, doi: 10.1086/498639
  * Lobanov (1998) Lobanov, A. P. 1998, A&A, 330, 79.  https://arxiv.org/abs/astro-ph/9712132
  * MAGIC Collaboration et al. (2018) MAGIC Collaboration, Acciari, V. A., Ansoldi, S., et al. 2018, A&A, 619, A159, doi: 10.1051/0004-6361/201833618
  * Mannheim & Biermann (1992) Mannheim, K., & Biermann, P. L. 1992, A&A, 253, L21
  * Marscher & Gear (1985) Marscher, A. P., & Gear, W. K. 1985, ApJ, 298, 114, doi: 10.1086/163592
  * 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
  * Max-Moerbeck et al. (2014a) Max-Moerbeck, W., Richards, J. L., Hovatta, T., et al. 2014a, MNRAS, 445, 437, doi: 10.1093/mnras/stu1707
  * Max-Moerbeck et al. (2014b) Max-Moerbeck, W., Hovatta, T., Richards, J. L., et al. 2014b, MNRAS, 445, 428, doi: 10.1093/mnras/stu1749
  * Meier et al. (2001) Meier, D. L., Koide, S., & Uchida, Y. 2001, Science, 291, 84, doi: 10.1126/science.291.5501.84
  * Meyer et al. (2019) Meyer, M., Scargle, J. D., & Blandford, R. D. 2019, ApJ, 877, 39, doi: 10.3847/1538-4357/ab1651
  * Mücke & Protheroe (2001) Mücke, A., & Protheroe, R. J. 2001, Astroparticle Physics, 15, 121, doi: 10.1016/S0927-6505(00)00141-9
  * Orienti et al. (2013a) Orienti, M., Koyama, S., D’Ammando, F., et al. 2013a, MNRAS, 428, 2418, doi: 10.1093/mnras/sts201
  * Orienti et al. (2013b) —. 2013b, MNRAS, 428, 2418, doi: 10.1093/mnras/sts201
  * Oshlack et al. (2002) Oshlack, A. Y. K. N., Webster, R. L., & Whiting, M. T. 2002, ApJ, 576, 81, doi: 10.1086/341729
  * Otero-Santos et al. (2022) Otero-Santos, J., Acosta-Pulido, J. A., Becerra González, J., et al. 2022, MNRAS, 511, 5611, doi: 10.1093/mnras/stac475
  * Paatero & Tapper (1994) Paatero, P., & Tapper, U. 1994, Environmetrics, 5, 111
  * Park & Trippe (2017) Park, J., & Trippe, S. 2017, ApJ, 834, 157, doi: 10.3847/1538-4357/834/2/157
  * Petry et al. (2000) Petry, D., Böttcher, M., Connaughton, V., et al. 2000, ApJ, 536, 742, doi: 10.1086/308955
  * Prince et al. (2019) Prince, R., Gupta, N., & Nalewajko, K. 2019, ApJ, 883, 137, doi: 10.3847/1538-4357/ab3afa
  * Pushkarev et al. (2012) Pushkarev, A. B., Hovatta, T., Kovalev, Y. Y., et al. 2012, A&A, 545, A113, doi: 10.1051/0004-6361/201219173
  * Pushkarev et al. (2010) Pushkarev, A. B., Kovalev, Y. Y., & Lister, M. L. 2010, ApJ, 722, L7, doi: 10.1088/2041-8205/722/1/L7
  * Rakshit (2020) Rakshit, S. 2020, A&A, 642, A59, doi: 10.1051/0004-6361/202038324
  * Ramakrishnan et al. (2015) Ramakrishnan, V., Hovatta, T., Nieppola, E., et al. 2015, MNRAS, 452, 1280, doi: 10.1093/mnras/stv321
  * Ramakrishnan et al. (2016) Ramakrishnan, V., Hovatta, T., Tornikoski, M., et al. 2016, MNRAS, 456, 171, doi: 10.1093/mnras/stv2653
  * Rebillot et al. (2006) Rebillot, P. F., Badran, H. M., Blaylock, G., et al. 2006, ApJ, 641, 740, doi: 10.1086/500653
  * Richards et al. (2011) Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29, doi: 10.1088/0067-0049/194/2/29
  * Roy et al. (2022a) Roy, A., Sarkar, A., Chatterjee, A., et al. 2022a, MNRAS, 510, 3641, doi: 10.1093/mnras/stab3701
  * Roy et al. (2022b) Roy, A., Chitnis, V. R., Gupta, A. C., et al. 2022b, MNRAS, 513, 5238, doi: 10.1093/mnras/stac1287
  * Sandrinelli et al. (2016) Sandrinelli, A., Covino, S., Dotti, M., & Treves, A. 2016, AJ, 151, 54, doi: 10.3847/0004-6256/151/3/54
  * Sarkar et al. (2019) Sarkar, A., Chitnis, V. R., Gupta, A. C., et al. 2019, ApJ, 887, 185, doi: 10.3847/1538-4357/ab5281
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Schulz & Mudelsee (2002) Schulz, M., & Mudelsee, M. 2002, Computers & Geosciences, 28, 421
  * Smith et al. (2009) Smith, P. S., Montiel, E., Rightley, S., et al. 2009, arXiv e-prints, arXiv:0912.3621.  https://arxiv.org/abs/0912.3621
  * Tanner et al. (1996) Tanner, A. M., Bechtold, J., Walker, C. E., Black, J. H., & Cutri, R. M. 1996, AJ, 112, 62, doi: 10.1086/117988
  * Tavecchio et al. (2000) Tavecchio, F., Maraschi, L., Sambruna, R. M., & Urry, C. M. 2000, ApJ, 544, L23, doi: 10.1086/317292
  * Teraesranta et al. (1998) Teraesranta, H., Tornikoski, M., Mujunen, A., et al. 1998, A&AS, 132, 305, doi: 10.1051/aas:1998297
  * Thompson et al. (1990) Thompson, D. J., Djorgovski, S., & de Carvalho, R. 1990, PASP, 102, 1235, doi: 10.1086/132758
  * Torrence & Compo (1998) Torrence, C., & Compo, G. P. 1998, Bulletin of the American Meteorological society, 79, 61
  * Tripathi et al. (2021) Tripathi, A., Gupta, A. C., Aller, M. F., et al. 2021, MNRAS, 501, 5997, doi: 10.1093/mnras/stab058
  * Türler et al. (2000) Türler, M., Courvoisier, T. J. L., & Paltani, S. 2000, A&A, 361, 850\.  https://arxiv.org/abs/astro-ph/0008480
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * Urry et al. (1993) Urry, C. M., Maraschi, L., Edelson, R., et al. 1993, ApJ, 411, 614, doi: 10.1086/172864
  * Weaver et al. (2022) Weaver, Z. R., Jorstad, S. G., Marscher, A. P., et al. 2022, ApJS, 260, 12, doi: 10.3847/1538-4365/ac589c
  * Xie et al. (2002) Xie, G. Z., Liang, E. W., Zhou, S. B., et al. 2002, MNRAS, 334, 459, doi: 10.1046/j.1365-8711.2002.05528.x
  * Xie et al. (2005) Xie, G. Z., Liu, H. T., Cha, G. W., et al. 2005, AJ, 130, 2506, doi: 10.1086/497163
  * Xie et al. (2008) Xie, G. Z., Yi, T. F., Li, H. Z., Zhou, S. B., & Chen, L. E. 2008, AJ, 135, 2212, doi: 10.1088/0004-6256/135/6/2212
  * Zechmeister & Kürster (2009) Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577, doi: 10.1051/0004-6361:200811296
  * Zhang et al. (2017) Zhang, B. K., Zhao, X. Y., Zhang, L., & Dai, B. Z. 2017, ApJS, 231, 14, doi: 10.3847/1538-4365/aa7ed9