# Transient Quasi-Periodic Oscillations at $\gamma$-rays in the TeV Blazar PKS
1510-089

Abhradeep Roy,1 Arkadipta Sarkar,1 Anshu Chatterjee,1 Alok C. Gupta,2 Varsha
Chitnis,1 and P. J. Wiita3
1Department of High Energy Physics, Tata Institute of Fundamental Research,
Homi Bhabha Road, Mumbai-400005, India
2Aryabhatta Research Institute of Observational Sciences (ARIES), Manora Peak,
Nainital, 263002, India
3Department of Physics, The College of New Jersey, PO Box 7718, Ewing, NJ
08628-0718, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present periodicity search analyses on the $\gamma$-ray lightcurve of the
TeV blazar PKS 1510-089 observed by the Fermi Large Area Telescope. We report
the detection of two transient quasi-periodic oscillations: a 3.6-day QPO
during the outburst in 2009 that lasted five cycles (MJD 54906–54923); and a
periodicity of 92 days spanning over 650 days from 2018 to 2020 (MJD
58200–58850), which lasted for seven cycles. We employed the Lomb-Scargle
periodogram, Weighted Wavelet Z-transform, REDFIT, and the Monte Carlo
lightcurve simulation techniques to find any periodicity and the corresponding
significance. The 3.6-day QPO was detected at a moderate significance of
$\sim$3.5$\sigma$, while the detection significance of the 92-day QPO was
$\sim$7.0$\sigma$. We explore a few physical models for such transient QPOs
including a binary black hole system, precession of the jet, a non-
axisymmetric instability rotating around the central black hole near the
innermost stable circular orbit, the presence of quasi-equidistant magnetic
islands inside the jet, and a geometric model involving a plasma blob moving
helically inside a curved jet.

###### keywords:

Galaxies: individual (PKS 1510-089) – Galaxies: active – Galaxies: jets –
Radiation mechanisms: non-thermal – Gamma rays: galaxies

††pubyear: 2021††pagerange: Transient Quasi-Periodic Oscillations at
$\gamma$-rays in the TeV Blazar PKS 1510-089–Transient Quasi-Periodic
Oscillations at $\gamma$-rays in the TeV Blazar PKS 1510-089

##  1 Introduction

(a)

(b)

(c)

Figure 1: Top: Fermi 12-year weekly aperture photometry light curve. The blue
shaded regions are the epochs where QPO analyses were carried out (i.e., MJD
54890–54935 and MJD 58000–59400). Middle: Fermi 3-hour binned light curve
within MJD 54890–54935. The pink shaded region is the epoch MJD 54906–54923
(EP1) where the final QPO analyses were carried out. Bottom: Fermi 7-day
binned light curve within MJD 58000–59400. The pink shaded region is the epoch
MJD 58200–58850 (EP2) where the final QPO analyses were carried out.

Active Galactic Nuclei (AGN) are known to exhibit highly variable emission
across the electromagnetic spectrum. Radio-loud AGN contain a pair of highly
collimated relativistic plasma jets emanating from the central super-massive
black hole (SMBH). The jets are powered by the accretion process of dense
ionized gases onto the SMBH. Radio-loud AGNs, with jets oriented close to our
line of sight, form a subclass called the blazars (Urry & Padovani, 1995).
Blazars are classified further into BL Lacertae objects (BL Lacs) and flat
spectrum radio quasars (FSRQs) based on the strength of emission lines present
in their optical-UV spectra. BL Lac spectra contain very weak and narrow
emission lines, whereas FSRQs show broad and strong emission lines. The
Doppler-boosted jet radiation dominates the blazar non-thermal emission from
radio to very high energy (VHE) $\gamma$-ray wavebands. Various observational
studies show that the blazars display flux variability of the order of minutes
to years at $\gamma$-ray waveband, as observed by the Fermi-LAT and ground
based atmospheric Cherenkov telescopes (Aleksić et al., 2011; Shukla et al.,
2018).

Although the nature of blazar variability is mostly non-linear, stochastic and
aperiodic (Kushwaha et al., 2017), many studies have claimed detections of
strong quasi-periodic oscillations (QPO) in blazar lightcurves of different
electromagnetic wavebands. In accordance to the variability timescales, the
reported QPOs range from a few tens of minutes to hours to days and even years
of timescales (Raiteri et al., 2001; Liu et al., 2006; Gupta et al., 2009;
Lachowicz et al., 2009; King et al., 2013; Gupta, 2014, 2018; Gupta et al.,
2019; Ackermann et al., 2015; Zhou et al., 2018; Bhatta, 2019; Sarkar et al.,
2020a; Sarkar et al., 2020b, and references therein). However, most of the
QPOs claimed in older studies are marginal detections in that they lasted for
only a 2–4 cycles, while their significances were overestimated (Gupta, 2014).

Continuous monitoring of blazars by Fermi-LAT in the last 12-years has led to
a few recent highly significant $\gamma$-ray QPO detections in different
blazars, including a $\sim$34.5-day transient QPO in PKS 2247-131 (Zhou et
al., 2018), a $\sim$71-day transient QPO in B2 1520+31 (Gupta et al., 2019), a
$\sim$47-day QPO in 3C 454.3 (Sarkar et al., 2020a), a fast periodicity of
$\sim$7.6 days in CTA 102 during an outburst (Sarkar et al., 2020b), and a
periodicity of $\sim$314 days in OJ 287 (Kushwaha et al., 2020a). Although
Covino et al. (2019) did the periodogram analysis on the Fermi-LAT aperture
photometry lightcurves of 10 blazars and claimed the absence of any global
significant periodicity in $\gamma$-rays, a more systematic approach to
$\gamma$-ray QPO detection by Peñil et al. (2020), involving multiple
independent techniques applied on about 2300 AGNs, revealed the presence of
global periodicities of $>4\sigma$ significance in 11 sources, along with 13
more sources with moderately significant (3–4$\sigma$) QPOs.

According to the leptonic models for jet dominated blazar emission, the source
of the radio through optical-UV photons from blazars is the synchrotron
emission by the dense population of ultra-relativistic electrons inside the
magnetized jet. External photon fields from the accretion disk, the broad-line
region (BLR), and the dusty torus can enter the jet. Along with the
synchrotron photons, these external photons get Compton-upscattered by the
same relativistic electron population and produce the high energy
$\gamma$-rays, often up to TeV energies. A few recent studies claimed probably
related QPOs in optical and $\gamma$-rays for several sources, although the
significances were low (Sandrinelli et al., 2016a; Sandrinelli et al., 2016b).
Ackermann et al. (2015) reported a significant $\sim$2.2-year periodicity in
highly-correlated optical and $\gamma$-ray light curves of PG 1553+113.
Several recent studies have reported simultaneous oscillations in both optical
and $\gamma$-ray light curves of the blazars BL Lac (Sandrinelli et al.,
2017), 3C 454.3 (Sarkar et al., 2020a) and CTA 102 (Sarkar et al., 2020b).

PKS 1510-089 (R.A. = 15h 12m 52.2s, Dec. = -09∘ 06’ 21.6") belongs to the FSRQ
subclass of blazars, and it is situated at a cosmological redshift of
$z=0.361$. It is one of the most well-studied blazars and shows high
variability across all electromagnetic wavebands (Zacharias, 2018). The
H.E.S.S. telescope detected PKS 1510-089 in 2010 at TeV energies (H. E. S. S.
Collaboration et al., 2013). It exhibits occasional huge multi-wavelength
outbursts crossing the daily $\gamma$-ray flux level of 10-5 photons cm-2 s-1
(Barnacka et al., 2014; Prince et al., 2017; Meyer et al., 2019), as well as
orphan $\gamma$-ray flaring episodes (Patel et al., 2021). These correlated
multi-wavelength flares can be explained using a shock-in-jet model that
indicates the formation of an emission component in the compact core region,
resulting in optical and $\gamma$-ray flares (Beaklini et al., 2017).
Recently, PKS 1510-089 has been subject to several QPO studies at different
wavelengths. Sandrinelli et al. (2016a) carried out QPO searches on about
3000-day long quasi-simultaneous multi-wavelength lightcurves (MJD
54000–57000) and reported modestly significant periodicities of 115 days
(<3$\sigma$) in $\gamma$-rays, 206 and 490 days ($\sim$3$\sigma$) in the
optical-R band, and 207 and 474 days (>3$\sigma$) in the IR-K band. Wu et al.
(2005) detected periodic deep flux minima of around 1.84 years in its optical
observations of past few years and inferred the presence of a binary black
hole system at the centre. Based on 27-year long UMRAO data of PKS 1510-089 at
4.8, 8 and 14.5 GHz, Fan et al. (2007) claimed a possible $\sim$12-year
periodicity. The 15-year long 22 and 37 GHz radio data from 1990 to 2005
revealed two periodicities of 0.9 years and 1.8 years, which agrees well with
the deep flux minima periodicity of 1.84 years and thereby, strengthens the
binary black hole assumption (Xie et al., 2008; Zhang et al., 2009; Fu et al.,
2014). According to Castignani et al. (2017), QPO analysis of hard X-ray data
of PKS 1510-089 taken by RXTE-PCA from 1996 to 2011 reveals no obvious
periodicity. The most recent QPO study on radio data of PKS 1510-089 spanning
over 38 years detects quite significant ($>4\sigma$) QPOs of $\sim$570,
$\sim$800 and $\sim$1070 days in the 8 GHz and 14.5 GHz light curves (Li et
al., 2021).

In this paper, we report the detection of a probable fast $\gamma$-ray
periodicity of 3.6 days on top of the outburst in 2009 with $>$3.5$\sigma$
significance and a strong $\sim$3-month periodicity during a moderate state of
activity between 2018 and 2020 with a significance above 7$\sigma$. We look
for the plausible scenarios among several QPO models proposed in the
literature to explain these transient QPOs. We consider that the 3.6-day QPO
might have resulted from a hotspot rotating close to the innermost stable
circular orbit around the central SMBH, or from enhanced emission from a few
quasi-equidistant magnetic islands inside the jet. The source of the 92-day
QPO could well be a helical motion of a plasma blob inside a curved jet. We
state the Fermi-LAT data analysis procedure in section 2, and describe the QPO
finding algorithms in section 3. Then we summarize our major results in
section 4. We try to interpret our key results on the basis of various
physical models in section 5 and summarize our conclusions in section 6.

(a)

(b)

Figure 2: left: WWZ map of PKS 1510-089 $\gamma$-ray light curve in the
interval of MJD 54890–-54935. The bright yellowish red patch around the
frequency of 0.3 days-1 indicates the possible presence of a QPO in the
interval of MJD 54900–-54923. right: WWZ map of PKS 1510-089 $\gamma$-ray
lightcurve in the interval of MJD 58000–59400. The bright red patch indicates
a probable QPO in the interval of MJD 58200–58850.

(a)

(b)

(c)

(d)

Figure 3: QPO analysis in EP1 (MJD 54906–54923): Upper-left: WWZ map of PKS
1510-089 $\gamma$-ray light curve during EP1. On the upper-left sub-panel, the
light curve is shown; the lower left panel shows the WWZ map; the lower-right
sub-panel shows the time-averaged WWZ (black) as well as the LSP (red). The
blue dashed line represents the 3.5$\sigma$ significance line against the
power-law red-noise spectrum and the red band on the WWZ map indicates a
strong periodicity of $\sim$3.6 days. Upper-right: Power-spectrum of PKS
1510-089 $\gamma$-ray light curve using REDFIT during EP1. The black line
represents the power-spectrum, the red line is the theoretical AR1 spectrum,
and the yellow, green and blue dashed lines represent FAP levels of 10%, 5%
and 1%, respectively. A strong periodicity of $\sim$3.7 days crosses the 5%
FAP level. Lower-left Flux distribution of the PKS 1510-089 $\gamma$-ray light
curve (PDF in black) fitted with a log-normal model (red) that is used as an
input in light curve simulation. Lower-right Result of light curve simulation
of the PKS 1510-089 $\gamma$-ray light curve during EP1. The black line
represents the LSP of the original light curve and the red line is the mean
LSP of the simulated light curves. The dominant period of $\sim$3.6 days
crosses the 0.01% FAP level (cyan dashed line) and the 3.5$\sigma$
significance curve (green dashed curve).

(a)

(b)

(c)

(d)

Figure 4: QPO analysis in EP2 (MJD 58200–58850): Upper-left: WWZ map of the
PKS 1510-089 $\gamma$-ray light curve during EP2. On the upper-left sub-panel,
the light curve is shown, the lower left sub-panel displays the WWZ map and
the lower-right sub-panel shows the time-averaged WWZ power (black) on top of
the LSP (red). The green and blue dashed curves represent the 5.0$\sigma$ and
7.0$\sigma$ significance respectively against the power-law red-noise spectrum
and the red band on the WWZ map indicates a strong periodicity of $\sim$92
days. Upper-right Power-spectrum of PKS 1510-089 $\gamma$-ray light curve
using REDFIT during EP2. The black line represents the power-spectrum, the red
line is the theoretical AR1 spectrum, and the yellow, green, and blue dashed
lines represent FAP level of 10%, 5%, and 1%, respectively. A strong
periodicity of $\sim$92 days crosses the 1% FAP level. Lower-left: PDF of the
flux distribution of PKS 1510-089 $\gamma$-ray light curve (black) fitted with
a log-normal model (red) that is used as an input in the light curve
simulations. Lower-right: Result of light curve simulations of PKS 1510-089
$\gamma$-ray light curve during EP2. The black line represents the LSP of the
original light curve, and the red line is the mean LSP of the simulated light
curves. The dominant period of $\sim$92 days crosses the 0.2% FAP level (cyan
dashed line) and the 5$\sigma$ significance level (green dashed curve) and
nearly touches the a significance of 7.0$\sigma$ (blue dashed curve).

##  2 Observations and Data Analysis

A detailed recent study on $\gamma$-ray periodicity in AGN by Peñil et al.
(2020) showed that any long-lived QPO in PKS 1510-089 is not significant
(<3$\sigma$). But on visual inspection, the weekly Fermi-LAT aperture
photometry light curve from MJD 58000 to MJD 59400 seemed to have some
periodicity during a long moderate activity state. Moreover, the flare state
light curve modelling by Prince et al. (2017) indicated the presence of fast
periodicity ($\sim$ days) during the flare between MJD 54890 to MJD 54935.

###  2.1 Fermi-LAT Data

We obtained the $\gamma$-ray data from the Large Area Telescope (LAT)
facility, on-board the Fermi observatory. The Fermi-LAT is an imaging space-
based telescope that detects $\gamma$-rays using the pair-production technique
within the 30 MeV to 1 TeV energy range. LAT has a large angular field of view
of about 2.3 sr and covers the entire sky every three hours (Atwood et al.,
2009). We collected the PASS8 (P8R3) processed events’ data of PKS 1510-089
between MJD 54890–54935 and MJD 58000–59400 from the Fermi-LAT data archive
111https://fermi.gsfc.nasa.gov/ssc/data/access/. The PASS8 data provides a
significant improvement in the data quality using an improved reconstruction
of the entire LAT events (Abdollahi et al., 2020).

###  2.2 Data Reduction

We used the standard software package FERMITOOLS-v2.0.8 recommended by the
Fermi-LAT collaboration (Fermi Science Support Development Team, 2019) and the
user-contributed python script ENRICO (Sanchez & Deil, 2013). Following the
recommendations of the Fermi-LAT
collaboration222https://fermi.gsfc.nasa.gov/ssc/data/analysis/, we chose the
events belonging to the SOURCE class (evclass=128, evtype=3) within the energy
range of 0.1–300 GeV from a circular region of interest (ROI) having a radius
of 15∘ centred at the source PKS 1510-089. To get rid of the $\gamma$-ray
contribution from the Earth’s albedo, we selected the events having zenith
angle less than 95∘ followed by the good time interval (GTI) selection using
the standard filter "(DATA$\\_$QUAL$>0$)$\&\&$(LAT$\\_$CONFIG==1)". We
generated an XML file containing the spectral shapes of all the sources lying
within ROI+10∘ radius around the source location according to the fourth
Fermi-LAT (4FGL) catalogue, including the $\gamma$-ray background emission
templates "gll_iem_v07.fits" and "iso_P8R3_SOURCE_V3_v1.txt" for the Galactic
and extra-galactic contributions respectively. We carried out an unbinned
maximum likelihood analysis over the input XML spectral file using the GTLIKE
tool to obtain the source spectrum using the instrumental response function
P8R3_SOURCE_V3. Except for scaling factors, We kept all the spectral
parameters free to vary during the optimization process for the sources lying
within 5∘ from PKS 1510-089. The iterative likelihood analysis removed the
sources having significance less than 1$\sigma$ after each fitting pass. As
mentioned in the 4FGL catalogue, the final source spectrum was modelled using
a log-parabola given as,

$\frac{dN}{dE}=k\left(\frac{E}{E_{b}}\right)^{-\alpha-\beta\log(E/E_{b})},$
(1)

where $\alpha$ is the spectral index at the break energy ($E_{b}$). We kept
$E_{b}$ fixed during the likelihood fitting process.

The full 12-year Fermi unfiltered aperture photometry light curve (1(a)) was
obtained from the Fermi monitoring source list
webpage333https://fermi.gsfc.nasa.gov/ssc/data/access/lat/msl_lc. To extract
the $\gamma$-ray light curve, we divided the whole epoch in a number of time-
bins of our required size and carried out the entire above-mentioned procedure
in each time-bin. For the epoch of MJD 54890–54935, we made a light curve with
3-hour long bins (1(b)), and for the epoch of MJD 58000–59400, we set the bin
size to 7 days (1(c)). For the time-bins where the test-statistics of the flux
estimation were less than 16 (i.e., essentially $<4\sigma$ detection
significance), we estimated flux upper limits at 95% confidence level using
the profile-likelihood method.

##  3 QPO Analysis Methods

On visual inspection, the light curves indicated possible quasi-periodic
modulations. To estimate the time-period of the modulation and the
corresponding significance, we applied four different methods to analyse the
light curves: the Generalized Lomb-Scargle periodogram (GLSP), Weighted
Wavelet Z-transform (WWZ), REDFIT, and light curve simulations. These methods
follow different approaches to detect periodicities and their corresponding
significances in unevenly sampled time-series. Although we have used evenly
binned Fermi-LAT light curves, they became unevenly sampled due to the
omission of flux upper-limits in the analysis processes. The details of these
methods are discussed below.

###  3.1 Generalized Lomb-Scargle Periodogram

The periodogram is one of the most common methods to find periodicities in
light curves, and it gives the power of flux modulations at different
frequencies. For an evenly sampled light curve, the square of the modulus of
its discrete Fourier transform gives the periodogram. But for irregular
sampling, the Lomb-Scargle periodogram (LSP) method iteratively fits sinusoids
with different frequencies to the light curve and constructs a periodogram
from the goodness of the fit (Lomb, 1976; Scargle, 1982). In this work, we
used the Generalized LSP sub-package of the PYASTRONOMY python
package444https://github.com/sczesla/PyAstronomy (Czesla et al., 2019). Unlike
the classical LSP, the GLSP fits a sinusoid plus a constant to the light curve
and takes the errors associated to the measured fluxes into account
(Zechmeister & Kürster, 2009). This code also provides the significance of a
peak in the periodogram in terms of the False Alarm Probability (FAP) given
as,

$FAP(Pn)=1-(1-prob(P>Pn))^{M},$ (2)

where the $FAP$ denotes the probability that at least one out of $M$
independent power values in a given frequency band of a white-noise
periodogram is larger than or equal to the power threshold, $Pn$. In this
work, a peak in a periodogram was considered to be significant when it crossed
the 1% FAP line. The peak position and its corresponding uncertainty were
estimated by fitting a Gaussian to the dominant periodogram-peak. GLSP is an
effective tool to find persistent periodicities. But it cannot usually detect
transient periodicities, as the non-periodic part of the light curve decreases
the goodness of GLSP sinusoid fit. Therefore, the power of the transient
periodicity is reduced.

###  3.2 Weighted Wavelet Z-transform

We used the Weighted Wavelet Z-transform method to detect transient quasi-
periodicities (Foster, 1996). The WWZ method convolves a light curve with a
time- and frequency-dependent kernel and decomposes the data into time and
frequency domains to create a WWZ map. In this work, we used the Morlet kernel
(Grossmann & Morlet, 1984) having the following functional form,

$f[\omega(t-\tau)]=\exp[\iota\omega(t-\tau)-c\omega^{2}(t-\tau)^{2}].$ (3)

Then the WWZ map is given as,

$W[\omega,\tau;x(t)]=\omega^{1/2}\int x(t)f^{\ast}[\omega(t-\tau)]\text{d}t,$
(4)

where $f^{\ast}$ is the complex conjugate of the Morlet kernel $f$, $\omega$
is the frequency, and $\tau$ is the time-shift. This kernel acts as a windowed
discrete Fourier transform having a frequency dependant window size of
$\exp[-c\omega^{2}(t-\tau)^{2}]$. The WWZ map has the advantage of being able
to detect both any dominant periodicities and the time spans of their
persistence.

###  3.3 REDFIT

The REDFIT555https://www.manfredmudelsee.com/soft/redfit/index.htm software
calculates the bias-corrected power-spectrum of a time-series and provides the
significance of the peaks in the spectrum (Schulz & Mudelsee, 2002). REDFIT
fits the light curve with a first-order Autoregressive process (AR1) to
estimate the underlying red-noise spectrum, which is the characteristic of
variable $\gamma$-ray emission from blazars (Covino et al., 2018). The
Autoregressive (AR) models assume that the present observation in a time-
series is related to the past observations. Thus, large fluctuations in the
light curve become less likely (Robinson, 1977). A discrete AR1 process
$\mathscr{F}$ for times $t_{i}$ ($i=1,2,...,N$) with uneven spacing is given
as,

$\displaystyle\mathscr{F}(t_{i})=\theta_{i}\mathscr{F}(t_{i-1})+\epsilon(t_{i}),$
(5) $\displaystyle\theta_{i}=\exp((t_{i-1}-t_{i})/\tau),$

where $\tau$ is the characteristic timescale and $\epsilon$ denotes white-
noise with zero mean. The power-spectrum of an AR1 model has the following
analytical form (Percival & Walden, 1993),

$G_{rr}(f_{j})=G_{0}\frac{1-\theta^{2}}{1-2\theta\cos(\pi
f_{j}/f_{Nyq})+\theta^{2}},$ (6)

where $f_{j}$ denotes the discrete frequency up to the Nyquist frequency
($f_{Nyq}$) and $G_{0}$ is the average spectral amplitude. The “average
autoregression coefficient” ($\theta$) is related to the arithmetic mean of
the sampling intervals $\Delta t=(t_{N}-t_{1})/(N-1)$ as, $\theta=\exp(-\Delta
t/\tau)$, while the $\tau$ comes from the Welch-overlapped-segment-averaging
(WOSA, Welch, 1967) of the LSP.

REDFIT estimates the significance of the peaks in the power-spectrum using FAP
level up to the minimum of 1%.

###  3.4 Lightcurve Simulation

Another way to estimate the peak significance in a periodogram is to simulate
light curves using a Monte Carlo method, following the power spectral density
(PSD) and the flux distribution (PDF) of the original light curve
(Emmanoulopoulos et al., 2013). Simple power-laws ($P\propto\nu^{-\alpha}$)
give reasonably good approximations of the underlying red-noise PSDs of blazar
lightcurves (Vaughan, 2005). Hence, we used a power-law and a log-normal model
to respectively fit the PSD and PDF of the original light curve. The log-
normal model for PDF has the form,

$PDF(x)=\frac{1}{x\sigma\sqrt{2\uppi}}\exp\left[-\frac{(\ln
x-\mu)^{2}}{2\sigma^{2}}\right].$ (7)

Then we simulated 2000 light curves with the fitted PSD and PDF models as
inputs using the
DELightcurveSimulation666https://github.com/samconnolly/DELightcurveSimulation
code (Connolly, 2016). The mean and standard deviation of the simulated light
curve GLSP at each frequency allowed us to estimate the significance of the
dominant periodicities.

##  4 Results

1(a) shows the epochs on the weekly Fermi-LAT full aperture photometry
lightcurve that were initially considered for QPO analysis. But after
analysing the Fermi-LAT data, we selected two sub-intervals (1(b) and 1(c)) to
carry on detailed QPO analysis: MJD 54906–54923 (EP1) and MJD 58200–58850
(EP2).

###  4.1 EP1: 16 March 2009 – 2 April 2009 (MJD 54906–54923)

From 2(a) it is evident that there is a QPO in MJD 54900–54923 around the
frequency of 0.3 days-1. But in 1(b), it can be clearly seen that there is a
big gap in the data and in addition, the available flux points are mostly
upper-limits during MJD 54890–54906. So, we avoided these intervals and
carried out our QPO analyses in EP1 (MJD 54906–54923).

3(a) shows the WWZ map of PKS 1510-089 during EP1. The strong horizontal red
patch denoting the 3.6-day QPO spans the entire 18-day EP1 light curve,
indicating the presence of 5 cycles. The time-averaged WWZ plot indicates a
QPO of 3.6 days with 3.5$\sigma$ significance. 3(b) shows the power-spectrum
and the corresponding peak significance obtained using REDFIT. It indicates a
QPO of $\sim$3.7 days, crossing the 5% FAP level. 3(c) represents the log-
normal fitted flux distribution in this epoch. 3(d) shows the LSP of PKS
1510-089 $\gamma$-ray light curve during EP1 and the result of significance
check using our light curve simulation. The dominant periodicity of
3.63${}^{+0.07}_{-0.07}$ days crosses the 0.01% FAP level. We note that this
FAP significance estimation only holds for Gaussian random noise, whereas
blazars show red-noise type variability. However, the light curve simulation,
which does not rely on that noise assumption, indicates that the QPO of 3.6
days has a significance of 3.5$\sigma$.

In recent studies (Kushwaha et al., 2020b; Sarkar et al., 2020a; Peñil et al.,
2020; Sarkar et al., 2020b), the reported blazar QPOs generally have
$>$3$\sigma$ significance and cross the REDFIT FAP level of 1%. In this case,
the QPO spans only 18 days and the apparent periodicity is very fast
($\sim$3.6 days) which should be the reason for low significance in the REDFIT
output. But from the light curve simulation, we find a $\sim$3.5$\sigma$
significance, which makes the QPO sufficiently significant to be reported. The
QPO peak significances remained the same when we applied the same analysis
procedures on light curves in EP1 with different bin-sizes, such as 4 hours or
5 hours.

###  4.2 EP2: 23 March 2018 – 2 January 2020 (MJD 58200–58850)

1(c) shows the 7-day binned Fermi-LAT lightcurve between MJD 58000 to MJD
59400. Although the GLSP did not reveal any significant dominant QPO in this
interval, the WWZ map showed a bright red patch with a span of $\sim$650 days
from MJD 58200 to MJD 58850 (2(b)). This led to the selection of EP2 (MJD
58200–58850) for detailed analysis.

4(a) shows the WWZ map of PKS 1510-089 during EP2. The strong horizontal red
patch, spanning all the 650-day EP2 light curve, denotes a strong 92-day QPO,
indicating the presence of 7 cycles. The time-averaged WWZ also show a QPO of
92 days, with almost 7.0$\sigma$ significance. 4(b) shows the power-spectrum
and the corresponding peak significance obtained using REDFIT. It indicates a
QPO of $\sim$92 days, crossing the 1% FAP level. 4(c) represent the log-normal
fitted flux distribution during EP2. 4(d) shows the LSP of PKS 1510-089
$\gamma$-ray light curve during EP2 and the significance check of the observed
QPO using light curve simulation. The dominant periodicity of
91.5${}^{+1.2}_{-1.2}$ days crosses the 0.2% FAP level, if one considers an
underlying Gaussian-type noise. Light curve simulation shows that the QPO of
92 days seems to have a significance of $>6\sigma$. The 92-day peak touches
the blue dashed 7.0$\sigma$ significance line. Thus, all the applied tools
indicate this 92-day QPO to be highly significant. We have also checked that
our key results remain the same when light curves produced when we employed
different bin-sizes of 3 days or 5 days during EP2.

##  5 Discussion:

We obtained $\gamma$-ray data from the Fermi-LAT archive and analysed then to
generate the light curves in two different epochs: EP1 (MJD 54906–54923) and
EP2 (MJD 58200–58850). We employed LSP, WWZ and REDFIT methods to detect
significant transient periodicities in the light curves. We generated 2000
light curves in each epoch using Monte-Carlo simulation to account for the
underlying red-noise spectrum while estimating the QPO significances. All of
these methods revealed two highly probable transient QPOs: (1) a fast QPO of
3.6 days in EP1 and; (2) a 92-day QPO in EP2.

Due to lack of good coverage of optical observations during the selected
epochs, it was not possible to employ QPO analysis tools on the available
optical data. In case of an FSRQ, both synchrotron emission from the jet and
thermal emission from the accretion disk contribute to the optical-UV
emission. The $\gamma$-ray emission is dominated by the inverse-Compton
scattering of seed photons by the charged particle population inside the jet
responsible for the synchrotron emission. Thus, simultaneous optical and
$\gamma$-ray QPO can lead to strong inferences about possible reasons behind
such phenomena and the underlying disk-jet connection. Only having
$\gamma$-ray data here, it is quite hard to conclusively investigate the
probable physical reasons behind such periodicities. Still, the $\gamma$-ray
QPO timescale and persistence can lead to a few insights about the jet
structure or emission processes.

Earlier studies have proposed several possible models, such as supermassive
binary black-hole systems (Valtonen et al., 2008; Villforth et al., 2010;
Ackermann et al., 2015; Li et al., 2021), persistent jet precession model
(Romero et al., 2000; Rieger, 2004; Liska et al., 2018), and Lense-Thirring
precession of accretion disks (Stella & Vietri, 1998; Liska et al., 2018) to
explain the long-term quasi-periodicities in different blazar’s emissions.
Although PKS 1510-089 is a candidate to contain a binary black-hole system at
the centre (Wu et al., 2005; Li et al., 2007), all these models exhibit
persistent QPOs with at least year-long periods. So, we can probably discard
these models as an explanation for the QPOs discussed here.

Optical QPOs with tens of days period can be explained by hot-spots rotating
at or near the innermost stable circular orbit (ISCO) around the central SMBH
(Zhang & Bao, 1991; Gupta et al., 2009). In this case, the rotation of the
hot-spot will modulate the seed photon field of the external inverse-Compton
scattering inside the jet, causing a modulation in the $\gamma$-ray emission.
The $\gamma$-ray emission from blazar jet is Doppler boosted and should have a
faster quasi-periodicity ($\sim$days) in this model. Thus, rotation of a hot-
spot near the ISCO could be a possible scenario behind the $\sim$3.6-day QPO
on top of the flare during EP1. Availability of an optical counterpart of this
short QPO could actually emphasize the applicability of this model in EP1.
Assuming that the QPO is related to orbital rotation of a hotspot, presence of
a spiral shock or any other non-axisymmetric instabilities close to the ISCO,
we can estimate the central SMBH mass using an expression given by Gupta et
al. (2009),

$\frac{M_{\text{BH}}}{M_{\sun}}=\frac{3.23\times 10^{4}P}{(r^{3/2}+a)(1+z)},$
(8)

where $P$ is the QPO period in seconds, $z$ is the cosmological redshift of
the source ($z$ = 0.361 for PKS 1510-089), $r$ is the radius of ISCO in units
of $GM_{\text{BH}}/c^{2}$, and $a$ is the SMBH spin parameter.

For a Schwarzschild black hole, $r=6.0$ and $a=0$, and for a maximally
rotating Kerr black hole, $r=1.2$ and $a=0.9982$. Using the period of $\sim
3.6$ days, the estimated mass of a Schwarzschild BH is $\sim$5.0$\times
10^{8}M_{\sun}$, and the mass of a maximally rotating Kerr BH is
$\sim$3.2$\times 10^{9}M_{\sun}$. While using the period of 92 days, we get
the mass estimate of Schwarzschild BH to be $\sim$1.2$\times 10^{10}M_{\sun}$,
and that of a Kerr BH to be $\sim$8.2$\times 10^{10}M_{\sun}$. There are
various methods to estimate the mass of the central SMBH. Gu et al. (2001)
used the data of the Hβ line width and the optical continuum luminosity and
reported a mass of $\sim$1.3$\times 10^{9}M_{\sun}$ for the central SMBH of
PKS 1510-089. The reverberation mapping technique is one of the most accurate
methods to estimate the mass of a primary black hole. Rakshit (2020) reported
a mass of $\sim$5.71$\times$10${}^{7}M_{\sun}$ for the SMBH of PKS 1510-089
using a stereoscopic reverberation mapping technique. Xie et al. (2005) used
reverberation mapping and short-timescale optical variability to estimate the
SMBH mass to be $\sim$2.0$\times 10^{8}M_{\sun}$ and $\sim$1.6$\times
10^{8}M_{\sun}$, respectively. The masses estimated using the 92-day QPO are
much higher, thereby essentially excluding the possibility that disk
instabilities directly yield it. But the masses estimated using the 3.6-day
QPO are comparable to the SMBH mass estimated by Gu et al. (2001). If this is
the case, it favours the presence of a non-maximally rotating SMBH at the
centre. However, this rotating hot-spot scenario suffers a major disadvantage
in explaining blazar fluctuations. Blazar disks have almost face-on
orientation with respect to the observer, so the motion of the rotating
hotspot should be almost azimuthally symmetric, which implies that this
situation is unlikely to generate enough flux variability. Rotating hotspots
are more likely to generate QPOs if the observer’s line of sight is close to
the plane of the accretion disk (Rani et al., 2009).

Another interesting model that might explain fast quasi-periodicity in the jet
emission involves magnetic reconnection events in almost equidistant magnetic
islands inside the jet (Huang et al., 2013). These equispaced magnetic islands
periodically enhance the flux, thereby producing a rapid transient QPO. Shukla
et al. (2018) used such a magnetic reconnection process to model the extremely
fast variability ($\sim$5 minutes) in the FSRQ CTA 102 during its outburst in
2016. This model does appear to be capable of producing the 3.6-day QPO in
$\gamma$-rays for PKS 1510-089, and suffers no obvious difficulties.

We can also attempt to attribute both the 3.6-day and 92-day QPO of PKS
1510-089 to a very reasonable model that involves a plasma blob moving
helically down the jet (Mohan & Mangalam, 2015; Sobacchi et al., 2016; Sarkar
et al., 2020a). For the simplest leptonic one-zone model (where the bulk of
the synchrotron and inverse-Compton emission comes from a single region), such
a plasma blob contains higher particle and magnetic energy densities and is
responsible for occasional enhanced emission from blazars. Due to the
postulated helical motion of the blob, the viewing angle of the blob with
respect to our line of sight ($\theta_{\text{obs}}$) changes periodically with
time as,

$\cos\theta_{\text{obs}}(t)=\sin\phi\sin\psi\cos(2\pi
t/P_{\text{obs}})+\cos\phi\cos\psi,$ (9)

where $\phi$ is the pitch angle of the helical path, $\psi$ is the angle of
the jet axis with respect to our line of sight, and $P_{\text{obs}}$ is the
observed periodicity in the light curve (Sobacchi et al., 2016; Zhou et al.,
2018). The Doppler factor ($\delta$) varies with the viewing angle as
$\delta=1/[\Gamma(1-\beta\cos\theta_{\text{obs}})]$, where
$\Gamma=1/\sqrt{1-\beta^{2}}$ is the bulk Lorentz factor of the blob motion
with $\beta=v_{\text{jet}}/c$. Then the periodicity in the blob rest frame is
given as,

$P_{\text{rf}}=\frac{P_{\text{obs}}}{1-\beta\cos\psi\cos\phi}.$ (10)

This model can naturally explain the transient nature of any periodicities as
the QPO starts when the blob is injected into the jet and lasts until the blob
dissipates. One limitation of this model is that it can only explain a QPO
having almost constant amplitude. It is evident from 3(a) and 4(a) that both
the QPOs in EP1 and EP2 have varying amplitudes. We note that Sarkar et al.
(2020a) used a curved jet scenario to model a transient QPO with varying
amplitude in 3C 454.3. In this model, the angle between the jet axis and our
line of sight ($\psi$) becomes time-dependent ($\psi(t)$). We find that the
"blob moving helically in a curved jet" model might be able to explain them
both. Assuming $\phi\simeq 2^{\circ}$ (Zhou et al., 2018),
$\langle\psi\rangle=2^{\circ}.2$, $\Gamma=20.0$ (Roy et al., 2021) and
$P_{\text{obs}}=3.6$ days, the periodicity in the blob rest-frame is
$P_{\text{rf}}\simeq 3.8$ years. The blob traverses about a distance $D=c\beta
P_{\text{rf}}\cos\phi\simeq 1.16$ parsec down the jet during one period. But
only a very high jet curvature can explain the rapid changes in QPO amplitudes
in EP1. However, the presence of high curvature within a few parsecs is highly
unlikely in an extremely long, well-collimated, and powerful FSRQ jet. On the
other hand, for $P_{\text{obs}}=92$ days, the periodicity in the blob rest-
frame is $P_{\text{rf}}\simeq 97.1$ years. The blob travels $\sim$30 pc during
one period, i.e., it travels $\sim$200 pc during EP2. 4(a) shows a slow
increasing trend in the first 5 cycles and then a faster attenuation in the
QPO amplitudes. This situation can be explained with a much lower curvature in
the jet, which is physically more likely. Thus, we tentatively conclude that
the 3.6-day QPO probably resulted from flux enhancements by magnetic
reconnection events at quasi-equispaced magnetic islands, while most probably
the origin of the 92-day QPO is a plasma blob moving helically inside a curved
jet.

##  6 Conclusions

In this work, we report two transient QPOs in $\gamma$-rays exhibited by the
TeV blazar PKS 1510-089. We employed several standard tools such as Lomb-
Scargle periodogram, REDFIT, WWZ and lightcurve simulation on the Fermi-LAT
light curve to detect the significant periodicities. Our key results are as
follows.

  1. 1.

PKS 1510-089 showed a fast periodicity of $\sim 3.6$ days during a flare in
2009, from MJD 54906 to MJD 54923. It lasted for only 18 days. Lightcurve
simulation indicates the significance of this QPO to be $\sim 3.5\sigma$
against the underlying red-noise spectrum. To our knowledge, this is the
shortest period QPO so far reported in the $\gamma$-ray emission of a blazar.

  2. 2.

Multiple QPO analyses show the presence of a QPO around 92 days in the recent
Fermi-LAT observations that lasted for about 650 days from 2018 to 2020 (MJD
58200–58850). From light curve simulation, it seems that this periodicity has
a significance of about 7$\sigma$, and would be the most significant blazar
QPO ever reported.

## Acknowledgements

This research has made use of the Fermi-LAT data, obtained from the Fermi
Science Support Center, provided by NASA’s Goddard Space Flight Center (GSFC).
The Fermi LAT Collaboration acknowledges generous on-going support from a
number of agencies and institutes that have supported both the development and
the operation of the LAT as well as scientific data analysis. These include
the National Aeronautics and Space Administration and the Department of Energy
in the United States, the Commissariat l’ Energie Atomique and the Centre
National de la Recherche Scientifique / Institut National de Physique
Nucléaire et de Physique des Particules in France, the Agenzia Spaziale
Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry
of Education, Culture, Sports, Science and Technology (MEXT), High Energy
Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency
(JAXA) in Japan, and the KȦ. Wallenberg Foundation, the Swedish Research
Council and the Swedish National Space Board in Sweden. Additional support for
science analysis during the operations phase is gratefully acknowledged from
the Istituto Nazionale di Astrofisica in Italy and the Centre National d’
Études Spatiales in France.The data and analysis software were obtained from
NASA’s High-Energy Astrophysics Science Archive Research Center (HEASARC), a
service of GSFC. We used a community-developed Python package named ENRICO to
make Fermi-LAT data analysis easier (Sanchez & Deil, 2013). Finally, We
acknowledge the support of the Department of Atomic Energy, Government of
India, under Project Identification No. RTI 4002.

## Data Availability

  1. 1.

The Fermi-LAT data used in this article are available in the LAT data server
at https://fermi.gsfc.nasa.gov/ssc/data/access/.

  2. 2.

The Fermi-LAT data analysis software is available at
https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/.

  3. 3.

We agree to share data derived in this article on reasonable request to the
corresponding author.

## References

  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Ackermann et al. (2015) Ackermann M., et al., 2015, ApJ, 813, L41
  * Aleksić et al. (2011) Aleksić J., et al., 2011, ApJ, 730, L8
  * Atwood et al. (2009) Atwood W. B., et al., 2009, ApJ, 697, 1071
  * Barnacka et al. (2014) Barnacka A., Moderski R., Behera B., Brun P., Wagner S., 2014, A&A, 567, A113
  * Beaklini et al. (2017) Beaklini P. P. B., Dominici T. P., Abraham Z., 2017, A&A, 606, A87
  * Bhatta (2019) Bhatta G., 2019, MNRAS, 487, 3990
  * Castignani et al. (2017) Castignani G., et al., 2017, A&A, 601, A30
  * Connolly (2016) Connolly S. D., 2016, DELightcurveSimulation: Light curve simulation code (ascl:1602.012)
  * Covino et al. (2018) Covino S., Sandrinelli A., Treves A., 2018, MNRAS, 482, 1270
  * Covino et al. (2019) Covino S., Sandrinelli A., Treves A., 2019, MNRAS, 482, 1270
  * Czesla et al. (2019) Czesla S., Schröter S., Schneider C. P., Huber K. F., Pfeifer F., Andreasen D. T., Zechmeister M., 2019, PyA: Python astronomy-related packages (ascl:1906.010)
  * Emmanoulopoulos et al. (2013) Emmanoulopoulos D., McHardy I. M., Papadakis I. E., 2013, MNRAS, 433, 907
  * Fan et al. (2007) Fan J. H., et al., 2007, A&A, 462, 547
  * Fermi Science Support Development Team (2019) Fermi Science Support Development Team 2019, Fermitools: Fermi Science Tools (ascl:1905.011)
  * Foster (1996) Foster G., 1996, AJ, 112, 1709
  * Fu et al. (2014) Fu J. P., Zhang H. J., Zhang X., Xiong D. R., Guo F., 2014, AcASn, 55, 1
  * Grossmann & Morlet (1984) Grossmann A., Morlet J., 1984, SJMA, 15, 723
  * Gu et al. (2001) Gu M., Cao X., Jiang D., 2001, MNRAS, 327, 1111
  * Gupta (2014) Gupta A. C., 2014, JApA, 35, 307
  * Gupta (2018) Gupta A., 2018, Galaxies, 6, 1
  * Gupta et al. (2009) Gupta A. C., Srivastava A. K., Wiita P. J., 2009, ApJ, 690, 216
  * Gupta et al. (2019) Gupta A. C., Tripathi A., Wiita P. J., Kushwaha P., Zhang Z., Bambi C., 2019, MNRAS, 484, 5785
  * H. E. S. S. Collaboration et al. (2013) H. E. S. S. Collaboration et al., 2013, A&A, 554, A107
  * Huang et al. (2013) Huang C.-Y., Wang D.-X., Wang J.-Z., Wang Z.-Y., 2013, RAA, 13, 705
  * King et al. (2013) King O. G., et al., 2013, MNRAS, 436, L114
  * Kushwaha et al. (2017) Kushwaha P., Sinha A., Misra R., Singh K. P., de Gouveia Dal Pino E. M., 2017, ApJ, 849, 138
  * Kushwaha et al. (2020a) Kushwaha P., Sarkar A., Gupta A. C., Tripathi A., Wiita P. J., 2020a, MNRAS, 499, 653
  * Kushwaha et al. (2020b) Kushwaha P., Sarkar A., Gupta A. C., Tripathi A., Wiita P. J., 2020b, MNRAS, 499, 653
  * Lachowicz et al. (2009) Lachowicz P., Gupta A. C., Gaur H., Wiita P. J., 2009, A&A, 506, L17
  * Li et al. (2007) Li J., Fan J.-H., Yuan Y.-H., 2007, Chinese Physics, 16, 876
  * Li et al. (2021) Li X.-P., Zhao L., Yan Y., Wang L.-S., Yang H.-T., Cai Y., Luo Y.-H., 2021, JApA, 42, 92
  * Liska et al. (2018) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2018, MNRAS, 474, L81
  * Liu et al. (2006) Liu F. K., Zhao G., Wu X.-B., 2006, ApJ, 650, 749
  * Lomb (1976) Lomb N. R., 1976, Ap&SS, 39, 447
  * Meyer et al. (2019) Meyer M., Scargle J. D., Blandford R. D., 2019, ApJ, 877, 39
  * Mohan & Mangalam (2015) Mohan P., Mangalam A., 2015, ApJ, 805, 91
  * Patel et al. (2021) Patel S. R., Bose D., Gupta N., Zuberi M., 2021, JHEAp, 29, 31
  * Peñil et al. (2020) Peñil P., et al., 2020, ApJ, 896, 134
  * Percival & Walden (1993) Percival D. B., Walden A. T., 1993, Spectral Analysis for Physical Applications. Cambridge University Press, doi:10.1017/CBO9780511622762
  * Prince et al. (2017) Prince R., Majumdar P., Gupta N., 2017, ApJ, 844, 62
  * Raiteri et al. (2001) Raiteri C. M., et al., 2001, A&A, 377, 396
  * Rakshit (2020) Rakshit S., 2020, A&A, 642, A59
  * Rani et al. (2009) Rani B., Wiita P. J., Gupta A. C., 2009, ApJ, 696, 2170
  * Rieger (2004) Rieger F. M., 2004, ApJ, 615, L5
  * Robinson (1977) Robinson P. M., 1977, Stochastic Processes and their Applications, 6, 9
  * Romero et al. (2000) Romero G. E., Chajet L., Abraham Z., Fan J. H., 2000, A&A, 360, 57
  * Roy et al. (2021) Roy A., Patel S. R., Sarkar A., Chatterjee A., Chitnis V. R., 2021, MNRAS, 504, 1103
  * Sanchez & Deil (2013) Sanchez D. A., Deil C., 2013, Enrico : a Python package to simplify Fermi-LAT analysis (arXiv:1307.4534)
  * Sandrinelli et al. (2016a) Sandrinelli A., Covino S., Dotti M., Treves A., 2016a, AJ, 151, 54
  * Sandrinelli et al. (2016b) Sandrinelli A., Covino S., Treves A., 2016b, ApJ, 820, 20
  * Sandrinelli et al. (2017) Sandrinelli A., et al., 2017, A&A, 600, A132
  * Sarkar et al. (2020a) Sarkar A., Gupta A. C., Chitnis V. R., Wiita P. J., 2020a, MNRAS, 501, 50
  * Sarkar et al. (2020b) Sarkar A., Kushwaha P., Gupta A. C., Chitnis V. R., Wiita P. J., 2020b, A&A, 642, A129
  * Scargle (1982) Scargle J. D., 1982, ApJ, 263, 835
  * Schulz & Mudelsee (2002) Schulz M., Mudelsee M., 2002, Computers & Geosciences, 28, 421
  * Shukla et al. (2018) Shukla A., et al., 2018, ApJ, 854, L26
  * Sobacchi et al. (2016) Sobacchi E., Sormani M. C., Stamerra A., 2016, MNRAS, 465, 161
  * Stella & Vietri (1998) Stella L., Vietri M., 1998, ApJ, 492, L59
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, PASP, 107, 803
  * Valtonen et al. (2008) Valtonen M. J., et al., 2008, Nature, 452, 851
  * Vaughan (2005) Vaughan S., 2005, A&A, 431, 391
  * Villforth et al. (2010) Villforth C., et al., 2010, MNRAS, 402, 2087
  * Welch (1967) Welch P., 1967, IEEE Transactions on Audio and Electroacoustics, 15, 70
  * Wu et al. (2005) Wu J., Zhou X., Peng B., Ma J., Jiang Z., Chen J., 2005, MNRAS, 361, 155
  * Xie et al. (2005) Xie G. Z., Liu H. T., Cha G. W., Zhou S. B., Ma L., Xie Z. H., Chen L. E., 2005, AJ, 130, 2506
  * Xie et al. (2008) Xie G. Z., Yi T. F., Li H. Z., Zhou S. B., Chen L. E., 2008, AJ, 135, 2212
  * Zacharias (2018) Zacharias M., 2018, in High Energy Astrophysics in Southern Africa (HEASA2018). p. 33 (arXiv:1903.02274)
  * Zechmeister & Kürster (2009) Zechmeister M., Kürster M., 2009, A&A, 496, 577
  * Zhang & Bao (1991) Zhang X. H., Bao G., 1991, A&A, 246, 21
  * Zhang et al. (2009) Zhang H., Zhao G., Zhang X., Dong F., Xie Z., Yi T., Zheng Y., Yu Y., 2009, ScChG, 52, 1442
  * Zhou et al. (2018) Zhou J., Wang Z., Chen L., Wiita P. J., Vadakkumthani J., Morrell N., Zhang P., Zhang J., 2018, NatCo, 9, 4599