Detecting stellar flares in the presence of a deterministic trend and stochastic volatility

TL;DR

The paper introduces a two-stage framework to detect stellar flares in time series with deterministic baseline trends and stochastic volatility. It first removes a time-varying harmonic baseline, then models the residuals with an ARMA-GARCH process, identifying flares as large deviations that cannot be explained by the established correlation structure, with rigorous control of false discoveries via Holm-Bonferroni and Benjamini-Hochberg procedures. Applied to three TESS stars, the method detects hundreds of flares across sectors and cadences, and finds power-law distributions for flare energies and peak fluxes with indices that vary by stellar type. The approach enhances flare detectability in challenging baselines, enabling robust population studies and insights into flare onset physics and coronal heating mechanisms.

Abstract

We develop a new and powerful method to analyze time series to rigorously detect flares in the presence of an irregularly oscillatory baseline, and apply it to stellar light curves observed with TESS. First, we remove the underlying non-stochastic trend using a time-varying amplitude harmonic model. We then model the stochastic component of the light curves in a manner analogous to financial time series, as an ARMA+GARCH process, allowing us to detect and characterize impulsive flares as large deviations inconsistent with the correlation structure in the light curve. We apply the method to exemplar light curves from TIC13955147 (a G5V eruptive variable), TIC269797536 (an M4 high-proper motion star), and TIC441420236 (AU Mic, an active dMe flare star), detecting up to $145$, $460$, and $403$ flares respectively, at rates ranging from ${\approx}0.4$--$8.5$~day$^{-1}$ over different sectors and under different detection thresholds. We detect flares down to amplitudes of $0.03$%, $0.29$%, and $0.007$% of the bolometric luminosity for each star respectively. We model the distributions of flare energies and peak fluxes as power-laws, and find that the solar-like star exhibits values similar to that on the Sun ($α_{E,P}\approx1.85,2.36$), while for the less- and highly-active low-mass stars $α_{E,P}>2$ and $<2$ respectively.

Figures (14)

• Figure 1: Illustration of the similarity between financial returns and astronomical time series. The time series of the log returns of the Dow Jones Utility Average (sampled daily) during the years 2021 and 2022 (DJUA; top left) bears remarkable similarity to the TESS light curve (at a 10 minute cadence) of the star TIC 13955147-s0031 after removing the harmonic baseline ( top right). Both time series have $\approx$500 points. The fitted conditional volatility for both time series are shown in the corresponding lower panels, with the DJUA ( bottom left) showing variations that are similar to those seen for the star ( bottom right). This similarity motivates our adoption of techniques developed for financial data analysis to TESS light curves.
• Figure 2: Illustrating the detrending of the baseline flux with a harmonic model. Top panels: The first halves of TESS light curves of TIC 13955147 sector 31 ( left) and TIC 269797536 sector 29 ( right) are shown, with the harmonic detrended curve overplotted as a dashed curve. Bottom panels: The residuals after harmonic detrending are shown for the corresponding datasets.
• Figure 3: Illustrating the analysis of the residuals after harmonic detrending for the light curves in Figure \ref{['fig:compar']}, with the detected flares passed BH procedure marked as solid dots and passed HB procedure as hover squares. Top panels: Residuals $Z_t$ after detrending with the harmonic fit and bypassing the ARMA-GARCH filter. Bottom panels: Standardized residuals $\hat{\varepsilon}_t=\hat{Z}_t / \hat{\sigma}_t$, where $\hat{Z}_t$ is the ARMA($r,s$) residual and $\hat{\sigma}_t$ is the estimated conditional volatility.
• Figure 4: Distributions of the standardized residuals $\varepsilon_t$ for the data sets shown in Figure \ref{['fig:quatile_detect_demo']}. The vertical solid line marks the empirical 95% upper quantile of $\varepsilon_t$, the line dashed vertical line marks the BH threshold under $\alpha = 0.05$ while the dot dashed line marks the HB threshold, and the dashed curve is the fitted normal distribution shown for comparison. Time bins which exceed the 95% quantile are chosen as candidate flares.
• Figure 5: Illustration of screening the flare detections by correcting for multiple-testing. The two panels represent the candidate flares identified at the 95% threshold of the distribution in $\varepsilon_t$ in Figure \ref{['fig:epsilont_dist']} for TIC 13955147-s0031 ( left) and TIC 269797536-s0029 ( right). The $p$-values of each candidate flare is shown as a function of its rank order, with $p$ bounded below at $10^{-4}$ for the sake of visibility. The curved lines represent the estimated threshold at each rank, and the horizontal lines mark the adopted critical $p$-value thresholds. The dashed and dotted lines represent the BH (Benjamini-Hochberg) and the HB (Holm's Bonferroni) procedures respectively. Flare candidates that are found via BH and HB (i.e., rejected as arising from the standard distribution of $\varepsilon_t$) are marked with a filled circle and open square respectively, and candidates discarded as flares are plotted as open circles.