# Analysis plan: how often do Sun-like stars produce superflares? A test with TESS of the Kepler rate

Written 2026-10-24, before any light curve of the confirmatory set was read. The code in `code/` and the environment in `env/` beside this file are the search and the analysis this plan runs, frozen.

## The question

Vasilyev et al. (2024, *Science* 386, 1301) searched Kepler's long-cadence light curves of 56,450 stars with Sun-like temperatures (5000–6500 K), absolute G magnitudes of 4–6, RUWE below 1.4, and no known rotation period under 20 days, and found superflares of bolometric energy 10^34–10^35 erg at (8.63 ± 0.20) × 10^-3 per star per year and above 10^35 erg at (3.90 ± 0.60) × 10^-4 per star per year: about one above 10^34 erg per century. Earlier Kepler studies of Sun-like stars found such flares once per 2000–3000 years (Notsu et al. 2019; Okamoto et al. 2021). The two differ by a factor of about 20, and the difference matters: a superflare on the Sun of 10^34 erg would release about twenty times the energy of the Carrington flare of 1859.

This study measures the rate with a different telescope, camera, cadence, and set of stars: TESS's 2-minute light curves of nearby stars selected as Vasilyev et al. selected theirs.

## Hypotheses

The ratio ρ is the rate of superflares at or above 10^34 erg on the confirmatory stars divided by the rate Vasilyev et al.'s Table 1 implies.

- H_V: ρ = 1, the once-per-century rate.
- H_E: ρ = 1/22.5 = 0.044, the earlier rate (one per 2500 years above 10^34 erg), with the same energy distribution.

## Data

- **Light curves:** every TESS 2-minute SPOC light curve (PDCSAP flux) of sectors 1–107 that MAST's bulk download scripts list, read from MAST's public AWS bucket. `data/sample_lcs.csv` lists each file.
- **Stars** (`code/select_sample.py`, writing `data/sample_stars.csv`): TESS 2-minute targets with a median Tmag below 10.5 in TESS's target lists; TIC 8.2 effective temperature 5000–6500 K and a radius; Gaia DR3 parallax over error above 10; absolute G magnitude 4–6, from Gaia DR3 G and parallax less A_G = 2.74 E(B−V) from the TIC; Gaia DR3 RUWE below 1.4; not classified `ECL` by Gaia DR3; not in the TESS eclipsing binary catalog of Prša et al. (2022); no Gaia DR3 rotation period under 20 days; Gaia DR3 vbroad below 15 km/s or none. This leaves 42,166 stars with 128,042 light curves.
- **Development set:** stars whose TIC number is divisible by 20 (2,229 stars, 6,852 light curves). The search was built and its rules set on these alone, and while they were set their counted flares and dips were seen. They are reported separately and never pooled with the confirmatory set.
- **Confirmatory set:** every other selected star (39,937 stars, 121,190 light curves).

## The search (`code/sflare.py`, `code/inject.py`, `code/run_search.py`)

For each light curve, as the code does:

1. Cadences with PDCSAP flux and none of the quality flags 1, 2, 4, 8, 16, 32, 128 set (impulsive outliers, 512, are kept: SPOC flags many flare peaks so). Segments split at gaps over 0.25 days.
2. Trend: the median of each 0.05-day bin, a running median over 5 bins, interpolated; noise: the running median over 11 bins of each bin's median absolute deviation, times 1.4826. The trend is fitted twice, the second time without the cadences of the first pass's candidates.
3. Candidates: three or more consecutive cadences at 3σ or more, extended while the flux stays above 1σ; equivalent duration (ED) is the sum of the residual over the extent times 120 s.
4. Vetting: a candidate counts only if it is at least 0.02 days from a segment's edge and from a momentum dump; no other quality flag falls within 5 cadences of it; its peak is at 5σ or more; it rises from half maximum to its peak in 4 cadences or fewer and decays from the peak to half maximum in no fewer cadences than it rose; the Davenport et al. (2014) template fits it with a FWHM of 1 hour or less; the source its centroid shift implies (shift over fractional brightening) is no more than 0.5 pixel from the centroid, or the shift is under 3 standard errors; the background changes by no more than a quarter of the flux, or by under 4 standard errors; and the uncorrected (SAP) flux shows the same event, at 0.5–2 times the PDCSAP amplitude after crowding.
5. Energy: E = σ T_fl^4 π R² [∫ R(λ) B_λ(T★) λ dλ / ∫ R(λ) B_λ(T_fl) λ dλ] × ED, with T_fl = 9000 K, TESS's response R(λ), and the TIC's effective temperature T★ and radius R: the equation Vasilyev et al. used.
6. The same search on the light curve turned upside down gives "dips": how often noise and systematics pass the same rules.
7. Completeness: flares of the Davenport shape at 14 energies (10^33.0 to 10^35.6 erg in steps of 0.2 dex, six of each, in six runs at least a day apart) are injected into each light curve, the whole search is run again, and each injection is recovered if a vetted candidate's extent holds its peak, with its measured energy. FWHM = 11.5 minutes × (E / 10^34 erg)^0.39, with a log-normal scatter of 0.27 dex: the e-folding decay times of superflares on solar-type stars in Kepler's 1-minute data (Maehara et al. 2015, their Fig. 6), which match the template's FWHM. Seeds follow each file's name. In one light curve in five (by the hash of its name), injections are repeated with flares 0.5 and 2.5 times as long.

## Stars set aside after the search (rules fixed now)

- **TESS fast rotators:** any sector whose trend's Lomb–Scargle periodogram peaks between 0.2 and 12 days with power above 0.3 and semi-amplitude above 0.1 per cent (Vasilyev et al. set aside known periods under 20 days; TESS can't see longer ones).
- **Eclipsers:** any sector with two or more dips of ED 20 s or more (eclipsing binaries and transit hosts the catalogs missed).

## Counting

A flare counts if it passes vetting, its energy is 10^34 erg or more, and fewer than two other stars of the sample in the same sector and camera have a candidate of either sign within 3 cadences of its peak.

## Primary analysis (`code/analyze.py`)

- N_obs: counted flares in the confirmatory set.
- N_V: the number Vasilyev et al.'s rate predicts. Their two Table 1 numbers fix a cumulative power law, F(>E) = 9.02 × 10^-3 (E/10^34 erg)^-1.364 per star per year. Each 0.2-dex energy bin's share of it is multiplied by the exposure at that energy: the sum over light curves of their years of good data times the fraction of injections at that energy recovered with measured energy of 10^34 erg or more. The top bin takes the whole tail above it.
- N_E: the same for H_E.
- ρ̂ = N_obs / N_V, with a 95% interval from the exact Poisson interval on N_obs over N_V; the rate above 10^34 erg is ρ̂ × 9.02 × 10^-3 per year, and the probabilities of N_obs or fewer under H_V and of N_obs or more under H_E.

Reading it, fixed now:
- If the interval holds 1, the TESS data are consistent with the once-per-century rate; if its upper end is below 1, they show a lower rate, and by how much.
- If its lower end is above 0.044, they rule out the earlier rate; if the interval holds 0.044, they are consistent with it.
- If both hold, the test failed to separate them, and the paper says so.

## Secondary analyses, all reported whatever they show

1. Stars of 5500–6000 K whose median TESS Rvar is below 0.3 per cent, against Vasilyev et al.'s Sun-like row (9.31 × 10^-3 and 3.59 × 10^-4).
2. Gaia DR3 activity index below, and at or above, the confirmatory stars' median.
3. Tmag below 8.5, 8.5–9.5, and 9.5–10.5: ρ̂ should not depend on brightness if counted flares are on their targets and the completeness is right.
4. Thresholds of 10^34.5 and 10^35 erg.
5. Background-subtracted: N_obs less the dips counted by the same rules.
6. Completeness with flares 0.5 and 2.5 times as long, from the light curves injected so.
7. Energies with flare temperatures of 8000 K and 10,000 K, Vasilyev et al.'s numbers unchanged.
8. Each counted flare's target pixel file: the difference image of the flare against the cadences around it, located against the target. A recount without the flares that locate off their targets.
9. Strict single stars: without Gaia DR3 non-single-star solutions and without stars whose radial velocity varies (rv_chisq_pvalue below 0.01).

## Known differences from Kepler's sample, stated in advance

TESS can't see rotation periods over about 12 days, so stars with periods of 12–20 days, which Vasilyev et al. set aside, stay in; they flare more, which would raise ρ̂, and analysis 2 shows how much activity matters. TESS's stars are nearer (tens to hundreds of parsecs) and brighter than Kepler's (around a kiloparsec), and its 21-arcsecond pixels take in more of their neighbors. The flare durations the injections assume set the completeness; analysis 6 bounds it.

## Reporting

The study reports all of the above for the confirmatory set, and the development set's numbers separately, by 2026-12-15.