Observed Rate Variations in Superflaring G-type Stars

Flare catalogues for highly flaring stars in the TESS Input Catalog (TIC) are available, e.g. Günther et al. (2020); Tu et al. (2021). As we are attempting to identify G-type stars which exhibit flare rate variation, we select the most flare-producing stars from these previous studies. The effective temperature $T_{\text{eff}}$ and surface gravity $g$ are obtained from the TIC if available. If a star does not have parameters listed in the TIC, we fall back to using data from GAIA Data Release 2 and note these stars for inspection. We ensure each chosen star is a G-type star by requiring $5100$ K $\leq T_{\text{eff}}\leq 6000$ K and $\log g>4.0$, consistent with the criteria used by Tu et al. (2021).

TESS provides data in the form of a Target Pixel File (TPF), which encodes the flux of each pixel over time, or as Simple Aperture Photometry (SAP) data, an integrated light curve containing the flux from the pixels within a defined aperture mask. TESS provides observations for around 20,000 stars per sector in 2-minute cadence, and following the start of the extended mission with Sector 27, also provides 20-second cadence for up to 1000 stars per sector (Barclay, 2020). Each observation sector is viewed for approximately 27.4 days, divided into two segments corresponding to the ecliptic orbit of the satellite.

For each target star, we download the TPFs in 2-minute cadence and perform a visual inspection to confirm the star is not a binary star. Then, data points which have quality flags 1, 2, 4, 8, 16, 32, 64, 128, 512, 4096, 16384 are removed from the data set. These flags correspond to systematic issues such as altitude tweaks and scattered light incidents, unexpected brightening events (e.g. argabrightening, cosmic rays), as well as impulsive outliers and anomalous data points (Twicken et al., 2020; Fausnaugh et al., 2019). These events are undesirable and can affect the detection of superflares by introducing similar large, impulsive structures into the lightcurve, and hence they are removed. Finally, we create the integrated light curve of the star by integrating over each pixel specified in the provided aperture mask.

In order to identify flares, we first need to model the variability present in the light curve, which includes both systematic and astrophysical components. By doing this, we can construct a model of the quiescent luminosity of the star and separate the remaining outliers from this model such as flares. In TESS, systematic non-physical variations are due to instrumental effects, including the differing efficiency of each camera. Astrophysical variability present in the light curve most commonly consists of rotational modulation – quasi-periodic variation due to the presence of starspots (Notsu et al., 2019).

We begin by identifying the gaps in the data, which separate periods of continuous observation consisting of at least 720 consecutive data points with no gaps larger than an hour between any data points. The data in each of the continuous observation periods is normalised by dividing by the median value of the flux in that period. Large excursions from the mean are masked out before calculating the median via an iterative ‘sigma-clipping’ process: values above 2.5 standard deviations away from the mean are masked out of the data set, over 10 repeated iterations. Long term systematic variation is then modelled by fitting a third-order polynomial to each continuous observation period. This polynomial fit is then subtracted from the light curve. Using this systematics-corrected light curve, we identify the presence of periodicity (e.g. associated with rotational modulation) using a Lomb-Scargle periodogram. A star is considered to have periodicity if a significant signal is detected (Lomb-Scargle power $>$ 0.1).

To model the remaining variability, we use the nonlinear iterative filter previously used by Aigrain & Irwin (2004) for processing light curves to identify exoplanet transits. An example of the application of this filter is shown in Figure 1. This filter has two parameters, nmed and nlin, which control the window sizes of the median filter and box-car filter used in the algorithm. From manual experimentation, we find that due to the variable length of flaring events (as opposed to the fixed time scales of exoplanet transit events), scaling the length of the median filter window by the detected period of the star generally results in an accurate model of the quiescent background, and avoids including parts of the flare in the background model. Specifically, we choose $\mathtt{nmed}=100\times T$ and $\mathtt{nlin}=30$ for stars with period $T$, and otherwise use $\mathtt{nmed}=200$ for stars with no detected period.

Using the detrended light curve, flares can be identified as tranisent increases in flux, as shown in Figures 1 and 2. We use the flare finding algorithm by Chang et al. (2015), requiring that flares consist of at least six consecutive detrended flux values which are all at least 1.5 standard deviations above the mean (and 1.6 standard deviations when accounting for the photometric error of each data point).

Detected flares may not necessarily have occurred on the target star – flux contamination from a flare on an unresolved neighbouring star might also lead to flare signatures in the lightcurve. By measuring the flux-weighted average position (the ‘photometric centroid’; Bryson et al. (2013)) over time, we can compare the measured centroid during a flare with the centroid when the star is not active. Large shifts in the centroid correlated in time with flares would provide evidence that a flare did not originate on the target star.

For each continuous observation period, we calculated the median centroid row and column of the star, and then subtracted this from the centroid time series. The standard deviation in the row and column centroid offsets were calculated, excluding datapoints marked as flares. Any flare candidates which exhibited an offset larger than 2 standard deviations in either their row or column centroids were marked for inspection.

For each flare candidate, we calculate the energy and estimate the full-width at half-maximum (FWHM) duration. Energy calculations are done following the steps of Günther et al. (2020), using temperature and radius estimates and uncertainties from the Tess Input Catalog (Stassun et al. (2019), TIC v8) if available, and GAIA-DR2 (Gaia Collaboration et al., 2016, 2018) measurements otherwise. To estimate the FWHM, we fit the classical flare profile developed by Davenport et al. (2014) to each flare, and then calculate the FWHM using the fitted model. In some cases this is not possible, in which case the FWHM is estimated by fitting a second-order polynomial to the rise and decay phases separately.

We test the robustness of our detrending and flare finding algorithms by injecting artificial flares into the light curve, and then recovering them. Using the classical flare profile from Davenport et al. (2014), we generate flares with a specified amplitude and FWHM. The FWHM of each flare is drawn from a uniform distribution between 1 minute and 1 hour, and the amplitude of each flare is drawn from a log-normal distribution with the same mean and standard deviation as measured from the detected distribution of flares. We require the standard deviation of this log-normal distribution to be at least 0.3 relative flux units. For distributions with low variability (e.g. not many flares or flares with similar amplitudes) this ensures a more complete sampling of the amplitude-space. In practice, we find that a majority of observed flares on stars have estimated FWHMs within our sampled range, with the few flares with FWHMs longer than 1 hour having a large amplitude as well.

For each continuous observation period, we inject 10 flares into the light curve and attempt to recover them using the detrending and flare-finding algorithms described in Sections 2.2 and 2.3. The generated flares are spaced with a buffer of 10 data points between any other flares, previously detected or previously injected, to avoid overlap. We considered an injected flare to be detected if the peak of any flaring candidate is within 10 data points of the time period where the injected flare was inserted.

After a sufficient number of injection and recovery cycles (at least 50,000 flares), we can estimate the recovery probability for flares based on their amplitude and FWHM. The recovery probability can be visualised in the form of a FWHM-amplitude-recovery probability surface, or as separate plots of the recovery probability as a function of FWHM or amplitude. Figure 3 shows an example of this recovery probability surface, and highlights the characteristic sharp decrease in the recovery probability at low amplitude and short FWHMs.

For the individual recovery probability-amplitude and recovery probability-FWHM plots, we model the distributions using a logistic function. The choice of a logistic function follows from both the observed form of the distributions, and also from the expectation that at a sufficiently high amplitude or long FWHM, flares should always be recovered, and similarly at a sufficiently low amplitude or short FWHM the flare should be indistinguishable from noise.

Using a fitted logistic function, we identify the amplitude at which the recovery probability is at least 75%. All flares which are above this threshold are considered to be real flares and are marked for use in the analysis of the rate of flare occurrence. The FWHM of the flares is not factored into this filtering process – we find that most flares which would be excluded from FWHM filtering are already removed due to the amplitude filtering.

Our statistical analysis focuses on identifying the rate of flaring for each star, and any variability in the rate. To achieve this we use the Bayesian Blocks algorithm (Scargle et al., 2013), which optimally determines a piecewise constant-rate model for a set of times of flare occurrence, consisting of a set of blocks (intervals with constant rate) and rates for each block. The two free parameter in the algorithm – ncp_prior, the prior for the number of change points, and fp_rate, the percentage of detections that may be false positives – are kept at their default values (fp_rate$=0.05$, ncp_prior calculated using Equation 21 from Scargle et al. (2013)). The flare data passed into the Bayesian Blocks data is treated as having no gaps in observation. This is achieved by calculating the length of each gap in observation, and shifting any time data after this gap by the length of the gap, essentially giving a measure of time in terms of the active observation time rather than in real time.

To assign uncertainties to the Bayesian Blocks model we follow Scargle et al. (2013), and perform a bootstrap analysis. We generate 1000 bootstrap samples of the flare data by sampling with replacement an equivalent number of flare event times. Then, the Bayesian Blocks algorithm is applied to each sample, and the rates and change points (the edges of each block) for each run are saved. Since both the number of blocks, their positions and their size can vary, we sample the rate of each run at 4096 equally spaced intervals in observation time. Then, we calculate the average (which we refer to as the ‘mean bootstrap rate’) and the standard deviation (the ‘uncertainty’) of the rates over the bootstrap runs at these intervals.

The difference between the mean bootstrap rate and the optimal block representation gives an indication of how frequently that block is identified in the bootstrap samples, as well as information on when blocks are typically detected in the samples and their relative sizes. The uncertainty gives a measure of the spread in the rates – how much the size of detected blocks varies between the bootstrap samples.

It is of interest to investigate the waiting-time distributions (WTDs) for the flares. For each star, we generate the WTD for all flares, and also for flares in each of the identified blocks. In order to be consistent with the Bayesian Blocks algorithm, waiting times are calculated using the active observation time used by the flare data which the Bayesian Blocks algorithm is applied to, where gaps in observation are ignored. The WTD constructed in this way may differ from the actual WTD for flares from the star – ignoring gaps in the data reduces longer waiting times which span the gaps, and can also introduce additional waiting times which are not present in the true distribution. For stars with a sufficiently high flaring rate (more than 2 flares per observation period), the constructed WTD should closely match the real WTD as the number of waiting times affected by ignoring gaps is only a small proportion of all the waiting times.

Following Wheatland & Litvinenko (2002), we use the Bayesian Blocks model to construct a model WTD for the flares, using equation (4).

We also construct the flare frequency distribution (FFD) for each star. To investigate whether the FFD varies with flare rates for superflares, we also constructed the FFD for only flares in each of the blocks identified by the Bayesian Blocks algorithm. We fit a simple power-law model to each FFD using a Bayesian method with a uniform prior, which is equivalent to maximum likelihood estimation (Wheatland, 2010). For some stars, there was a rollover at lower energies due to the incomplete recovery of lower energy flares, and for these stars we fit the simple power-law model above a minimum threshold energy, with the threshold being estimated by visual inspection of the distribution.

We also examine the rotational phase distribution of flares. Previous analyses of the phase distributions focused on the overall phases of all flares on a star, however we investigate also whether there is a relationship between flare rate and phase dependence. Similarly to the FFD and WTD, we plot the phase distribution not only for the whole set of flares, but also for each block identified by the algorithm. We use the period from the Lomb-Scargle periodogram to fold the times for the star, with the position of 0 phase being arbitrarily chosen as the first data point in the time series.

We construct FFDs for each star using the whole set of flares, and then FFDs individually for each identified block from the Bayesian blocks algorithm. Uncertainties are calculated for both the energy estimates and for the counting statistics, the latter using Equations 11 and 14 from Gehrels (1986) with a confidence limit of 1$\sigma$. We find that four of our rate variable stars have a rollover at lower energies, indicative of incomplete recovery of lower energy flares. This is expected due to the lower signal-to-noise ratio of these smaller flares due to their lower amplitudes and shorter durations. This is consistent with the recovery rates from the artificial flare injection and recovery procedure reducing with smaller amplitudes and durations. At high energies we also observe departures from a power-law. Similar breaks from the power-law distribution at high energies have been seen previously in G dwarfs (Davenport et al., 2019).

Figure 8 shows the FFDs for the seven stars which exhibited rate variation. A simple power-law model matched the distributions of three of the stars (TIC 43472154, TIC 85487971 and TIC 284358867), with the left panel showing the FFDs and power-law fits for these stars. The other four stars with rate variation (TIC 258771803, TIC 364588501, TIC 382575967 and TIC 394030788) exhibited significant departure from a simple power-law model, as shown in the right panel of Figure 8. These stars were fit using the simple power-law model with a minimum energy threshold. The power-law indices for these stars are higher than those fit with a simple power-law distribution, and are high compared to indices for solar flare distributions (for total energy, $1.62\pm 0.12$, Aschwanden et al. (2016)) as well as indices reported for superflares (e.g. $\gamma\approx 1.5$ (Maehara et al., 2015), $\gamma=1.76\pm 0.11$ (Tu et al., 2021)). In the cases of TIC 258771803 and TIC 382575967, which have power-law indices $\gamma\approx 2.7$, the large indices may be due to the energy thresholds being too high and capturing the tail-end of the power-law distribution with a departure at high energies due to incomplete sampling from limited observation time.

For the FFDs for each block, we attempted to fit the power-law models similarly to the FFDs for all flares. However, some of the distributions exhibit sharp departures from the power-law model. For the blocks we are able to fit successfully with power-laws, we find that the power-law indices for each block match the power-law index for the overall FFD. The fits for individual blocks have much larger uncertainties due to the lower number of flares.

Figure 9 shows the WTDs for the seven stars with rate variation, together with the piecewise-exponential model implied by the Bayesian Blocks results. There is a good qualitative agreement between the observed WTDs and the models, in particular for stars with more flares.

The WTD for each of the identified blocks and the exponential (Poisson) models corresponding to that block also showed agreement for each of the stars. Figure 10 shows the WTD of all flares and the individual blocks for both TIC 382575967 and TIC 364588501. The exponential models for each of the identified blocks (shown as the non-blue distributions) are in good agreement with the distributions, particularly at shorter waiting times. On stars with smaller numbers of flares such as TIC 258771803 and TIC 284358867, deviations between the distribution and models can be attributed in part to the larger uncertainties associated with smaller numbers of flares. In Figure 10 we have indicated the uncertainties for TIC 258771803. These deviations from the piecewise-exponential model may also be attributed to the Bayesian blocks representation not capturing all rate variability in the data. This is particularly prevalent for stars with a smaller number of events, as the detection of rate changes is more difficult.

We calculate the phase distribution for the flares of each star (both the overall phase distribution and the distribution in each identified block). The phase is calculated by taking the modulus of the time using the period detected from the Lomb-Scargle periodogram. The flare phases are grouped into 10 bins and plotted as a histogram. For most stars, the distributions of flares with phase are roughly uniform, both overall and for each block. However, TIC 364588501 and TIC 382575967 show possible phase dependence, as seen in Figures 11 and 12.

Our choice of phase is chosen arbitrarily based on the first data point, however we can further investigate the phase dependence of flares and any evolution in the rotational modulation by plotting the background flux model folded over each rotational cycle using the period detected with the periodogram. We refer to these diagrams as phase maps. These phase maps are similar to the diagrams used by Davenport et al. (2015) and Davenport et al. (2020). On these phase maps, we also plot the positions of flares in terms of their phase and their rotational cycle.

Figure 11 shows the phase map for TIC 382575967, and reveals a consistent darkening across many cycles at the same phase. This is indicative of a persistent structure on the stellar surface rotating with the same period as measured by the periodogram. The times of flares are indicated by the dark spots on the phase map, and the phase histogram for the flares is the right-hand panel. The histogram indicates some tendency for flares to occur at a phase corresponding with maximum brightness in the rotational modulation. Comparison of the phase distribution to a uniform distribution using a Kolmogorov-Smirnov (K-S) test gave a $p$-value of 0.076. This $p$-value is very low compared to other phase distributions, and provides some evidence that the phase distribution is not well fit with a uniform distribution.

For several stars (including TIC 364588501) the histogram plots did not match up with the expected phase distribution when checking the light curve directly. We observed for these stars that the rotational modulation occurred with a different periodicity than the dominant periodicity detected through the periodogram, and this was evident through the phase map. Structures on the surface which are not rotating with the same periodicity as the dominant period used for folding will appear to shift forward or backwards in phase in each consecutive cycle. This leads to diagonal streaks in the phase map. For particularly stable structures, it is possible to fit a linear model to calculate the rotational period (Davenport et al., 2015, 2020).

Figure 12 shows the phase map for TIC 364588501. During the identified block between cycles 57 and 69, the flare phase distribution histogram suggests a possible dependence on phase. However, in the top panel for Figure 13, we show an excerpt of the light curve during this block. It can be seen that unlike TIC 382575967, the oscillations in the light curve do not match the folding period, resulting in diagonal streaks. Visual analysis of the light curve shows a complicated oscillation pattern up to cycle 62 which appears to be two out of phase sinusoids with different amplitudes. By the end of cycle 62, one of these sinusoidal oscillations has decayed, and starting with cycle 64 we see a clear single sinusoidal oscillation. We interpreted this as two groups of starspot structures, with one decaying by cycle 62.

In order to measure the phase for each flare with respect to the observed oscillations, we first identify the start and end of each cycle. In Figure 13, we show an excerpt of the light curve. Each cycle is identified by locating a minimum in the quiescent luminosity model, and then stepping forward by the measured period and identifying the minimum in a one day window around this time. For cycles at the start and end of a period of continuous observation, we assume their length is the measured period. The length of each cycle identified in this way varies, with the cycles having an average duration of $2.38\pm 0.55$ days. Calculating the phase of each data point and flare using these identified cycles gives the phase map shown in Figure 13. We calculate the phase distribution separately for the flares before and after the fifth manually identified cycle, corresponding to the aforementioned change in the oscillation pattern. Individually, each phase distribution appears to be roughly uniform, with the phase distribution of the starspot structure after the fifth cycle (shown in orange) showing a possible preference to flaring towards the peak of the oscillation. We individually compared each of the distributions to a uniform distribution using a K-S test. For both distributions, the $p$-values from the K-S test were high ($p$-values of 0.93, 0.53 for before and after the fifth cycle, respectively), suggesting that both of these phase distributions can be fit well with a uniform distribution.