Optical Quasi-periodic Oscillations in the TESS light curves of three blazars

We report the time series analysis of TESS light curves of three blazars, BL Lacertae, 1RXS J111741.0+254858, and 1RXS J004519.6+212735, obtained using a customized approach for extracting AGN light curves. We find tentative evidence for quasi-periodic oscillations (QPOs) in these light curves that range from 2 to 6 days. Two methods of analysis are used for assessing their significance: generalized Lomb-Scargle periodograms and weighted wavelet Z-transforms. The different approaches of these methods together ensure a robust measurement of the significance of the claimed periodicities. We can attribute the apparent QPOs to the kink instability model which postulates that the observed QPOs are related to the temporal growth of kinks in the magnetized relativistic jet. We confirm the application of this model to BL Lacertae and extend the kink instability model to the other two BL Lac objects.


INTRODUCTION
Active Galactic Nuclei (AGNs) are among the brightest non-transient objects in the Universe and are understood to possess supermassive black holes (SMBHs) at their centers.Objects of such high luminosity (∼ 10 47 erg s −1 ) are thought to be powered by the accretion from their surroundings onto the SMBHs.AGNs are often classified based on their optical spectroscopic properties, which depend on their inclination to our line of sight, and on the basis of their radio emission, as radio-loud and radio-quiet objects (e.g.Kellermann et al. 1989).Radio-quiet objects are believed to have either weak radio jets emerging close to the central object and not propagating very far or a lack of radio jet emission.By contrast, blazars are radio-loud AGNs with low inclinations, and thus with jets aligned close to the line of sight which make relativistic effects significant (e.g.Urry & Padovani 1995).Blazars are usually considered to be the union of BL Lacertae (BL Lac) objects, which have very weak optical and UV emission lines, and flat spectrum radio quasars (FSRQs), for which the broad emission lines that are a defining characteristic of all quasars can be seen.Blazars have a rapidly varying continuum spanning the radio to -ray wavebands and exhibit strong polarization.The Doppler-boosted emission enhances the variability throughout the electromagnetic spectrum and blazars can be distinguished by these and other extreme properties (Scheuer & Readhead 1979;Blandford & Königl 1979).Much of the radio and optical variability is thought to be due to instabilities that develop in the jet (e.g.Marscher & Gear 1985;Marscher et al. 2008).
★ E-mail: ashutosht@tamu.eduQuasi-periodic Oscillations (QPOs) have been observed in X-ray light curves of stellar-mass black hole and neutron star binaries (see Remillard & McClintock 2006, for a review).The QPOs discovered in these systems are classified into low-and high-frequency QPOs.Low-frequency QPOs (LFQPOs) are very common and happen at frequencies ranging from mHz to 30 Hz.The High-frequency QPOs (HFQPOs) are much rarer and are only found in certain states of hardness and luminosity, indicating that thermal disk emission is the dominant component.The associated frequencies often occur in a 3:2 ratio (e.g.Remillard & McClintock 2006;Motta 2016).In neutron star and stellar-mass black hole systems, these features are believed to originate from the vicinity of the black hole, probably within a few gravitational radii.Thus QPOs in these sources can be used to study strong gravity and its effects on the surrounding plasma.Various models have been proposed to explain the origin of LFQPOs (Motta 2016), including Lens-Thirring precession (Ingram et al. 2009) and unstable spiral density waves (Tagger & Pellat 1999).For HFQPOs, other theoretical models could explain the periodicities; these include warped accretion disks (Kato 2005), disk-jet coupling (Li & Narayan 2004), diskosiesmology (Wagoner et al. 2001), and the relativistic precession model (Stella et al. 1999).Despite all these efforts, the physical mechanisms for the origins of QPO are still debatable.Although rare, QPOs have also been apparently observed in systems possessing SMBHs (Halpern et al. 2003;Gierliński et al. 2008;Lin et al. 2013;Zhang et al. 2017;Bhatta 2017;Smith et al. 2018;Gupta et al. 2018Gupta et al. , 2019;;Tripathi et al. 2021, and references therein).The frequencies of the HFQPOs are roughly 40 times greater than those of the LFQPOs, in both stellar mass and SMBH systems (Zhou et al. 2015;Smith et al. 2021).
QPOs have been reported in AGNs across the whole electromagnetic spectrum, from radio to -rays, and claimed timescales vary from a few tens of minutes to days, to weeks, and even decades.Primarily based on the timescales of periodicity, different models have been proposed to explain these QPOs.Valtaoja et al. (1985) reported a periodicity of 15.7 min in the radio light curve of OJ 287 taken in 1981.Carrasco et al. (1985) claimed the periodicity of 27 min in the optical observation of the same object taken 2 years later which they suggested could be explained by the presence of "hotspots" in a thick accretion disk.Quasi-periods of the order of a few days have been reported for several AGNs (e.g.Halpern et al. 2003;Smith et al. 2018;Jorstad et al. 2022;Kishore et al. 2023).Models such as diskoseismic modes, Lens-Thirring precession, and kink instabilities were employed in those papers to explain these day-like QPOs.However, the nature of such periodicities remains uncertain and thus it is important to search for more such fast putative QPOs.
On longer timescales, Sillanpää et al. (1996) reported the 11.7year periodicity in the century-long optical light curve of OJ 287.Periods on the order of a few years have also been reported in several other (e.g.Raiteri et al. 2001Raiteri et al. , 2003;;Foschini et al. 2006;Villata et al. 2009;Wang et al. 2014;Graham et al. 2015;Sandrinelli et al. 2017;Roy et al. 2022;Li et al. 2023).Those papers suggested various mechanisms that might explain such year-long periodicities; these include the presence of a binary supermassive black hole system, a plasma blob moving helically inside the jet, and precession of the jet.
The frequencies associated with the accretion disk-based oscillations are inversely proportional to the mass of the central object.This remarkable relationship has been observed in both stellar mass and SMBHs with masses spanning many orders of magnitudes from a few solar masses to billions of solar masses (e.g.Abramowicz et al. 2004;Zhou et al. 2015).This relation suggests that the origin of such periodicity may be the same for both kinds of black holes.Moreover, this relation can be used to deduce the mass of an object if periodicity is detected in its light curve.Although SMBH masses can be plausibly estimated for some FSRQs (e.g.Foschini et al. 2015;Xiao et al. 2022), for BL Lacs, there are very few reasonable measurements of SMBH mass and they are model dependent (e.g.Woo & Urry 2002;Nguyen et al. 2021).
Optical QPOs are rarely detected in ground-based surveys due to their uneven and irregular sampling.If periodicity is searched for in such data, there is a serious risk that the stochastic red noise characteristic of AGNs can either mimic or mask a quasi-periodic signal (Vaughan et al. 2016).The Transiting Exoplanet Survey Satellite (TESS, Ricker 2015), similar to its predecessor Kepler (Borucki et al. 2010), has substantially overcome these issues with improvement in precision over that yielded by the modest aperture ground-based optical telescopes normally used in these studies and with high-cadence regular sampling they cannot provide.Thanks to TESS, space-based regularly sampled light curves with high precision and high cadence can be produced and used for studying periodicity at higher confidence than ground-based surveys.In this work, we search for QPOs in the high-cadence and regularly sampled TESS observations and see that we can interpret the periodicities seen in three blazars in terms of the kink instability model.
Our team has developed a specialized approach to extracting TESS light curves of stochastically-varying sources like AGN.Smith et al. (2018) detected a 44-day low-frequency QPO in a Kepler light curve using a preliminary version of our current approach, which has been significantly improved to avoid over-fitting of intrinsic AGN variability while handling both background light and instrumental systematics in a careful and flexible way.Weaver et al. (2020) analyzed the multi-wavelength observations including the observation from TESS of BL Lacerate (BL Lac afterward) and suggested a 13-hour timescale.A few other papers also have examined the variability of individual blazars using TESS observations (Raiteri et al. 2021a,b;Kishore et al. 2023).
In this paper, we report candidate QPO detections in the TESS light curves of BL Lacertae itself, as well as in those of two other BL Lac objects: 1RXS J004519.6+212735 and 1RXS J111741.0+254858(which has two non-contiguous TESS observations).In Section 2, we describe the sample selection and our method of reducing data taken from TESS.In Section 3, we briefly explain the various methods used for QPO detection.Results of those analysis are reported in Section 4 and conclusions are given in Section 5.

Sample Selection
The candidate periodicities presented in this work were discovered as part of an ongoing large-scale investigation of the TESS light curves of Seyfert galaxies.We are extracting TESS light curves of 5971 Type 1 Seyferts from the SDSS-IV SPIDERS catalog: a database of X-ray selected AGN from the ROSAT All-Sky Survey and the XMM-Newton Slew Survey with either eBOSS or SDSS DR12 spectral coverage (Comparat et al. 2020).Of these objects, the preliminary periodogram analyses indicated 22 had possible periodicities.Of these, three objects survived the more rigorous test of periodic significance described below in Section 3. Two of them, described in this work, were BL Lac objects, and one is a Seyfert, which will be discussed in terms of the broader Seyfert sample in an upcoming work.
Because the BL Lac objects were contaminants in the SPIDERS sample of X-ray selected Seyferts, it is not possible to use these numbers to draw statistical conclusions about the occurrence rates of periodicities in blazars.In this large sample, there were only 28 blazars that exhibited detectable variations.The two presented here, alongside the TESS light curve of BL Lacertae itself, are the only ones showing possible indications of QPOs and are serendipitous discoveries that can provide support for the physical interpretation of the QPO in BL Lacertae reported recently by Jorstad et al. (2022).

The Transiting Exoplanet Survey Satellite (TESS)
The TESS instrument, designed primarily to search for transiting exoplanets using long-term monitoring at high photometric precision of stars across nearly the entire sky (Ricker 2015), provides very high cadence optical monitoring, with a cadence of 30 minutes in the early cycles (2018 -2019) and 10 minutes or 2 minutes in later and current cycles (2020 -present).Coverage is nearly continuous, as there are no gaps due to daylight or seasonal sunlight.TESS light curves are therefore significantly better-sampled than ground-based optical light curves, with longer continuous baselines.The length of nearly continuous monitoring for a TESS source depends upon its ecliptic latitude, ranging from 27 days nearest the ecliptic plane to approximately 1 year at the ecliptic poles.After each year-long Cycle, the spacecraft swaps hemispheres and performs the same survey in the other half of the sky.Thus, sources will occasionally have two "Sectors" of data, separated by one year.The mission began in Cycle 1, in the southern hemisphere, then went to the northern sky in Cycle 2, and returned to the south for Cycle 3, then the north for Cycle 4. The photometric precision of TESS is dependent upon the source magnitude; at the magnitudes of our sources, variability at the ∼ 1 − 10% level can be confidently detected.

TESS Data Reduction
Because the TESS spacecraft was designed to detect a periodic signal in a point source, special treatment is required to derive light curves of stochastic signals that may be in resolved sources like bright host galaxies.Many of these techniques to study AGN were first developed using very similar data from the Kepler spacecraft (Smith et al. 2018).We have developed these techniques significantly, to avoid over-fitting of intrinsic AGN variability by carefully modeling the background additive light curve effects and the subtler multiplicative systematics in separate phases using principal component analysis (PCA), then removing these effects using linear regression1 (Smith & Satori 2023).
We first download a "postage stamp" of the TESS full-frame images around our sources, 25×25 pixels on a side.From this image, we select a custom extraction aperture for each target, allowing us to avoid the inclusion of nearby sources and to build an aperture appropriate for any spatially resolved host galaxies.
Next, we define a flux threshold of 1.5  above the background to determine what pixels in the downloaded postage stamp are likely to contain only background sky or very, very faint sources at the background level.These faint pixels are fit by three principal components.The results of the fit are stored in a design matrix, representing the additive systematics dominant in the light curve.This method makes use of the matrix regression methods in the publicly-available TESSCut (Brasseur et al. 2019) and lightkurve (Lightkurve Collaboration 2018) software packages.
Next, all of the pixels in the field excluded from the previous mask, i.e., the "bright" pixels, but excluding the aperture of the source itself, are corrected by these additive background systematics.Once corrected, the variability of these bright pixels is again fit by three principal components, which are stored in a design matrix; this matrix describes the higher-order and more subtle multiplicative systematics.Armed with these two matrices, we combine them into a hybrid design matrix that contains principal component vectors for additive and multiplicative effects.
Finally, the raw light curve of the source is extracted from its aperture and linearly regressed against the hybrid design matrix.This method was inspired by the various corrector techniques available from lightkurve and makes heavy use of its RegressionCorrector class2 .
Figure 1 shows the light curves of the sources analyzed in this work.The upper left panel shows the light curve of BL Lacertae.TESS observed this source beginning on 12 September 2019 for around 24.6 days in sector (Sec.)16 of Cycle 2 observations.The gray curves represent the error associated with the fluxes which are shown in blue.The red curve is the running average of the light curve taken with the window=20 to highlight the possibly periodic features present in the data (window defines the number of measurements for which the running average is calculated).The upper right panel shows the light curve of 1RXS J004519.6+212735.This BL Lac object was observed beginning on 8 October 2019 for about 23.7 days in Sec. 17 of Cycle 2 observations.Both lower panels belong to 1RXS J111741.0+254858.The Sec. 22 observation of Cycle 2 began on 21 February 2021 and continued for about 25.6 days; it is shown on the lower left panel of Fig. 1.The lower right panel shows the light curve of the Sec.49 observation of Cycle 4 of 1RXS J111741.0+254858 starting on 28 February 2022 and continuing for about 24.8 days.Figure 2 shows the light curve of two inactive stars from the same postage stamps as 1RXS J004519.6+212735 and 1RXS J111741.0+254858,reduced using the same pipeline.The reduced light curves are flat, indicating our successful removal of systematics.

DATA ANALYSIS
In this section, we discuss the timing analysis techniques used in this work to search for any periodicities or quasi-periodicities and to assess their significance in the four light curves shown in Fig. 1.These methods are the generalized Lomb-Scargle periodogram and weighted wavelet Z-transform analyses.

Generalized Lomb-Scargle Periodogram
The Lomb-Scargle periodogram (LSP) is a powerful Fourier transform frequently used to detect periodicity in unevenly sampled data (Lomb 1976;Scargle 1982).It uses likelihood minimization to fit sinusoidal components throughout the data.The Generalized LSP (GLSP) is an advanced version of the 'classical' LSP which includes the error of fluxes in the periodogram calculations and fits the data with both sinusoidal and constant components.
In our work, we employ the GLSP routine provided by pyastronomy3 .It applies the method described in Zechmeister & Kürster (2009).In the code, the frequency resolution is governed by the oversampling factor which is needed for the calculation of the power spectrum and its significance.Here, we took the oversampling factor to be 2.0, which corresponds to the number of frequencies considered to be equal to the number of data points in a single observation.We have confirmed that the result is insensitive to different oversampling factors, by also employing oversampling factors of 1.0 and 10.0.We found that the values of the frequencies of the peaks of the signals exceeding 3 significance vary by no more than 3% for that range of oversampling factors and the value of the powers at those peaks are within just a few percent of each other.

Weighted wavelet Z-transform analysis
In addition to detecting any possible periodicity in a signal, it is also imperative to assess its persistence throughout the observation.The wavelet method decomposes the signal into both frequency and time domains simultaneously.Instead of the sinusoidal components, this method uses wavelet functions to fit the observation.In this work, we employed the weighted wavelet Z-transform (WWZ) method which can be used effectively for unevenly sampled and sparse data (Witt & Schumann 2005).It uses the z-statistics for the calculations (Foster 1996) and is based on the Morlet wavelet function which is the function of both frequency and time.In this work, we use the publicly available WWZ4 software that has been used extensively in the literature (see Bhatta 2017;Zhang et al. 2017;Tripathi et al. 2021, and references therein).
The major assumption in wavelet calculations is the type of wavelet function used.The finite length of the observation can induce the edge effects which could be observed in the density plots.We also note that the frequency and temporal resolutions can have important effects on the calculation of the significance of peaks.

Significance Calculations
It appears that the stochastic processes occurring in both turbulent jets and accretion disks usually manifest themselves in the light curves of blazars as is seen in observations taken throughout the electromagnetic spectrum (e.g.Max-Moerbeck et al. 2014;Nilsson et al. 2018;Bhatta & Dhital 2020).A blazar Power Spectrum Density (PSD), (  ), as a function of frequency  , often follows a simple power-law of the form   where  is the spectral index for a significant range of frequencies.Any QPO, on the other hand, is believed to arise from coherent disk or jet processes.So, a strong QPO feature should project itself in the blazar PSD above the power-law.It is imperative to estimate the significance of the peak relative to the red noise of the PSD.As  usually has a negative value, meaning that the spectrum at lower frequencies fluctuates with larger amplitudes, it is possible to misinterpret these fluctuations as quasi-periodic signals.In addition, there may be some fluctuations introduced in the signal by processes arising from the observational constraints which could be misinterpreted as periodic peaks.Therefore, it is very important to estimate the significance of these peaks.
To perform significance calculations, we first simulate light curves having the same properties as the observation (mean, variance, duration) and power spectrum.Using the standard discreet Fourier Transform (DFT) for the unevenly sampled data as used in this work would be problematic as it results in amplitude fluctuation and unreasonable frequency shifting (e.g.Fan et al. 2017).Although the TESS data used in this work are far more uniformly spaced than most observations, they do contain some gaps, so the observations are not taken at strictly regular intervals.To avoid these issues, we employ the GLSP method.The light curve analyzed in this work is essentially two segments separated by a rather wide ∼4-day gap (except for BL Lac, where the gap is only ∼1 day).To simulate these light curves realistically, we essentially treat segments individually.The resultant simulated light curves consisted of the two ∼10-day segments with ∼4-day gap in the middle.We chose the frequency spacing such that the number of frequencies is equal to the number of data points in the observation.We fit the power spectrum () at a given frequency  with a bending power-law plus a constant representing the instrumental white noise that dominates at high frequencies where  is the normalization,   is the bending frequency,  is the power-spectrum index, and  is the constant representing the instrumental white noise.We used the curve fit routine in scipy package of python to estimate the parameters of the red-noise spectrum.We also fit the power spectrum with simple power law defined as where  0 is the constant.Although a simple power-law and bending power-law are chosen to model the 'red-noise' in the spectra, we note that it is possible to explain the stochastic process using different mathematical models; e.g., short memory ARMA models are capable of producing periodic patterns for AR(p) with p>2 (see Schulz & Mudelsee 2002, and references therein).
Using the parameters of each power spectrum obtained using GLSP and the other properties of the source, we simulate 10000 light curves employing the stingray (Huppenkothen et al. 2019) timing analysis package5 which implements the algorithm described in Timmer & Koenig (1995).We calculate GLSP for the simulated light curves.The significance is estimated from the PSD distribution at each frequency.
We also perform Monte Carlo Markov Chain (MCMC) simulations to estimate errors on the best-fit parameters of the best-fit bending power law.Table 1 and Table 2 show the best-fit parameters and the associated errors for the bending power-law model and the simple power law, respectively.
To estimate the significance of the signal in the WWZ method, we calculate the time-averaged WWZ, which is essentially the WWZ power marginalized over time.This time-averaged WWZ is a periodogram that follows the  2 distribution with 2 degrees of freedom in the limit of even sampling (Foster 1996).We calculate time-averaged WWZ for 10,000 simulated light curves.The significance, as in the case of GLSP, is derived from the distribution of the simulated timeaveraged WWZ at each frequency.

RESULTS
Fig. 1 shows the light curves we obtained from the TESS observations for BL Lac (upper left panel), 1RXS J004519.6+212735(upper right panel), and 1RXS J111741.0+254858(lower panels).Visually, there seem to be possible quasi-periodic components that were searched for and confirmed by the GLSP and WWZ methods.These methods were chosen as they provide distinctive approaches to the time series data.Explicitly, GLSP gives the Fourier transform of the data, probing the frequency domain while the wavelet analysis decomposes the data into frequency and time simultaneously.Only if a peak shows at least 3  significance in both the GLSP and time-averaged WWZ methods do we consider it to be a solid QPO signal.Fig. 3 shows the normalized power (black curve) of the GLSP plotted against frequency for BL Lac.The blue curves denote the 3 (99.73%) significance obtained by simulations assuming simple power law as the underlying model.The green curve corresponds to 3 significance obtained by simulations assuming bending power law.The most prominent peak, with a significance of more than 3 , occurs at the period of 4.94 days.The peak around 1.40 days is also prominent (> 3  significance) for both bending and simple power law.
Note that quoting a period of more than 6 days would correspond to having fewer than 4 cycles for the timescales of our observations and therefore are not suitable for claiming a QPO.This limitation is supported by early work which showed that it is possible to apparently see 3 cycles in many light curves based purely on flicker noise (Press 1978).Vaughan et al. (2016) found that greater numbers of cycles that follow a periodic pattern more strongly support the presence of periodicity.Even a pure sinusoidal component could not be distinguished from a simple stochastic process if there are only ∼ 2 cycles.However, the presence of ∼ 5 cycles could be easily distinguished from a stochastic process.
As the time-averaged WWZ is also effectively a periodogram that contains similar information as the GLSP, showing only one underlying model in the WWZ plots is sufficient.We note that the claimed frequencies in the Fourier space lie in the "red-noise" regime where the PSD can be adequately described by both simple power-law and bending power-law models.So we show the WWZ significance plots only for bending power laws.
The upper left panel of Fig. 4 shows the result of our WWZ analysis for BL Lac.The left plot shows the color-color diagram of WWZ power plotted against time and frequency.The maximum power is concentrated around the period of 5.26 days, and it persists throughout the observation.The second most prominent concentration, which also persists throughout the observation, is around 2.94 days.There are other weaker features that seem to be present in different parts of the observation.For instance, the feature around 1.54 days becomes stronger in the later half of the observation whereas the feature around 2.0 days is only visible in the first half.The black curve in the right sub-panel of Fig. 5 for BL Lac shows the time-averaged WWZ and the red dotted curve denotes the 3 (99.73%) confidence curve.The feature at 5.26 days, which also appears in the WWZ plot, is found to be above 3.The signal around 5 days is found to be significant in both methods.However, as it corresponds to fewer than 5 cycles during the 24.6-dayTESS observations, it cannot be considered to be very convincing.
The upper right plot of Fig. 3 shows the GLSP for Cycle 2 data of 1RXS J004519.6+212735.The peaks at the period of 5.8 days are marginally consistent with 3  significance.The color-color diagram of the WWZ analysis of this observation (upper right panel of Fig. 4) shows that the peak around 5.88 days is persistent throughout the observation and is the only peak having significance more than 3 in the time-averaged WWZ plot on the right-hand side of the panel.Although the peak at 5.9 days is nominally the most significant measurement from both methods, it corresponds to only 4.0 cycles during the length of the observations, so cannot be fully trusted.
We now discuss the results obtained for 1RXS J111741.0+254858.There are two observations of this source, one from sector 22 of Cycle 2 and one from sector 49 of Cycle 4. The lower left panel of Fig. 3 shows the GLSP for Cycle 2 data.The most prominent period (at less than 6 days) is found at 3.68 days.The feature around 3.68 days has a significance of ≈ 3 and has a double-horned structure which is confirmed by the other method.Note that we report the period of the 'horn' of this broad peak that has the higher significance.The double horned structure is also seen in the color-color diagram of the wavelet of this Cycle 2 observation (lower left panel of Fig. 4 at 3.78 days (higher peak).Interestingly, only one of the peaks is present in the first half of the observation but as the observation evolves with time, the second one appears.So, the behavior of this double-horned peak could only be distinguished when the observation is decomposed into time and frequency domains simultaneously which is achieved by using WWZ.The feature at the period of 5 days is relatively stronger at the beginning of the observation but gets weaker with time and becomes negligible by the end of the observation.All these methods show the most significance for a possible QPO at around 3.8 days.
The GLSP for the Cycle 4 data of 1RXSJ111741.0+254858 is plotted in the lower right panel of Fig. 3.The most prominent peak is found at the period of 6.20 days (although that is over 25% of the length of the observations, so not convincing) followed by the peaks at 4.38 and 2.86 days.The perseverance of these peaks can be assessed by the wavelet analysis which is plotted in the lower right panel of Fig. 4 for Cycle 4 data.The feature at 6.25 days, which is more than 3 significant, seems to be strongly persistent throughout the observation, though it corresponds to too few cycles to be strongly supported.The features at 4.64 days continue to present strongly throughout the observation which is also the case for the somewhat weaker signal at 2.87 days.So, the most promising possible QPOs indicated by both GLSP and WWZ methods are around 6.2, 4.5, and 2.9 days, and the latter would correspond to over 9 cycles.

Caveats
One of the main challenges encountered while analyzing TESS light curves is the presence of gaps between light curves which are not negligible in comparison to the duration of the light curves.For instance, a gap of 4 days is present in the 25-day light curve of Cycle 2 of 1RXS J111741.0+254858.To assess the effect of gaps on the power spectrum, we have simulated 1000 light curves each, with and without gaps, and with initial conditions similar to that of Cycle 2 observation of 1RXS J111741.0+254858.The left panel of Fig. 5 shows the 1 confidence bound for the periodogram for the 20-day light curves both without gaps and with a gap of 4 days, assuming a bending power law as the underlying model.We could deduce from the plot that the periodograms in both cases are consistent with each other at higher frequencies.However, as we go towards lower frequencies (  ≈ 0.1d −1 ), the periodogram distributions no longer closely agree with each other.The gap between the light curve segments can affect its power spectrum at frequencies lower than 0.1 d −1 .So, we restrict our analysis to frequencies above 0.1 d −1 to avoid this issue.We note that the effect of a gap is more pronounced in the case of simple power law model for the PSD (left panel of Fig. 6).The confidence intervals for the light curves with gaps are more stringent as compared to continuous light curves at low frequencies (  ≈ 0.07d −1 ) whereas at higher frequencies, the confidence bounds for continuous light curves are less stringent.
A second issue is the possibility of non-stationarity of the data which affects the power spectrum and the resultant underlying model.To address this, we performed wavelet analysis to assess not only the significance of the signal but also its perseverance.So we have simulated 25-day light curves with a 4-day gap in-between its halves for two cases: stationary, with the same statistical properties throughout the light curve; and non-stationary, where the statistical properties are different for the data before and after the gap.For these simulations, we took the same light curves as before (Cycle 2 observation of 1RXS J111741.0+254858).The middle panels of Figs. 5 and 6 show the 1 error bound for stationary and non-stationary light curves with gaps for simple power law and bending power law respectively.In both cases, the bounds are quite similar at frequencies above  ≈ 0.2d −1 and only start to significantly diverge below  ≈ 0.1d −1 .
The third challenge is the choice of the underlying model.In this work, we employed both simple power-laws and bending power-laws which could provide incomplete descriptions of the data.By limiting the analysis to frequencies above 0.2 d −1 , we assumed that the power spectrum behaves like the simple power law (or bending power law in some cases) in this regime.
It is, however, also possible that a single model could not explain the behavior throughout the light curve.So, we fit the segments of the light curves separately and then simulate the light curves using the best-fit parameters for each segment.In this case, the confidence intervals with the same spectral index in the underlying bending power law (right panel of Fig. 5) are somewhat less tightly bound as compared to the different spectral index cases at higher frequencies; however, at lower frequencies, the confidence intervals are actually tighter for the same case.In the case of simple power laws (right panel of Fig. 6), the error bounds agree with each other quite 1RXS J004519.6+212735Sec. 17 closely throughout, though the different power-law case is somewhat broader.It can be concluded from these simulations that if we use the same spectral index throughout the observation, the confidence bound is underestimated as compared to the case where different spectral indices are used for different segments of light curves.So, it is preferable to simulate the segments within light curves separately if the gap is large in comparison to the resolution of the data.
The timescales found in this work are in the range of 3-5 days, corresponding to frequencies of 0.20-0.34d −1 so the caveats discussed in this section should have little impact.

Segment-wise analysis
The stochastic processes occurring in the jets and accretion disks of a blazar may be reflected in observations taken throughout the whole electromagnetic spectrum.In general, they may be non-stationary as their statistical properties can change with time.It is possible that the light curves of blazars may be non-stationary in nature even on the relatively short timescales investigated in this work.Indeed, evidence of non-stationarity has been observed in X-ray light curves of blazars and Seyferts with similar durations as the TESS light curves explored here (e.g.Alston 2019;Bhattacharyya et al. 2020).However, the simulation procedure we have used to this point to generate significance curves assumes the simple or bending powerlaw description of the PSD is constant throughout the observation and therefore, should be considered cautiously.While the observations  In each plot, the black curve is the GLSP.The green and blue curves show the 3  significance level for the simple power law and the bending power law, respectively.
of BL Lac have only a short (∼ 1d) gap, the observations of 1RXS J004519.6+212735 and 1RXS J111741.0+254858 have substantial gaps.It is important to examine whether any QPO signal is present both before and after the gap or if it is only present in one of the segments.To investigate this issue, we have divided the light curve into two segments separated by a gap in the middle and performed the GLSP analysis assuming both simple and bending power laws as the underlying red noise model.The portion of each light curve before the 4-day gap is termed as Segment 1 and that after the gap is termed as Segment 2.
The left panels of Fig. 7 show the segment-wise analysis for the Sector 17 observation of 1RXS J004519.6+212735.In Segment 1, a peak at around 6.0 days is found with at least 3  significance using the simple power-law but no peak is found significant employing a bending power-law.For Segment 2, a similar signal at around 5.9 days is significant at > 3 using the simple power law.While these peaks are not found to be significant using the bending power law model, their GLSP significance increases to 3 when both segments are considered together (see Fig. 3) and is further supported by the wavelet analysis (Fig. 4).
We analyzed the Sec.22 observation of 1RXS J111741.0+254858 by segments and the GLSP results are shown in the middle panels of Fig. 7. Signals of at least 3  significance are found at 4.5 days for Segment 1 and 2.8 days for Segment 2, assuming either simple or bending power laws as the underlying model.In the combined analysis (Fig. 3), these two peaks blend into a double-horned peak.
This transition from one frequency to another is also evident in the color-color wavelet analysis of Sec.22 (Fig. 4).
For the Sec.49 observation of 1RXS J111741.0+254858, the segment-wise GLSP results are plotted in the right panels of Fig. 7.In Segment 1, a period of around 1.8 days is found to be significant at 3  with the simple power law, though not with the bending one.This signal is also seen in the combined wavelet analysis of this observation (Fig. 4) but with less significance, as it is relatively strong only in the first half of the observation.In Segment 2, a very strong signal with at least 3  with both simple and bending power law is found at a period of around 6 days.This evolution of the possible ∼6 day period is shown in the combined wavelet analysis (Fig. 4) where the power of the signal increases with time.

DISCUSSION
Emission from blazars is believed to be dominated by non-thermal processes in the relativistic jet.As the jet is viewed at a very small inclination angle (< 10 • ) (Urry & Padovani 1995), the emission from the jet is magnified enormously by the relativistic effects and outweighs the cumulative emission from the host galaxy and accretion disk (e.g.Jorstad et al. 2001).The emission, extending from radio to -rays, is observed to have variability on timescales of a few minutes (e.g.Ahnen et al. 2017;Ackermann et al. 2016;Jorstad et al. 2022) to years (e.g.Villforth et al. 2010;Acciari et al. 2014).Some variability timescales have also been claimed to be QPO signatures in blazars (Hayashida et al. 1998;King et al. 2013;Graham et al. 2015;Tripathi et al. 2021;Jorstad et al. 2022).The short timescales of much of the variability indicate that the emitting region must be compact.

Kink Instability
Very recently, Jorstad et al. (2022) reported a quasi-periodic feature with a 0.57 day period in BL Lacertae, one of the objects analyzed in this work, employing both optical and -ray observations that were taken between March and December 2020.They considered a kink instability as the reason for the origin of these QPOs based on the approach given by Dong et al. (2020).In this model, periodicity is the result of a type of current-driven plasma instabilities known as kink instabilities.Relativistic magnetohydrodynamic simulations have shown that these kinds of instabilities can arise in the jet having a strong toroidal magnetic field and result in both twisting of field lines and an increase in particle acceleration (Mizuno et al. 2009;Barniol Duran et al. 2017).This distorted magnetic field grows to produce a quasi-periodic kink in the jet that manifests into a region of compressed moving plasma with enhanced emission that could be observed to have a QPO signature.The timescales of such QPO signatures in blazars can be of the order of days and are related to the growth of kink instabilities over timescales of weeks to months.Jorstad et al. (2022) performed semi-analytical simulations to model the time-dependent emission from the plasma following Dong et al. (2020).This model assumes the dominant source of fluctuations arises from a blob in the kink region and has three components: a constant toroidal magnetic component; a sinusoidal poloidal component having temporal period  and amplitude     ; and a turbulence profile having amplitude    which averages to   ;   corresponds to the average emission power during the observation.They constrained these parameters along with the normalization using MCMC simulations and found the overall magnetic (and hence polarization) period to be 1.14 days, which corresponds to the observed flux variability QPO of 0.57 d.They attributed the 4-day QPO period also detected using the wavelet method to the growth of this one-day period kink with time.
Quantifying the transverse motion of the kinked region can be used to calculate its growth rate (Mizuno et al. 2009).Dong et al. (2020) estimated the kink growth time (   ) to be the ratio of the transverse displacement of the jet from its center (   ) to the average velocity of propagation (⟨  ⟩).They found the growth time of the kink instability to be consistent with the period of the QPOs estimated in their simulations.The period associated with the kink instability (  ) in the observer's frame is given by (Dong et al. 2020) where  is the Doppler factor and    is the size of the emission region in the co-moving frame and its value can be 10 16 -10 17 cm for a typical BL Lac object.In Dong et al. (2020), ⟨  ⟩ is estimated to be ≈ 0.16.So, the kink instability is related to the properties of the jet which include its inclination, velocity, hence the Doppler factor, and size of the emission region.Calculating   for BL Lacertae using  = 15 yields a time of 1-10 days which is consistent with the QPOs timescales found in this work.So it is reasonable to suggest that the apparent quasi-periods of 3-5 days we obtained here for all three sources using different timing analyses can be attributed to the growth of the quasi-periodic kinks originating in the inner portions of jets.
The longer QPO periods of around 3-5 days, indicated by both GLSP and wavelet approaches could be explained by the kink evolving with time.Any quasi-periods of more than ∼ 5 days cannot be firmly detected with the TESS light curves having a duration of 25-27 days as a period of more than 6 days corresponds to fewer than 4 cycles.It should also be noted that any QPOs of periods shorter than a few hours could not be reliably detected in these data sets as the sampling frequency of these TESS data is 30 minutes so they would be very close to the Nyquist frequency and lie in the white-noise portion of the power spectrum.
Still, there are other processes that might be able to explain any QPOs on the scale of a few days.Possible explanations include the normal modes of oscillations trapped in the inner part of the accretion geometry due to the strong gravity of the central compact object (e.g., Perez et al. 1997;Espaillat et al. 2008).These modes correspond to the internal oscillations of the accretion disk and could produce the observed QPOs.Magnetorotational instability (MRI) turbulence in the accretion flow (e.g., Abramowicz et al. 2004) might also result in such periodicity.The frame-dragging around a rotating SMBH leads to the Lens-Thirring precession near the inner portion of the accretion disk which could naturally produce a QPO with a period of the order of a few weeks (Stella & Vietri 1998).However, these processes are more likely to be more important in Seyfert galaxies and quasars where the emission from the disk dominates over that of the jets, whereas the jet emission is the most important in these blazars.

BL Lacertae
BL Lacertae is the prototype of the BL Lac objects characterized by very weak or absent broad emission lines and high variability observed throughout the electromagnetic spectrum.This object has been observed extensively in all wavebands, from radio to -rays (e.g.Sandrinelli et al. 2018;Raiteri et al. 2003;Ren et al. 2023, and references therein).It also has been claimed to show quasi-periodic  behavior in various wavebands with periods ranging from a few minutes to hours, weeks, and even years.Villata et al. (2004) claimed a periodicity of 8 years for radio outbursts for this source which was later supported by Villata et al. (2009).Sandrinelli et al. (2017) reported a periodicity of around 700 days in both optical and -ray observations.As mentioned above, Jorstad et al. (2022) analyzed the optical and -ray observations and claimed a periodicity of around 13 hours in both.They found a very strong correlation between the optical and -ray light curves (correlation coefficient: 0.62±0.04)with a temporal lag between them consistent with 0. The detection of essentially identical periodicites implies that both optical and -rays are produced from the same mechanism which is claimed to be the kink instability in the jet.They also found a periodicity of the order of a few days which is explained by the temporal growth of the kink.
In our work, the periodicity of ∼5 days is found in the TESS light curves of BL Lacertae, which we suggest is the optical counterpart of the kink developed in the jet.BL Lacertae was observed by TESS for 27 days beginning on 12 September 2019, approximately 5 months prior to the onset of the observations analyzed by Jorstad et al. (2022).Ren et al. (2023) analyzed the weekly and monthly binned Fermi-LAT observations of BL Lacertae taken between August 2008 and April 2021 and found no significant long-lived periodicity.Although the TESS observation analyzed in this work is subsumed within the span of these Fermi observations, the cadences are significantly different.The minimum cadence normally available with Fermi data is 1 day and they usually need to be binned over longer periods to get adequate statistics.Therefore, the Fermi observations with a 1 day cadence could not resolve periodic variability with timescales of just a few days, while the 30 minute cadence of the Cycle 2 TESS data can.
We have applied the kink instability model to the TESS observation of BL Lacerate and followed the approach of Jorstad et al. (2022) to fit our 26-day observation using their publicly available code6 .Fig. 8 displays the corner plots showing the posterior distributions of various parameters using the algorithm applied in the code.We used the Doppler factor  = 15 which is derived in Jorstad et al. (2001) for BL Lacertae using Very Large Baseline Array (VLBA) observations at 43 GHz.We used 100 walkers with 120000 iterations and fit 5 parameters used in the code.We obtained the period of 1.16 days which agrees nicely with the results of Jorstad et al. (2022).
Motivated by the successful application of the kink instability model to BL Lac's TESS data, we applied this formalism to these other two sources.For 1RXS J004519.6+212735, the period is found to be 1.20 days.For the Cycle 2 and Cycle 4 observations of 1RXS J111741.0+254858,periods of 1.09 and 1.1 days are found, respectively.The corresponding MCMC corner plots are given in Figs.8-9.For this analysis, we used the Doppler factor appropriate for BL Lac itself,  = 15 for both sources.We also examined whether using different values of  would yield different temporal periods.Hence we also used the average Doppler factor value of 9.2 derived in Liodakis et al. (2017) for BL Lac objects.For 1RXS J004519.6+212735, the new period is found to be 1.26 days, while periods of 1.19 and 1.21 days are found for the Cycle 2 and Cycle 4 observations of 1RXS J111741.0+254858,respectively.So, plausible variations in  have only modest effects on the observed period of the kink developed in the jet.

CONCLUSION
In this work, we have examined densely sampled optical light curves of three blazars, BL Lacertae, 1RXS J111741.0+254858, and 1RXS J111741.0+254858.These were extracted from TESS using a customized reduction approach designed to preserve actual AGN variability.We have discovered tentative evidence for QPOs in these light curves on the timescale of a few days.For each light curve, these possible QPOs mainly lie between ∼ 2 and ∼ 6 days.We also ran MCMC simulations to examine the claimed QPOs present in the observations, assuming they originate with a kink instability in the jet.We found the periods of the kinks to be consistent with the period reported for BL Lacertae in Jorstad et al. (2022).Our analysis thus supports the kink instability hypothesis where the apparently observed QPOs can be produced through the temporal evolution of a kink developed in the jet.

Figure 1 .
Figure 1.TESS light curves analyzed in the work.The blue dots represent the fluxes and the gray curves denote the associated errors, while the red curves show the running averages.Upper panel: The left plot is of BL Lacertae and the right one is of 1RXS J004519.6+212735.Lower panel: light curves of 1RXS J111741.0+254858,where the left and right plots respectively are for the Cycle 2 and Cycle 4 observations.

Figure 3 .
Figure 3. Generalized Lomb-Scargle Periodograms for the TESS light curves analyzed in the work.The black histogram is the power spectral density, the blue curve shows 3 significance for simple power law and green curves shows the same confidence for bending power law.Upper panel: BL Lacertae (left) and 1RXS J004519.6+212735(right).Lower panel: Cycle 2 (left) and Cycle 4 (right) observations of 1RXS J111741.0+254858.

Figure 4 .
Figure 4. Weighted wavelet Z-transform (WWZ) analysis results for the light curves analyzed in the work.In each plot, the left sub-panel corresponds to the color-color diagram of wavelet power (red denoting the maximum power and decreasing towards violet and black) plotted against time and frequency.The right sub-panel shows the time-averaged WWZ plotted against frequency in black; the red dotted lines represent 3 significances.Upper panel: BL Lac (left) and 1RXS J004519.6+212735(right).Lower panel: 1RXS J111741.0+254858 in Cycle 2 (left) and Cycle 4 (right).

Figure 5 .
Figure 5. Simulated results demonstrating the effects of caveats for a bending power law as the underlying model.Left: 1  confidence region for the periodogram obtained by simulating the 20-day light curves; with a gap (blue) and without a gap (red) of 4 days.Middle: 1  confidence region for the power spectrum obtained by simulating stationary (blue) and non-stationary (red) light curves.Right: 1  confidence region for the periodogram obtained by simulating light curves having the same (blue) and different (red) PSD index.For all of these simulations, we used the inputs from the Cycle 2 observation of 1RXS J111741.0+254858.

Figure 6 .
Figure 6.Same as Fig. 5 but for a simple power law as the underlying model.

Figure 7 .
Figure 7. Segment-wise GLSP analysis for the light curves of 1RXS J004519.6+212735(left panel) and 1RXS J111741.0+254858(middle and right panels).In each plot, the black curve is the GLSP.The green and blue curves show the 3  significance level for the simple power law and the bending power law, respectively.

Table 1 .
Best fit parameters and associated errors for the bending power-law sampled from MCMC realizations.

Table 2 .
Best fit parameters and associated errors for the simple power-law sampled from MCMC realizations.