Search for Quasi-Periodic Oscillations in TESS light curves of bright Fermi Blazars

In a previous paper, we reported evidence for quasi-periodicities in the \textsl{TESS} light curves of BL Lacerate and two other blazars found serendipitously in the SDSS AGN catalog. In this work, we find tentative evidence for quasi-periodic features in the \textsl{TESS} observations of five sources in the fourth catalog of the Fermi--LAT (4FGL) sources: J090453.4$-$573503, J2345$-$1555, B0422+004, J002159.2$-$514028, and B0537$-$441. We analysed the \textsl{TESS} light curves of these blazars that we extracted using a customized approach. The quasi-periodic oscillations (QPOs) are searched for using two timing analysis techniques: generalized Lomb-Scargle periodogram and weighted wavelet Z-transform. Their apparent periods lie in the range of 2.8--6.5 days and have at least 3$\sigma$ significance in both of these methods. QPOs at such timescales can originate from the kink instability model which relates the quasi-periodic feature with the growth of kinks in the magnetized relativistic jets. We performed MCMC simulations to obtain the posterior distribution of parameters associated with this model and found the kink period consistent with previous studies.


INTRODUCTION
Supermassive black holes (SMBH) are believed to reside at the centers of galaxies and have masses between 10 5 -10 10 M ⊙ .The gravity of such a massive central compact object can lead to accretion from its surroundings so that the central region becomes highly luminous and sometimes even outshines the whole galaxy.Such a system is referred to as an active galactic nucleus (AGN) which is classified further based on different properties including optical spectroscopy, inclination towards the line of sight, and radio properties (e.g.Kellermann et al. 1989;Urry & Padovani 1995).Blazars are radio-loud AGNs having jets pointed toward the line of sight, thereby amplifying the observed brightness through relativistic effects.Blazars are further classified into BL Lacertae objects and flat spectrum radio quasars (FSRQs).They exhibit extreme variability throughout the whole electromagnetic spectrum, from radio to  rays (see Marscher et al. 2008, and references therein).
Quasi-periodic oscillations (QPOs) have been detected rarely in AGN but quite frequently in stellar-mass black hole and neutron star binary systems (Remillard & McClintock 2006).The timescales of reported QPOs in AGNs come from observations taken throughout the entire electromagnetic spectrum and range from a few tens of minutes to weeks and even years (see Halpern et al. 2003;Gierliński et al. 2008;Espaillat et al. 2008;King et al. 2013;Graham et al. 2015; ★ E-mail: ashutosht@tamu.eduSmith et al. 2018;Gupta et al. 2019;Tripathi et al. 2021;Jorstad et al. 2022;Tripathi et al. 2023, and references therein).
These QPO features could originate in the vicinity of the black hole and thus might be used to study the effect of gravity of the central object on its surroundings.Several models have been proposed to explain the origin of such features, including diskoseismology, warped accretion discs, and disc-jet coupling (Stella et al. 1999;Tagger & Pellat 1999;Wagoner et al. 2001;Kato 2005;Ingram et al. 2009).Recently, Čemeljić et al. (2022) performed three-dimensional general relativistic magnetic hydrodynamical simulations of accretion flows and argued that quasi-periodic episodic jets could arise from flux ropes formed via magnetic reconnection in the flows.The frequencies of QPOs associated with accretion discs scale inversely with the mass of the central object.Interestingly, this relation is valid for both stellar mass and many AGN black holes and suggests that the origin of these periodicities could be the same in all classes of black hole-powered sources (Abramowicz et al. 2004).
The Transiting Exoplanet Survey Satellite (TESS) (Ricker 2015) is a space-based optical instrument.Unlike temporal studies made with ground-based telescopes, TESS is not affected by seasonal gaps, irregular sampling, and degraded photometric precision arising from the atmosphere.The variability properties of a few blazars have already been studied using TESS observations (Weaver et al. 2020;Raiteri et al. 2021).Recently, we (Tripathi et al. 2023) claimed the presence of quasi-periodicity in TESS light curves of BL Lacertae itself and two other blazars discovered serendipitously in the SDSS-IV SPIDERS catalog, which essentially consists of X-ray selected AGNs from the ROSAT All-Sky Survey (RASS).
In section 2 we explain the sample selection and custom TESS data reduction.Data analysis methods are explained in section 3. The timing analysis results are discussed in section 4. We discuss some theoretical models that could explain the claimed QPO periods in section 5. Markov Chain Monte Carlo (MCMC) simulations of the kink instability model are discussed in Appendix A.

TESS Monitoring and Blazar Sample
The Transiting Exoplanet Survey Satellite (TESS; Ricker 2015) conducts long-term monitoring with high photometric precision of stars across nearly the entire sky with a rapid cadence: every 30 minutes in the early cycles (2018 -2019) and every 10 minutes or 2 minutes in later and current cycles (2020 -present).Except for a roughly 1-day long gap every orbit (approximately every two weeks), the coverage is continuous.The monitoring baseline for a TESS source depends upon its ecliptic latitude.The majority of sources at low-to-middle ecliptic latitudes are monitored for 27 days before the spacecraft shifts to the next wedge of the sky.Each 27-day monitoring segment is referred to as a "Sector."At higher ecliptic latitudes, these sectors overlap, resulting in a longer monitoring baseline, with a maximum of approximately 1 year of continuous monitoring at the ecliptic poles.After this year of exposure, the spacecraft turns over and surveys the other ecliptic hemisphere for 1 year with the same pattern.
To determine our sample for this analysis, we selected all objects from the fourth Fermi Large Area Telescope catalog of -ray sources (4FGL; Abdollahi et al. 2020) that are classified as BL Lac, FSRQ, or "unknown blazar type."We then cross-matched this list with the TESS Input Catalog (TIC; Stassun et al. 2019).The photometric precision of TESS is exquisite for bright sources, officially diminishing to 1% for sources with  ∼ 16 in pre-launch specifications (although actual performance is often better)1 .We, therefore, require that the TESS Magnitude (i.e., the magnitude in the TESS bandpass as reported by the TIC) be < 17.5.This cross-match and magnitude cut resulted in a total sample of 156 blazars.

Light Curve Extraction with Quaver
The light curve extraction and systematics correction were done using the Quaver software, specially designed with TESS AGN science in mind.In particular, the code allows the user to custom-select their extraction aperture from a cutout of the full FFI without downloading the entire image, which assists in avoiding nearby contaminating sources and can accommodate extended host galaxies.This is achieved using the TESSCut package (Brasseur et al. 2019).The program then makes use of several Lightkurve (Lightkurve Collaboration 2018) tasks.Standard principal component analysis of faint background pixels and background-corrected bright source pixels (excluding the target) is used to build design matrices housing the components of the additive background light and the more subtle multiplicative systematics affecting sources.These matrices are then used by Lightkurve's RegressionCorrector to correct the target's light curve.The approach is flexible, allowing the user to determine docs/tess/observing-technical.html the number of principal components used, and whether to handle the background using regression or simple subtraction.The Quaver code (Smith & Sartori 2023) is publicly available, with a detailed User Guide describing the different extraction modes2 .
Of the 156 blazars in our sample, 15 were excluded due to excessive crowding in the source region (usually at low Galactic latitudes), which precluded isolating the counterpart on TESS's large 21 ′′ pixels.Of the remaining sample, 44 sources are found to be variable from visual inspection.Of this subset, 12 sources are found to have indications of periodicity in their light curves supported by one or more different timing analysis techniques.Finally, only 5 sources (≈ 3%) appear to have QPOs that survived our significance tests.Four of these light curves show at least 3  significance in both of the analysis approaches we used and the other one has such high significance in one of those methods.

Generalized Lomb Scargle Periodogram (GLSP)
The periodogram of a given time series can be estimated by a discreet Fourier transform and used to search for periodicities.The Lomb-Scargle periodogram (LSP) is a commonly used Fourier transform method employed in analyzing irregular and unevenly sampled observations (Lomb 1976;Scargle 1982) by fitting a time series with sinusoidal components.The generalized LSP (GLSP) is a modified version of the LSP where not only sinusoidal and constant components are used to fit a given time series but errors in the fluxes are also included in the calculations (Zechmeister & Kürster 2009).In this work, we use the GLSP code provided by pyastronomy3 .The employed code uses Leahy normalization which includes the number of trials in calculations so that the significance is intrinsically globally estimated.Here, we set the oversampling factor to be 2.0 which means that the number of frequencies is equal to the number of data points in the observations.We also tried oversampling factors of 1.0 and 10.0 to assess their effects on the significance level.The significances calculated using these different oversampling factors are consistent within 1-2% for all light curves analysed in this work.In this work, the analysed light curves are too short to make a clear determination of the underlying baseline model for the power spectrum.So we have employed models which assume that the light curves of AGNs typically follow red noise due to stochastic processes in accretion discs or jets.
We fit the power spectrum of a given data both with a simple power law and with a bending power law using the scipy python routine.We then estimate the uncertainties on the parameters through MCMC simulations using 500 walkers and 200,000 iterations.The simple power law used is where  is the normalization,  is frequency,  0 is a constant, and  is the power-law index.The bending power law is expressed as where   is the bending frequency,  is the power-law index, and  is the constant representing the instrumental white noise that dominates at the highest frequencies.
We then simulate 10,000 light curves, with the best-fit parameters and statistical properties of the observation obtained from the timing analysis package stingray 4 (Huppenkothen et al. 2019) which uses an algorithm based on Timmer & Koenig (1995).Each sector of TESS light curves have a gap in their middle due to the orbital motion of the instrument and have a duration in the range of ∼1-4 days.So, we 4 https://github.com/StingraySoftwaresimulate each light curve with a such a gap whose length is similar to that found in the real TESS observations.For QSO B0537−441, two consecutive sectors are analysed.In this case, we also include the sector gap in our calculation in addition to the gaps present in the middle of each sector's observations.We calculate GLSPs for all these simulated light curves and obtain the distribution of power spectra at individual frequencies.The spectral distribution at each frequency can be used to estimate the significance of a peak for a given confidence interval.For instance, the 99.73 % confidence interval is estimated by the 99.73 percentile of the power distribution of these 10000 simulated light curves.Weighted Wavelet Z-transform (WWZ) A periodic signal may not persist throughout an observation.Wavelet analyses allow for and assess the possible non-stationarity present in a given time series by decomposing the signal in the time and frequency domains simultaneously.WWZ is a commonly used wavelet technique used to analyse unevenly sampled and sparse time series (Witt & Schumann 2005).
We employed a WWZ code5 used commonly for such analysis (see Zhang et al. 2017;Tripathi et al. 2021Tripathi et al. , 2023, and references therein).
When we average the WWZ power over the whole period of duration, we get the time-averaged WWZ, which is a periodogram following a  2 distribution with two degrees of freedom (Foster 1996).As this time-averaged WWZ is a Fourier transform, we follow an ap- Periodogram (GLSP) and Weighted Wavelet Z-transform (WWZ) analysis methods.
Fig. 1(a) shows the timing analysis results obtained for the 25-day sector 37 (cycle 3) observation of ATPMN J090453.4−573503.A strong signal at a period of around 2.8 days is at least 3 significant in the GLSP and WWZ results regardless of which of the red noise models is employed.The signal is found to be persistent throughout most of the observation as shown in the colour density (violet indicating the least power and red the most) plot of WWZ power (left plot of bottom panel) which also shows the presence of non-stationarity in the observation.The time-averaged WWZ (right plot in bottom panel) supports the presence of the peak at the claimed frequency with 3 significance.There is also a peak around 4.5 days which is 3 only for the simple power law.This signal is also found in the WWZ analysis but has less power as compared to the peak around 2.8 days.
Fig. 1(b) shows the results obtained for the 25-day sector 29 (cycle 3) observation of QSO J2345−1555.In the GLSP result, a strong peak having significance at least 3 is detected at 3.5 days assuming a either simple or a bending power law as the underlying model.The wavelet analysis shows this peak is entangled with other lower significance peaks and together, they are much stronger in the first half of the observation.In the time-averaged WWZ plot, this peak at 3.5 days has a significance of 3, so we accept it as meeting our full set of criteria or a tentative claim of a QPO.
The light curves and the timing analysis results of QSO B0422+004 are plotted in Fig. 2.This source was observed during two different TESS cycles.The results for the 25-day sector 5 (cycle 1) observation are shown in Fig. 2(a).Using a bending power law, the peaks around 3.1, 4.5, and 6.0 days are found to be at least 3 significant, whereas with the simple power law, only the peak around 3.2 days is found to be significant.So this peak is at least 3 significant with both underlying models.The WWZ power of the 3.1-day peak is somewhat less than that of the peaks around 4.1 and 6.0 days.It persists throughout the observation with similar WWZ power; however, this is not the case for the 4-day signal which is strong only in the first half of the observation.In the time-averaged WWZ result, the peak around 6.0 days has more than 3.However since the length of the 25-day segment contains only 4 cycles of a putative 6-day QPO, it is difficult to make a strong claim for it based solely on this sector.The 3.1-day peak is found to be at least 3 in both GLSP and WWZ methods and also corresponds to 8 cycles of a putative QPO.The analysis of the other 25-day sector 32 (cycle 3) data for this source (Fig. 2(b)), shows a strong (3) signal at 6.0 days in the GLSP with both red noise models.This is also present in time-averaged WWZ analysis with 3 significance, respectively.With a simple power law, a peak around 3 days is also found to be more than 3 significant.The WWZ colour diagram indicates the persistence of this 6.3-day signal throughout the observation.The fact that both the widely separated cycle 1 and cycle 3 observations indicate a QPO of around 6 days strengthens our confidence in its presence.The ∼6-day signal could be harmonically related to the signal detected at ∼3-day but because of the lower significance of the latter using both GLSP and WWZ methods, we take the former to be the most likely QPO period.
In Fig. 3(a), our timing analyses for the 27-day sector 2 observation of 1RXS J002159.2−514028 are shown.A peak of around 4.5 days is found which is at least 3 in both the GLSP and time-averaged WWZ results using both simple and bending power laws as underlying red noise models.This 4.5-day signal is found to be stronger in the first half of the observation, as displayed in the WWZ colour diagram so the claim of a QPO for this source is less strong.From the timeaveraged WWZ plot one could conclude that the peak at around 11 days is more significant than that at 4.5 days.However, the period of 11 days corresponds to only 2 cycles for the 25-day observation which is not enough to be considered as quasi-periodic.We restrict our analysis in the frequency range of 0.1-2.0day −1 because the peaks with frequency less than 0.1 day −1 corresponds to a period of more than 10 days which would not be enough for a QPO to be convincing.Fig. 3(b) shows the results of the combination of consecutive observations of QSO B0537−441 during sectors 32 and 33.So, unlike other TESS light curves, this one has a duration of 53 days.In the GLSP plot, the signal at 6.5 days is found to be at least 3 significant using both red noise models which is also confirmed with time-averaged WWZ analysis.The colour density WWZ plot shows long-lived strong concentrations of power around 6.5 and 11 days.The tentative QPO of around 6.5 days is supported by both methods with at least 3 significance.

Effects of gaps on periodogram
The light curves used in this work contain a gap between the sectors due to the orbital motion of the instrument during which the instrument points toward Earth to deliver data.The duration of this gap varies from 1 to 5 days depending on the pointing anomalies associated with the observation and on the cadence masking during the light curve extraction from quaver which essentially removes the noise from the data.These gaps are taken into account in the calculation of the significance of the peaks in the GLSP and WWZ analyses.Here we will quantitatively estimate the effect of these gaps in the error estimate of the periodogram.We first simulate 10,000 26-day light curves with no gaps with the methods explained in the previous section and then estimate the periodogram and the associated 1 confidence interval.In these simulations, we assume both simple power law and bending power law as underlying red-noise models.We repeat these simulations for the light curves with a gap of 2, 4, and 6 days in the middle of the observation.Fig. 4 shows the 1- error associated with the periodogram for different durations of gaps in the light curves that are simulated assuming both a simple power law and a bending power law.For both underlying red noise models, the error estimates increase as the duration of the gap increases.It implies that the significance of the peaks calculated for the light curves in this work may be overestimated but the significance likely would have been greater if the gaps were absent or at least shorter, than those present in these light curves.

DISCUSSION AND CONCLUSIONS
The majority of the observed emission from blazars is contributed by non-thermal processes in jets and typically substantially exceeds the total emission from the accretion disc and host galaxy (e.g.Jorstad et al. 2001).This is due to the special relativistic effect of Doppler boosting produced by the close alignment of the jet to the line of sight of the observer (Urry & Padovani 1995).These dominant relativistic enhancements in jets manifest themselves in the light curve of blazars as strong variability occurring on shortened timescales and can be observed throughout the electromagnetic spectrum, from radio to rays; these variability timescales range from days (e.g.Ackermann et al. 2016;Ahnen et al. 2017) or weeks to even years (e.g.Villforth et al. 2010;Acciari et al. 2014).Although the temporal power spectra of these blazar variations are dominated by red noise, a small fraction of them have QPOs that are statistically significant (see Hayashida et al. 1998;King et al. 2013;Graham et al. 2015;Jorstad et al. 2022;Tripathi et al. 2023, and references therein).
Several processes can explain the apparent QPOs with periods of the order of a few days such as we have tentatively found here.Such periodicity could be related to transient oscillations occurring in various localized regions of a turbulent accretion disc (Abramowicz et al. 2004) that are advected into the jet.Perez et al. (1997) shows full relativistic calculations for adiabatic oscillations confined in inner regions of the accretion disc due to non-Newtonian gravitational effects of the black hole.QPOs can also be produced naturally by jet precession that produces periodic changes in Doppler boosting (e.g.Gopal-Krishna & Wiita 1992).Such jet precession could arise from by Lens-Thirring precession of the accretion disc (Stella & Vietri 1998).However, these jet precession models are better fits for QPOs with periods of several years.Since the periods of QPOs reported in this work are of the order of a few days, any jet precession model can probably be excluded.
One very attractive model for such quasi-periodic features that is completely internal to the jet involves the development of a kink instability (e.g.Dong et al. 2020).The twisting of the toroidal magnetic field and resulting particle acceleration leads to the formation of a kink which is quasi-periodic in distance along the jet.This periodic kink evolves with time in the region of compressed plasma in the magnetized jet and could be detected as a quasi-periodic signature with a period of a few days to weeks (see Barniol Duran et al. 2017;Jorstad et al. 2022, and references therein).The transverse motion of the kinked region can be used to estimate the growth time of the kink which can drive such quasi-periodic signatures.The kink growth time (   ) can be calculated as the ratio of transverse displacement of the jet from its center (   ) and average velocity of motion (⟨  ⟩)

Figure 1 .
Figure 1.Results for (a) ATPMN J090453.4−573503 during TESS sector 37 and (b) QSO J2345−1555 during sector 29.The upper panels show the light curves analysed in this work for these two sources, where blue dots show the actual data points, gray denotes the associated error, and the red curve corresponds to the running average.The middle panel for each source shows a generalized Lomb-Scargle periodogram (GLSP).The left plot of the bottom panel for each source shows the weighted wavelet Z-transform analysis and the right plot shows the time-averaged WWZ power.The red dotted-dashed curves in the GLSP and time-averaged WWZ plots denote the claimed QPO signal.Blue dashed curves represent the 99.73% (3 ) global significances, while the green dashed curve in the GLSP results corresponds to the power spectrum obtained from a bending power-law fit.

Figure 4 .
Figure 4. 1 error bounds for the periodogram estimated for different durations of gaps in the middle of observation assuming a bending power law (left) and a simple power law (right) as the underlying red noise model for simulating light curves.

Figure A. 1 .
Figure A.1.Corner plots for the posterior distribution of parameters, related to the kink instability model, sampled using 100,000 iterations and 100 walkers through MCMC simulations.Each panel is labelled with the AGN name (and observation Sector numbers for B0422+004).MNRAS 000, 1-?? (2023)