Study of Periodic Signals from Blazars †

: The search for periodic signals from blazars has become an actively pursued ﬁeld of research in recent years. This is because periodic signals bring us information about the processes occurring near the innermost regions of blazars, which are mostly inaccessible to our direct view. Such signals provide insights into some of the extreme conditions that take place in the vicinity of supermassive black holes that lead to the launch of the relativistic jets. In addition, studies of characteristic timescales in blazar light curves shed light on some of the challenging issues in blazar physics that include disk-jet connection, strong gravity near fast-rotating supermassive black holes and release of gravitational waves from binary supermassive black hole systems. However, a number of issues associated with the search for quasi-periodic oscillations (QPOs) in blazars e.g., red-noise dominance, modest signiﬁcance of the detection, periodic modulation lasting for only a couple of cycles and their transient nature, make it difﬁcult to estimate the true signiﬁcance of the detection. Consequently, it also becomes difﬁcult to make meaningful inferences about the nature of the on-going processes. In this proceedings, results of study focused on searching for QPOs in a number of blazar multi-frequency light curves are summarized. The time series analyses of long term observations of the blazars revealed the presence of year-timescale QPOs in the sources including OJ 287 (optical), Mrk 501 (gamma-ray), J1043+2408 (radio) and PKS 0219-164 (radio). A likely explanations, we discuss a number of scenarios including binary supermassive black hole systems, lense-thirring precession, and jet precession. 2 → weighted Wavelet transform, ACF → auto-correlation function, and EF → epoch folding.


Introduction
Active galactic nuclei (AGN) are the most luminous sources in the Universe. As accreting matter swirls inward to the supermassive black hole at its center [1], the gravitational potential energy of the material is processed into radiation energy resulting in a large output of high energy emission. As a result, the central core, mainly the accretion disk with black hole at the center, outshines the whole galaxy. A small fraction of AGN happen to eject powerful jets aligned along the line of sight. These objects are known as blazars which are dominant sources of relativistically boosted non-thermal sources in the universe. Blazars consist of flat-spectrum radio quasars (FSRQ), some of the most powerful sources, and BL Lacertae (BL Lac) objects, cosmic sources of highest energy emission. As blazars are bright and visible across cosmological distances, they could be very useful in mapping the large scale structure of the Universe. More recently, blazar TXS 0506+056 was linked to the first non-stellar neutrino emission detected by the IceCube experiment [2].
Blazars display flux variations on diverse spatial and temporal frequencies see [3] and references therein. The flux variability shown by AGN mostly appears aperiodic, but presence of quasi-periodic oscillations (QPO) in the multi-frequency blazar (both flat spectrum radio quasars and BL Lacertae objects) light curves has been frequently reported, especially owing to the availability of long term observations from various instruments in a wide range of frequency bands e.g optical radio and γ-ray see [4,5] and references therein. As a result, the hunt for periodic signals in blazar light curves has become an active field of research: blazar OJ 287 is famous for its ∼12-year periodicity detected in its historical optical light curve [6,7] and another famous source PG 1553 + 113 is found to show 2. 18-year periodic modulations in its γ-ray observations from the Fermi/LAT, which was first reported in [8]. In addition, several authors have reported the existence of QPOs of various characteristic timescales in numerous blazar multi-frequency light curves including the four cases of QPO detection by our research group see [4,[9][10][11] and the references therein.
QPO signals mostly likely originate at the innermost regions of blazars and propagate along the jets to reach us. Therefore, the signals carry crucial information about the processes occurring near the central engine. Performing an in-depth time series analysis of blazar light curves can provide insights into the processes near the AGN central engines and thereby add to our current understanding of the physics of supermassive black holes and origin relativistic jets. In these proceedings, I summarize our results of study of QPO analysis in blazars, and also discuss some of the methods that are mostly used to detect possible periodic components in blazar multi-frequency light curves and establish significance of the detection using Monte Carlo (MC) simulations. Besides, some of the AGN scenarios that could possibly result in QPOs are explored.

Analysis Methods
Some of the popular methods of time series analysis that are employed to search for periodicity in astronomical observations are discussed below.

Discrete Fourier Periodogram
For a given time series x(t j ) sampled at times t j with j = 1, 2, ..., N, the discrete Fourier power at a temporal frequency ν is estimated using the relation where T andx represent the total duration of the series, and the mean flux of a variable source, respectively [12]. If the observations sampled at discrete time step ∆t span a total duration of T, the powers are computed at the evenly spaced log-frequency grid between the lowest frequency 1/T to the highest 1/(2∆t), also known as Nyquist frequency. Thus computed discrete power is of ten known as discrete Fourier periodogram (DFP). The distribution of DFP over the considered temporal frequencies reveal variability power of the source light curve at the corresponding timescales, and thereby provide crucial information about the underlying variability structures and dominant timescales. In particular, if a significant periodogram peak centered around a frequency is observed, it is most likely due to the oscillations in the data that repeat after the characteristic timescale corresponding to the peak frequency. The computational simplicity makes the method first choice among astronomers searching for periodic signals in astronomical observations. However, the results of the method are largely affected by the irregularities and gaps in the data.

Lomb-Scargle Periodogram
While the DFP of a unevenly sampled finite light curve can produce many artifacts which could be mistaken as periodic components, the Lomb-scargle periodogram LSP; [13,14] is considered to be an improved method to compute the periodograms. For an angular frequency, ω, the LSP is given as where τ is given by The method attempts a fit to sine waves of the form X f (t) = A cos ωt + B sin ωt to the data. Therefore it is less affected by the gaps and irregularity in the sampling of data when compared to DFP. Besides, the least-square fitting process enhances the sinusoidal component in the periodogram. This makes the method more efficient in detecting periodic signals in light curves. Usually the periodogram is computed for the minimum and maximum frequencies of , ν min = 1/T, and ν max = 1/(2∆t), respectively. The total number of frequencies considered N ν also plays an important role in the evaluation of the periodogram. It is empirically given as N ν = n 0 Tν max , where n 0 can be chosen in the range of 5-10 see [15]. As with the DFP, a distinct peak in the LSP signifies the possible presence of an underlying periodic signal.

Weighted Wavelet Z-Transform
The DFP and LSP are efficient in detecting periodic signals especially when the signals span the entire observation duration. However, in real astronomical systems, periodic oscillations may develop and evolve in frequency and amplitude over time, hence quasi-periodic signals. In such cases, time series analysis based on the wavelet transforms are much more useful. Wavelet methods attempt to fit sinusoidal waves that are localized in both time and frequency domains. In the context of astronomical observations where data are full of irregularities and gaps, weighted wavelet z-transform (WWZ) is one of the methods widely applied. In this method, the wavelet functions can be viewed as weighted projections on the trial functions given as which further are weighted with a weight function w i = e −c ω 2 (t i −τ) 2 , where c ∼ 0.0125 is a fine tuning parameter. Then WWZ power can be written in terms of weighted variations of the data, V x and V y , given as and where N e f f is the effective number of the data points see [16] for further details. When WWZ analysis is performed, the distribution of color-scaled WWZ power of the source in the time-period (Note that the time and the timescale corresponding to a temporal frequency are two different measures.) plane can reveal QPOs e.g., [4,9]. In case where the data contains a variable period the maximum powers are seen to "meander" on the time-period plane.

Epoch Folding
In contrast to the DFP, LSP, and WWZ, which represent frequency domain based methods of time series analysis, epoch folding (see [17][18][19]) is a method which is based on time domain. The method has several advantages over traditional methods: it is capable of detecting the periodic signals of any shapes, it is relatively less affected to the uneven sampling and gaps in the observations, and it is free of all the artifacts associated with frequency domain analysis such as aliasing and red-noise leak. For periodicity tests, a given light curve is folded on a range of trial periods and phase bins, and consequently, χ 2 is computed according to where x i and σ i represent the mean and standard deviation, respectively, of each of M phase bins.
In the absence of any periodic signal in the data, or the observations distributed areas Gaussian noise, we find χ 2 ∼ M. However in the case when a light curve contains a periodic signal, χ 2 at the characteristic timescale becomes large and significantly different from it's mean value for details refer to [20]. χ 2 values are computed for the trial periods in the range of suspected timescales. The pulse prof iles obtained for the trial periods are tested using Equation (6). Usually, the maximum χ 2 deviation at a timescale signifies the underlying most probable period in the data.

Discrete Auto-Correlation Function
The discrete auto-correlation function (ACF) is another time domain based method frequently applied to search periods in astronomical time series. The method is based on the discrete correlation function (DCF), a method less affected to the sampling irregularity in the observations. The unbinned DCF can be computed as, wherex andȳ represent mean values of the two time series and σ 2 , e 2 correspond to variance and uncertainties in the light curves, respectively see [21]. When these discrete pairs are binned in a suitable time lag bin, the mean DCF of the bin including the M pairs is expressed as, However, in this method, the sampling distribution of DCF are of ten highly skewed such that the measure of DCF uncertainties represented by sample variances can be unreliable. In such a case, the DCFs can be z-transformed (i.e., ZDCFs) so that the distribution closely resembles normal distribution, and then standard deviation provides a robust estimation of ZDCF uncertainties for details see [22]. Now to compute ACF we calculate ZDCF using the same light curve in place of two time series data. If a time series observations contain periodic signals, the distribution of ACF over the time lag also shows periodic oscillations see [23] and the reference therein.
In all the of methods described above, the half-width at the half maximum of the Gaussian fit to the central peak, which signified the presence of QPO, was used as a measure of the uncertainty in the detected period.

Significance Estimation and Monte Carlo Simulation
An observed light curve x(t i ), sampled at discrete t i ( where i = 1 through N) times can be viewed as the a continuous light curve x(t) multiplied with a window function w(t) such that The Fourier transform of such a function can be expressed as the Fourier transform of continuous function convolved with that of the window function, i.e., From the above relation, it is clear that the true underlying PSD of the observed light curve is in fact distorted by the effects of sampling pattern, i.e., window function. Besides, blazar light curves which are mostly dominated by power-law type noise can sometimes show a few cycles of periodic flux modulation, especially in the low-frequency domain see [9,[24][25][26], for the discussion.
Power response method PSRESP; [12], one of the methods extensively applied in the characterization of PSDs of AGN, incorporates these effects in its significance estimation method through the use of a large number of simulated light curves that mimic several properties of real observations, such as mean, standard deviation, sampling pattern and duration. Simulations of light curves can be performed following the Monte Carlo (MC) method described in [27]. First, the source periodogram is modeled with a power-law PSD of the form P(ν) ∝ ν −β + C; where ν and β represent temporal frequency and spectral index, respectively. For a light curve spanning the total duration of time T and with a mean flux ofx, Poisson noise C level can be given by where ∆x represent the flux uncertainties. In order to obtain a measure of significance of the periodic feature, χ 2 values are computed for the observed and the simulated light curves as following where P obs represents observed periodogram, and P sim,i , P sim (ν) and ∆P sim (ν) stand for periodogram of each simulated light curve, the mean periodogram and standard deviation of all the simulated periodograms, respectively. Then the ratio of the number of χ 2 i s greater than χ 2 obs to the total number of χ 2 i s in all simulations provides a measure of the probability than can be used to quantify the goodness of the fit for a given model. These ratios are computed for various trial shapes and slope indexes such that the best fitting PSD model has highest probability see [28,29], for further details. Here it is important to make a distinction between local and global significance. The local significance tells us how likely a periodic signal at a particular frequency is significant; whereas the global significance accounts for the fact that we do not have an a priori knowledge of the location of the frequency at which the signal might occur, and hence make considerations for all the periodogram frequencies see [9].

Results
Under a broader program of characterizing long term, multi frequency statistical properties of blazar sources, we performed time series analysis on a sample of blazars. Below we present particularly the periodicity analysis on four well-studies sources. In addition, the summary of the results of the analyses are presented in Table 1 that lists source names, red-shift of the sources, the waveband of the observations, total observation duration in approximate years, the number of observation points, the analysis methods applied to derive the periods, the source rest-frame periods approximately in days, and the number of cycles, approximated to nearest whole number, observed in the light curve, and the references for the previously published works in column 1-9, respectively.

•
Mrk 501: Decade-long Fermi/LAT observations (100 MeV-300 GeV) of the famous TeV blazar Mrk 501 were analyzed using four widely known methods: Lomb-Scargle periodogram, weighted wavelet z-transform, epoch folding, and z-transformed discrete auto-correlation function. The analyses revealed an underlying sub-year timescale (332 d) γ-ray QPO that persisted nearly 7 cycles. Applying PSRESP analysis the significance of the detection was found to be above 99% against possible spurious detection for details see [11]. • J1043+2408: A ∼ 560-d QPO was detected in the radio (15 GHz) observations of the BL Lac source J1043+2408. Multiple methods of time series analysis, such as epoch folding, Lomb-Scargle periodogram, and discrete auto-correlation function, were carried out on the source light curve. All three methods consistently found the repeating signals which were identified as being identical in periods within the uncertainties. In order to estimate the significance of the detection, a large number of Monte Carlo simulations were performed. This resulted in a significance of 99.9% 99.4% local and global, respectively, for details see [10]. • PKS 0219-164: The long-term variability properties of the blazar PKS 0219-164 were investigated utilizing the 15 GHz radio observations from the period 2008-2017. The methods Lomb-Scargle periodogram and weighted wavelet z-transform detected a strong signal of 270-day QPO along with possible harmonics of 550 and 1150 days periods. The statistical significance, both local and global, was estimated above 99% over underlying red-noise processes for details see [9]. • OJ 287: we performed a time series analysis of ∼ 9.2 year-long, densely-sampled optical light curve of the BL Lac OJ 287. Both the LSP and the WWZ methods provided evidence for the possible ∼ 400 d QPOs in the source light curve with a high significance of >99%. In addition, hint for a possible harmonic of ∼800 d period was also found for details see [4].

Discussion and Conclusions
We reported the presence of year timescale QPOs in four blazars in different wavebands by applying multiple methods of time-series analysis. It should be pointed out that the fact that all four sources in our studies are of BL Lac type blazars should not give the impression that QPOs are less frequent in FSRQs. In fact, the study of the possible differences between the types of QPOs detected in these two sub-classes of the blazars is yet to be carried out.
Time-domain analysis of multifrequency blazar light curves provides an excellent tool to investigate the physical processes occurring at the innermost regions of blazar, which otherwise are hardly accessible to current instruments. Such an approach can of fer important insights into a number of challenging issues including nature of space-time around fast-spinning supermassive black holes (SMBH), disk-jet connection and release of gravitational waves from the binary supermassive black hole systems. A number of possible scenarios can be linked to the observed year timescale for QPOs in blazars, some of which are discussed below.
First, from the observed period (P obs ), the corresponding period in the source rest frame (P) at the redshift of z is estimated using the relation P = P obs /(1 + z).
As a simplest case, the observed period might represent Keplerian period τ k around central black hole-such as a bright hotspot revolving around the black hole. In such case, the timescale can be used to estimate the radius of the orbit using the relation where a is the length of the semi-major axis of the elliptic orbits. Assuming circular orbits, for a typical black hole of mass of 10 9 M the radius of the Keplerian orbit for a year timescale can be estimated to be a few tens of gravitational radius (r g ). Similarly, the timescale might also represent the period of the supermassive binary black hole system (SMBBH) e.g., [7]. In such a system, elliptical orbits are gradually turned into circular orbits by dynamical friction that acts over the long course of merging. Therefore, the radius of the Keplerian orbit can be calculated using the relation, For a typical SMBBH system, the total mass of the system can fairly be assumed to be 1 × 10 9 M . Inferred from the year timescale periodicity, the black holes can be revolving around each other at a separation of a few milli-parsecs. Such a close binary system undergoes orbital decay due to the emission of gravitational waves. For the binary mass ratios in the rage 0.1-0.01 see e.g., [6,30], the orbital decay timescale in the GW-driven regime as estimated using the relation can be less than a thousand years see [31]. In such a case, we should expect that the gravitational coalescence of the system within a few centuries along with the emission of low frequency (∼10 −2 µHz) gravitational waves. But on the other hand, if the orbits in (SMBBH) are significantly elliptical, the secondary BH periodically might perturb the primary jet at periastron causing magnetohydrodynamic instabilities. These events inducing particle acceleration via magnetic reconnection events can release electromagnetic emission via synchrotron and inverse Compton emission see [32]. Moreover, the secondary black hole may perturb the disk of primary black hole such that the jets of the primary starts processing across the line of sight to result in the apparent QPO see [33].
In a different scenario, the observed year timescale QPOs could be associated with the instabilities taking place at the accretion disks. For example, oscillations developed in globally perturbed thick accretion can also give rise to QPOs and their fundamental frequency of the oscillation can be approximated as, see [34,35]. Similarly, QPOs can be linked to the precession of warped accretion disks around Kerr black holes residing at the hearts of blazars, known as the lense-thirring precession e.g., [36,37].
The timescale of such precession, τ LT , can be expressed in terms of the distance from the BH r and the mass of the central black hole M as given by where a s represent dimensionless spin parameter. For a maximally spinning (a s = 0.9) central black hole with a mass 10 9 M , a timescale of 360 days places the emission region around ∼12 r g . Furthermore, in the case where the accretion disk is strongly coupled to the jet e.g., in [38] (see also [25]), the precession of the disk can drive the jets precession. In fact, recent studies of 3D general relativistic MHD simulations (see [39]) based on the similar scenario has shown that the jets exhibit precession with a period of ∼1000 r g /c, which for a 10 9 M black hole turns out to be of the order of a few years. Moreover, the modulations arising near the central engine could propagate down the jet where the emission and the timescales are altered due to relativistic effects such that true periodic timescales at the BH could be longer than the observed according to the relation P = δ/(1 + z)P obs and δ = (Γ (1 − βcosθ)) −1 , where for Γ, β = v/c and θ are the bulk Lorentz factor, relativistic speed and angle to the line of sight. Equation (19) also suggests that QPOs in blazar could be related to the periodic swings in the angle between the emission regions and the line of sight. For example, emission regions moving along helical path of the magnetized jets can lead to the periodic changes in the viewing angle e.g., [40]. This way the Doppler boosted emission from the jet can be modulated periodically resulting in QPOs. In such cases, the rest frame flux (F ν ) is related to the observed flux (F ν ) through the equations F ν (ν) = δ(t) 3+α F ν (ν) and δ(t) = 1/Γ (1 − βcosθ(t)) .
Following the above argument, intrinsic flux remaining unchanged, the ratio of the observed flux modulation for a given change in the angle ∆θ can be obtained from the relation ∆logF = − (3 + α) δΓβsinθ∆θ.
As an illustration, for blazars having typical radio spectral index (α = 0.6) and viewing angles in the 1-5 • range, a slight change in the viewing angle e.g., ∼1.5 • , is sufficient to produce observed maxima twice as bright as minima in the light curves refer to Figure 4 in [10].
Finally, it is concluded that in blazars, we observe QPOs from mainly the jet emission that might be driven by perturbations near the black hole. In such a context, recurring gravitational perturbation in SMBBH systems driving the QPO production mechanism in jets serve as a most natural model. However, several issues linked with searching for QPOs in blazars such as red-noise dominance, modest significance of the detection and transient nature of the QPOs hinder the efforts to definitively characterize them. Therefore, future studies should be directed to more comprehensive multi-frequency time series analysis and collective discussion on the topic.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: