Rapid Optical Flares in the Blazar OJ 287 on Intraday Timescales with TESS

We have analyzed the optical light curves of the blazar OJ 287 obtained with the Transiting Exoplanet Survey Satellite over about 80 days from 2021 October 13 to December 31, with an unprecedented sampling of 2 minutes. Although significant variability has been found during the entire period, we have detected two exceptional flares with flux nearly doubling and then nearly tripling over 2 days in the middle of 2021 November. We went through the light-curve analysis using the excess variance, generalized Lomb–Scargle periodogram, and continuous autoregressive moving average methods and estimated the flux halving/doubling timescales. The most probable shortest variability timescale was found to be 0.38 days in the rising phase of the first flare. We briefly discuss some emission models for the variability in radio-loud active galactic nuclei that could be capable of producing such fast flares.


INTRODUCTION
The blazar subclass of radio loud (RL) active galactic nuclei (AGN) displays flux, spectral and polarization variability throughout the electromagnetic (EM) spectrum on diverse timescales ranging from a few minutes to several years.Blazar emission is predominantly nonthermal.Probably the least well understood temporal variability is observed on timescales of minutes to several hours and commonly called as microvariability (Miller et al. 1989) or intraday variability (IDV) (Wagner & Witzel 1995).Microvariability in optical flux provides a strong route toward understanding the physical processes occurring in the most compact emitting regions of blazars.Miller et al. (1989) made the pioneering discovery of optical microvariability in the blazar BL Lacerate, and in the subsequent ∼3.5 decades, many of the brighter blazars have been observed for the study of microvariability during thousands of observing nights using var-amp700151@gmail.com (SK), acgupta30@gmail.com(ACG), wiitap@tcnj.edu(PJW) ious telescopes (e.g., Miller et al. 1989;Carini et al. 1992;Wagner et al. 1993;Heidt & Wagner 1996;Sagar et al. 1999Sagar et al. , 2004;;Qian et al. 2002;Montagni et al. 2006;Gupta et al. 2008;Poon et al. 2009;Gaur et al. 2012;Agarwal & Gupta 2015;Goyal et al. 2018;Pandey et al. 2020;Kalita et al. 2021;Raiteri et al. 2021;Dhiman et al. 2023;Pininti et al. 2023;Wehrle et al. 2019, 2023, andreferences therein).In a statistical study of optical microvariability properties of various classes of AGN, Gupta & Joshi (2005) found that if a blazar is observed continuously for less than 6 h, the chance of seeing micro-variations is ≈60−65%, but it rises to 80−85% for nightly observations that exceed 6 h.
Here we consider data taken on OJ 287 over a span of about 80 days (from 2021 October 13 to 2021 December 31) by the Transiting Exoplanet Survey Satellite (TESS) 1 .Our focus is on a continuous and uniformly sampled set of observations, more precisely the first segment of sector 45 that spans ∼12 days with time resolution of 2 minutes.These data provide us an excellent opportunity to study optical microvariability of a blazar with essentially uniform sampling made at the 1 http://tess.gsfc.nasa.govshortest time resolution and over a quite extended duration.We found a double-peaked strong flare in the light curve (LC) that spanned about 2 days.We are unaware of a previous clear case for such a well-resolved double-peaked flare on such timescales.OJ 287 was earlier observed from the Kepler satellite in its K2 mission phase for a continuous period of 75 days (2015 April 27 to 2015 July 10) and this uniformly sampled optical LC displayed several significant flares (Goyal et al. 2018;Wehrle et al. 2019), but none as fast with as large an amplitude as those we see in these TESS data.Wehrle et al. (2023) analyzed a subsequent K2 observation of the source taken for 51 days (from 2018 May 13 to 2018 July 2) and also reported multiple rapid, but not quite as strong, flares during this epoch.
In Section 2, we describe data acquisition and reduction.
In Section 3, we explain the data analysis techniques we used and present results.A discussion is provided in Section 4.

DATA ACQUISITION AND REDUCTION
Aside from some data gaps of 1-2 days related to telemetry, OJ 287 was observed essentially continuously by TESS; its detector bandpass spans the range 600 -1000 nm and is centered at the traditional Cousins I-band2 .These observations spanned ∼80 days across sectors 44, 45 and 46, with a cadence of 2 minutes.We have used the PDCSAP_FLUX (Jenkins et al. 2016) values for reduction as described in Kishore et al. (2023) and have followed the reduction procedure given there.Hence, we refer the readers to that paper for the nomenclature used in this section and the optimum values found for the three parameters employed during data reduction: the two goodness metrics (overfitting and underfitting) and the regularization factor (α).
Unfortunately there were no LC files for nearby comparision stars (which contain SAP and PDCSAP fluxes) from the TESS pipeline already available, because these values are made available only for pre-selected targets.TESS also created full-frame images (FFIs) of the portions of sky observed at a cadence of 20 minutes for those sectors, so we used these FFIs to extract the SAP LCs of comparison stars 4, 10, 11 (Smith et al. 1985) along with that of OJ 287 to confirm that our source's variability is intrinsic.Comparison stars 4 and 10 had substantially higher fluxes than OJ 287, so, in Figure 1 Sector Sector Sector  we used star 11, with a very comparable brightness, to illustrate the genuine variability of the source.
It should be emphasized that Fig. 1 provides a preliminary version of the LC that provides a useful visual comparison, but a more fully reduced version has been used in our actual analysis of the LC of OJ 287.
Table 1 includes the values of these fitting parameters obtained for reduction of PDCSAP_FLUX of our object for each of the Sectors.As discussed in Kishore et al. (2023), the values of both the underfitting and overfitting goodness metrics should be kept at or above 0.8, consistent with the lowest possible α value.For sector 45, the overfitting metric was below 0.8 at α = 0.1, so the optimum value of α was found to be 0.34 when the overfitting metric crosses the critical limit of 0.8.Fig. 2 illustrates these conditions for the sector 45 LC reduction.Fig. 3 includes two plots for each of the Sectors.The upper panel in each set contains the raw LC and the different co-trending basis vectors used for the reduction, while the lower panel compares the raw LC and reduced LC.Variability is seen in all three of these LCs, but it is much stronger in Sector 45.In Fig. 4 we show the zoomed in Sector 45 reduced LC, showing the double-flare in the first segment of this sector in more detail.

Excess variance and flares
Flux variability is one of the fundamental properties of a typical blazar across all EM bands; however, an addi- tional variance is present in all LCs from the measurement errors in the observations.The excess variance method incorporates the measurement errors as well in quantifying the variability.The fractional root-meansquare variability amplitude, that is the square root of normalized excess variance, and corresponding uncertainty, have been computed as in Vaughan et al. (2003).
In our analysis of F var for the two flares in the first segment of Sector 45, the average flux counts after removing the two flares from the segment (where the determination of the start and end of the flares was done visually, selecting the minima in the LC nearest to the flaring span) was used as the baseline flux.This value came out to be ∼193 e − /s.The maximum flux values during the two flares are ∼349 e − /s and ∼566 e − /s, respectively, corresponding to nominal increases of ∼81% and ∼194%, respectively.With this baseline, the LC shows an overall F var of 72% during the total flaring period.Table 2 includes the F var obtained for the two flares along with the other parameters describing the flares.

Variability timescales
We have used the halving/doubling timescale given as where τ and F (t) are the characteristic halving/doubling timescale and the flux value at time t.In our analysis, we first divided each flare into two parts: a 'rising phase' and a 'declining phase' and then examined each and every possible distinct pair of data points within each flare.We then imposed a selection criterion so that only those data point pairs, where the differences in flux were greater than 3σ, were considered for the estimation of τ (Foschini et al. 2011).The doubling/halving time was estimated for each obtained pair.Thus, for each of the two phases of the two flares, the obtained ensemble of distinct pairs led to a distribution of the timescale.Fig. 5 shows the distribution of the τ estimates obtained for the two pairs of rising and declining phases and Table 2 includes the most probable timescales with the corresponding uncertainties (after fitting the distribution with a Gaussian function) for each of the four flaring phases, as well as the times of the two peaks.It should be noted that there is a tiny positive fluctuation during the declining phase of the LC during the first flare between epochs 2529.921 to 2529.953; so we dropped this small rising fluctuation from the declining phase while evaluating τ .The two portions of the declining phase, from 2529.740 to 2529.921 and from 2529.953 to 2530.160, led to separate sets of τ values.These two sets of τ values were combined to give the composite timescale distribution of this declining phase, shown in the upper right panel of Fig. 5.

Periodograms
Apart from the variability timescales that may be associated with the size of the emitting region, the spectral index of the periodogram, or the power spectral density (PSD) slope, can yield information about the source of the variability.Various physical models naturally yield somewhat different ranges of spectral indices (e.g., Pollack et al. 2016;Wehrle et al. 2019).We have used the generalized Lomb-Scargle periodogram as in Kishore et al. (2023) for the analysis of all six segments of the  three sectors, as well as the first flare and the second flare of Sector 45 individually.As usual, the PSDs so obtained flatten to instrumental white noise in the high frequency regime, so we used cut-off frequencies (Gierliński et al. 2008;Lachowicz et al. 2009) for each segment, which varied somewhat, while conservatively fitting the PSDs in the red noise region with power laws (P (ν) = Aν α ).To test the power-law fit to the PSD, we followed the approach of Vaughan (2005).This involved considering twice the ratio of the periodogram to the fitting model at each frequency.The set of these ratios were used to form a cumulative distribution function (CDF).Ideally this CDF should follow that of a χ 2 2 (χ 2 distribution with two degrees of freedom) (Vaughan 2005).Under the null hypothesis that these two distributions are the same, Kolmogorov-Smirnov (KS) tests were performed for each PSD fitting.The p-values (probability of not discarding the null hypothesis) were evaluated for each comparison of segment-wise fits, and the high p-values we obtained (given in Table 3) do indicate good PSD fits.The scipy3 python package was used for the goodness of fit estimation, following the steps described in Vaughan (2005).An oversampling by a factor of 5 has been employed in the Lomb-Scargle periodogram calculation as there is a paucity of datapoints in the low frequency red-noise region.
Fig. 6 displays the PSD behaviors of the first segment of the Sector 45 LC and Table 3 gives the spectral indices we obtained for all six of the segments.The PSD slopes have typically high values (close to or greater than 2) up through the segment including the flares, but following that, there seems to be a considerable flattening of the PSD slopes in the later segments.The overall PSD slopes during the flares and during the entire segment of Sector 45 agree within the errors.
The PSD slope values listed in Table 3 for the several segments cover quite a wide range.Hence they unfortunately cannot be used to eliminate any of the rather small number of models for blazar variablity that have evaluated the resulting PSDs (e.g.Pollack et al. 2016;Wehrle et al. 2019;Kadowaki et al. 2021), each of which is, or appears to be, capable of yielding a span of slopes within this range.

CARMA Modeling
Although the great majority of the analyses of AGN LCs in the literature have focused on periodogram slopes, a more sophisticated approach to analyzing the structure of LCs involves considering autoregressive (AR), moving average (MA), models, which in their continu-ous version are given the acronym, CARMA (e.g. Kelly et al. 2009Kelly et al. , 2014;;Kasliwal et al. 2017;Goyal et al. 2018).CARMA models employ the plausible assumption that a LC is a realization of a Gaussian noise process.Specifically, a CARMA(p, q) model connects the LC and its first p time derivatives to the noise and its first q time derivatives (for the definitions used here, see Kelly et al. (2014), Equation (1) and accompanying text).A CARMA(1,0) model is equivalent to a damped random walk, or Ornstein-Ulenbeck process, and seems to better describe the long term LCs of many quasars than does a single periodogram slope (e.g. Kelly et al. 2009Kelly et al. , 2014)).
A physical interpretation of this approach is that the AR part of the model describes the short-term memory in the system, while the MA part indicates how the amplitudes of random perturbations behave on different timescales.Hence, both the correlation structure and degree of smoothness of noisy processes can be described by CARMA models (e.g.Moreno et al. 2019).Ryan et al. (2019) found that blazar γ-ray LCs were usually better fit by the modestly more complex CARMA(2,1) models than by CARMA(1,0) ones, though the number of objects for which this analysis could be done was modest.In a paper describing a very wide range of approaches to analysing Fermi-LAT γ-ray LCs of 11 blazars (Tarnopolski et al. 2020) perform CARMA modeling of LCs binned into 7-day, 10-day and 14-day intervals.Their results are basically consistent with those of Ryan et al. (2019); though both CARMA(1,0) and (2,1) models were most frequently optimal, occasionally (3,0) or (3,1) models were preferred.
We have performed CARMA analyses of these TESS LCs following the approach of Yu & Richards (2022)4 which provides a fast way to produce CARMA models.All (p, q) pairs with 1 ≤ p ≤ 5 and q < p were considered for each of the 6 segments of the 3 Sectors.After randomly generating the initial parameter values for each of the CARMA (p, q) models, the LCs were fitted by those models.The goodness of each fit was estimated by finding the log-likelihood value with respect to the initial LC, and the best fit was taken to be the one that maximized that quantity.The resulting best fitting (p, q) values are given in the last column of Table 3.We see that both Sectors 44 and 45 can be characterized by a CARMA(1,0), or a damped random walk model, though the last Sector 46 prefers models with p = 2.

DISCUSSION
By comparing all the TESS observations of OJ 287 from Fig. 3, it is clear that the double peaked flare in sector 45 illustrates a strong outburst not otherwise seen during these three Sectors.From the subplot in Fig. 4, it is evident that the blazar flux rises during the flares almost monotonically, but the decays involve some jerks.We also notice that the rises are faster than the decays during each of the two flare phases.Chiaberge & Ghisellini (1999) showed that symmetric light curves, with similar rise and decay timescales (Abdo et al. 2010), are expected when the cooling time of electrons t cooling is significantly shorter than the light crossing time, R/c, with R the size of the emission region.In this analysis it was assumed that the relativistic electrons are accelerated so that their energies obtain a power-law distribution.Asymmetric profiles, with decay times longer than rise times, result when t cooling > R/c because then the time scale for decline is longer than the time scale for rise of the flare (Chiaberge & Ghisellini 1999;Li et al. 2018).
OJ 287 was earlier observed by Kepler for a continuous period of 75 days (2015 April 27 to 2015 July 10).A uniformly sampled optical LC of OJ 287 from the first K2 observation detected several rapid flares, though none as significant on these short time scales as seen in Sector 45 (Goyal et al. 2018;Wehrle et al. 2019).Goyal et al. (2018) found that the PSD (power spectral density) of the total LC at that time was well fitted by a CARMA (4,1) model5 , so of a higher order than we find during the later TESS observations.It may be worth noting, though, that a CARMA(4,2) or (4,3) model was among the two runners-up to the best fits we found for 5 of the 6 segments.Wehrle et al. (2023) found the PSD slope during that period to be −2.28 ± 0.17 for long-cadence (30 minute bins) data, whereas the PSD slope was found to be somewhat steeper, −2.65 ± 0.05, for the short-cadence (1 minute bins) data.In the later K2 observation of OJ 287 that lasted 51 days (2018 May 13 to 2018 Jul 2), the flux variability shows a similar jagged behavior as during its 2015 observation.The PSD slopes for long-cadence (30 minutes bin) data and short-cadence (1 minute bin) data were found to be −1.96± 0.20 and −2.26 ± 0.06, respectively (Wehrle et al. 2023), and so in the same range as the first three segments of the TESS observations discussed here.
Although the variable flux from blazars is dominated by jet emission, accretion disks can contribute in the low-flux state of flat spectrum radio galaxies.Since OJ 287 is a BL Lac and during this time it was in an overall intermediate brightness state (Gupta et al., in preparation) we can rule out the strong flux variation as arising from the accretion disk.Much of the optical (and other) variability in blazars can be explained by turbulence behind the shocks in a relativistic jet (e.g., Marscher 2014;Pollack et al. 2016) or in turbulence produced by magnetic reconnection (e.g., Kadowaki et al. 2021;Guo et al. 2021).These magnetic reconnection structures in jets can lead to very fast emission changes that typically produce flares with somewhat longer decay than rise times (Kadowaki et al. 2021), as we saw in OJ 287.The very rapid and substantial flares on sub-day timescales seen here could also arise from extremely compact regions with very high Doppler factors, as occurring in the mini-jet or jet-in-jet scenarios (e.g., Ghisellini & Tavecchio 2008;Giannios et al. 2009).In such models, portions of the plasma in the relativistic jets are accelerated to Lorentz factors ∼ 100 through magnetocentrifugal or magnetic reconnection processes and the resulting extreme Doppler boosting can yield fast and strong flux changes.
We obtained the shortest variability timescale (τ min = 0.38 days) in the rising phase of the first flare.To estimate an upper limit for the size of the emission region, R, we apply the simple causality constraint, where δ represents the Doppler factor.Cohen et al. (2018) compiled values of δ for OJ 287 from the literature: δ = 18.9±6 and 17.0 were derived from 43 GHz radio flares of OJ 287 in 1998-2000and 2003, respectively (Jorstad et al. 2005;Hovatta et al. 2009).A more recent value of δ = 8.7 was derived from millimeterwave flares after 2007 (Jorstad et al. 2017;Liodakis et al. 2017).By using the complete range of Doppler factor values, 8.7 to 18.9±6, taking z = 0.306 (Sitko & Junkkarinen 1985), applying the shortest variability timescale (τ min = 0.38 days), and using Equation (2), we estimate the size of emission region to be in the range of 2.2 × 10 15 cm -6.3 × 10 15 cm.
We have reduced and analyzed the TESS LCs for OJ 287 spanning three consecutive Sectors (44-46), which correspond to October 13 through December 31, 2021.Each Sector had a ∼2-3 day gap near its middle, so we analyzed the resulting 6 segments separately.All of them showed significant variability which can be approximately characterized as having PSD slopes in the range ∼ −1.5 to ∼ −2.5 and usually being well fit by CARMA(1,0) models.Such light curves are typical of AGN, but an unexpected result was the observation of two consecutive strong flares (with flux increases of ∼81% and ∼194%) seen during the first segment of Sector 45.Both of them had most probable doubling rise times around 0.4 d and decay times that were nearly as fast, indicating that this optical emission arises from a very compact region in the relativistic jet.

Figure 1 .
Figure 1.Partial Sector 45 LCs of OJ 287 and comparison star 11 during the strongest flares.

Figure 2 .
Figure 2. Goodness metric scan plot for the TESS Sector 45 observation of OJ 287.

Figure 3 .
Figure 3. Raw and reduced LCs of OJ 287 observed in all three sectors.The upper panels include the original PDCSAP fluxes and the co-trending basis vectors that are used to correct it and the bottom panels show the original and corrected LCs corresponding to each labeled sector.

Figure 4 .
Figure 4.The main plot includes the complete reduced Sector 45 LC.The subplot zooms in on the flare period; for better visualization of the trend, a 0.5 h binned LC has been overplotted.

Figure 5 .
Figure 5.Each panel shows the distribution of doubling/halving timescales corresponding to the labeled time spans.Positive values of τ indicate the rising phases and negative ones the declining phases of the flares.The orange curve in each panel shows the Gaussian fit determining the best value of τ .

Figure 6 .
Figure 6.The panels show the Lomb-Scargle periodograms corresponding to the labeled pieces of the first segment of the Sector 45 LC with their PSD slopes.

Table 1 .
Flux calibration details

Table 3 .
Segment-wise PSD and CARMA results