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

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.

💡 Research Summary

The paper presents a novel statistical framework for detecting stellar flares in TESS light curves that simultaneously accounts for deterministic periodic trends and stochastic volatility. Traditional flare‑detection techniques often rely on simple outlier thresholds or detrending methods that either over‑fit the baseline or ignore autocorrelation and heteroscedasticity, leading to missed low‑amplitude events and high false‑positive rates. To overcome these limitations, the authors construct a two‑stage time‑series model.

In the first stage, the observed flux Yₜ is decomposed into a deterministic component μ(t) and a stochastic residual Xₜ. The deterministic trend μ(t) is modeled as a sum of sinusoidal terms with time‑varying amplitudes (Eq. 5‑6). Up to K = 20 harmonics are fitted, allowing the model to capture complex rotational modulation caused by evolving star‑spot distributions. Crucially, the fitting procedure iteratively excludes data points that are likely flare‑contaminated based on p‑value thresholds, preventing flare power from leaking into the harmonic amplitudes or phases.

The second stage treats the residual Xₜ as an ARMA(r,s) process whose innovations Zₜ follow a GARCH(p,q) conditional variance model. This mirrors the well‑established financial econometrics approach where returns exhibit volatility clustering, heavy tails, and asymmetry. Here, Zₜ = σₜ εₜ with εₜ ~ i.i.d. N(0,1) and σ²ₜ = a₀ + ∑aᵢ Z²_{t‑i} + ∑bⱼ σ²_{t‑j}. By estimating σₜ, the method captures periods of heightened stochastic variability that often precede or accompany flares.

Flare detection is performed on the standardized residuals ε̂ₜ = Ẑₜ/σ̂ₜ. Points where ε̂ₜ exceeds a user‑defined significance threshold (e.g., 5σ) are flagged as positive outliers; contiguous outliers are merged into single flare events, and start/end times are defined via a secondary, lower‑threshold scan. For each flare, peak flux, duration, and equivalent duration are measured, enabling the construction of cumulative distributions of flare energy (E) and peak flux (P). These distributions are fitted with power‑law forms N(>E) ∝ E^{−α_E+1} and N(>P) ∝ P^{−α_P+1}.

The methodology is applied to three representative stars observed by TESS:

1. TIC 13955147 (G5V eruptive variable) – 145 flares detected, flare rate ≈ 0.4–1.2 day⁻¹ across multiple sectors. The power‑law indices α_E ≈ 1.85 and α_P ≈ 2.36 are close to solar values, indicating a solar‑like flare energy distribution.

2. TIC 269797536 (M4 high‑proper‑motion dwarf) – 460 flares detected, flare rate ≈ 2–8.5 day⁻¹. Both α_E and α_P exceed 2, implying a steeper distribution with relatively fewer low‑energy flares.

3. AU Mic (TIC 441420236, M1Ve active dwarf) – 403 flares detected, flare rate ≈ 3–7 day⁻¹. Here α_E ≈ 1.6 and α_P ≈ 1.8 (both < 2), reflecting the classic “over‑abundance” of high‑energy flares in highly active stars.

Detection thresholds reach remarkably low amplitudes: 0.03 % of bolometric luminosity for TIC 13955147, 0.29 % for TIC 269797536, and 0.007 % for AU Mic, demonstrating the method’s sensitivity to faint events.

Extensive simulations validate the approach, showing higher true‑positive rates and lower false‑positive rates than sigma‑clipping, Savitzky‑Golay, or Gaussian‑process detrending, especially in low signal‑to‑noise regimes. The authors discuss the physical interpretation of the conditional volatility σₜ, noting that spikes in σₜ often coincide with flare onset, suggesting a potential precursor diagnostic.

The paper concludes by highlighting several implications: (i) the successful transfer of financial econometric tools to astrophysical time‑series analysis; (ii) the ability to construct homogeneous flare catalogs across the vast TESS archive; (iii) prospects for real‑time flare forecasting using σₜ as an early‑warning indicator; and (iv) extensions to multi‑wavelength data (e.g., X‑ray, UV) where volatility clustering may also be present. Overall, this work provides a robust, statistically grounded pipeline that markedly improves flare detection fidelity and enables more accurate characterization of stellar magnetic activity.