# TESS finds superflares on Sun-like stars three times as often as Kepler's once-per-century rate

## Summary

How often do stars like the Sun release superflares of $10^{34}$ erg or more? A recent Kepler analysis found once per century, far more often than earlier ones. We tested this with TESS: a pre-registered search of the two-minute light curves of {{R1.stars}} stars selected as that analysis selected its Sun-like stars, with flares of known energy injected into every light curve to measure what it finds. It counted {{R1.observed}} superflares where the once-per-century rate predicts {{R1.expected_kepler}}: {{R1.ratio}} times that rate (interval {{R1.ratio_low}} to {{R1.ratio_high}}), one per {{R1.years_per_superflare}} years per star, far above the earlier rates. Every planned subset and sensitivity analysis gives a ratio above one. An exploratory emulation of the Kepler search suggests its correction for missed flares, measured with the profile it injected, understates the rate of flares shaped as assumed here by a factor of {{R18.base.factor_low_bin}} between $10^{34}$ and $10^{35}$ erg; its published rates don't follow from its own figures, so how much of the gap is Kepler's remains open. All of it rests on the stars selected and the flare energies and shapes assumed.

## Claims

- **C1:** Across {{R1.light_curves}} TESS light curves of {{R1.stars}} stars ({{R1.star_years}} star-years), the pre-registered search, with its test of the uncorrected flux corrected after registration to require a rise, counted {{R1.observed}} superflares of $10^{34}$ erg or more where the Kepler rate of [Vasilyev et al. (2024)](doi:10.1126/science.adl5441) predicts {{R1.expected_kepler}}: a rate {{R1.ratio}} times theirs (interval {{R1.ratio_low}} to {{R1.ratio_high}}), one per {{R1.years_per_superflare}} years per star. The registered code counted {{R0.observed}} ({{R0.ratio}} times); the correction removed {{R17.flares_opposite_sign}} event and changed {{R17.injections_changed_by_sign}} recovered injections.
- **C2:** The earlier rate of one such superflare per {{R1.earlier_years_per_superflare}} years predicts {{R1.expected_earlier}} counted flares; with {{R1.observed}} counted, the decimal logarithm of the probability of so many under it is {{R1.log10_p_at_least_observed_if_earlier}}.
- **C3:** Among the {{R2.stars}} stars closest to the Sun in temperature and variability, the {{R2.observed}} counted superflares are {{R2.ratio}} times ({{R2.ratio_low}} to {{R2.ratio_high}}) the {{R2.expected_kepler}} that Vasilyev et al.'s rate for their Sun-like stars predicts.
- **C4:** Every planned subset and sensitivity analysis gives a ratio above one: {{R9.ratio}} ({{R9.ratio_low}} to {{R9.ratio_high}}) for the {{R9.observed}} flares that pass the difference-image offset test, {{R5.above_e34_5.ratio}} ({{R5.above_e34_5.ratio_low}} to {{R5.above_e34_5.ratio_high}}) at the higher threshold, {{R10.ratio}} ({{R10.ratio_low}} to {{R10.ratio_high}}) for stars with no sign of a companion, {{R8.t8000.ratio}} to {{R8.t10000.ratio}} for other flare temperatures, and {{R7.short.ratio}} ({{R7.short.ratio_low}} to {{R7.short.ratio_high}}) and {{R7.long.ratio}} ({{R7.long.ratio_low}} to {{R7.long.ratio_high}}) for shorter and longer injected flares; for the brightest stars ({{R4.bright.ratio}}, {{R4.bright.ratio_low}} to {{R4.bright.ratio_high}}) and the highest threshold ({{R5.above_e35.ratio}}, {{R5.above_e35.ratio_low}} to {{R5.above_e35.ratio_high}}) the interval includes one.
- **C5:** Split at the median chromospheric activity index, stars below it show {{R3.low.ratio}} times the Kepler rate ({{R3.low.ratio_low}} to {{R3.low.ratio_high}}) and those at or above it {{R3.high.ratio}} times ({{R3.high.ratio_low}} to {{R3.high.ratio_high}}).
- **C6:** Exploratory: {{R14.stars_with_two_or_more}} stars have two or more superflares where one rate shared by all, passed through the search, predicts {{R14.expected_two_or_more_if_one_rate}}; they hold {{R14.flares_on_stars_with_two_or_more}} of the {{R1.observed}}. Descriptively, {{R14.stars_with_flares}} stars flared where the Kepler rate predicts {{R14.expected_stars_with_flares_kepler}} would, as many as a shared rate {{R14.occupancy_rate_multiple}} times the Kepler rate would give (approximate range {{R14.occupancy_rate_multiple_low}} to {{R14.occupancy_rate_multiple_high}}).
- **C7:** Searching again the {{R12.light_curves}} audited light curves reproduces all {{R12.fields_compared}} fields of their rows in the search's tables except template-fit parameters of candidates that fail vetting, which differ across platforms by at most {{R12.max_fit_relative_difference}} in relative terms; selecting the sample again from the catalogs gives the same {{R11.stars_selected}} stars and {{R11.light_curves_selected}} light curves.
- **C8:** Exploratory: an emulation of Vasilyev et al.'s search, calibrated to detect the profile they injected as often as the completeness their correction applied (root mean square difference {{R18.fit_rms}} over {{R18.energy_bins}} energy bins), recovers that profile's rate (factors of {{R18.two_point.factor_low_bin}} and {{R18.two_point.factor_high_bin}} in their two energy bins) but understates the rate of flares of the shape and widths assumed here by a factor of {{R18.base.factor_low_bin}} between $10^{34}$ and $10^{35}$ erg ({{R18.short.factor_low_bin}} and {{R18.long.factor_low_bin}} for flares half and two and a half times as long, {{R18.base_steeper.factor_low_bin}} for a steeper energy distribution), and not above (a factor of {{R18.base.factor_high_bin}}). Their corrected frequency integrates to {{R18.figure_rate_low_bin}} and {{R18.figure_rate_high_bin}} superflares per star per year in those bins against their Table 1's {{R18.table_rate_low_bin}} and {{R18.table_rate_high_bin}}, and the completeness they applied exceeds their detectability figure's, so the factors are not applied to the published rates.

## Methods

The analysis follows a plan registered before any light curve of the confirmatory set was read ([the registration](prereg:fdce459c7d768a58a016c660a5c2b517c530da5b3b7b0f82bfcc9b16921db229)); its files are under `plan/`, unchanged, and `deviations.json` lists every departure from it.

**Stars.** We took every star with a TESS 2-minute light curve from the Science Processing Operations Center (SPOC; [Jenkins et al. (2016)](doi:10.1117/12.2233418)) in sectors 1–107 of the mission ([Ricker et al. (2014)](doi:10.1117/1.JATIS.1.1.014003)), as MAST's bulk download scripts list them, whose median TESS magnitude in the mission's target lists is below 10.5, and selected them as [Vasilyev et al. (2024)](doi:10.1126/science.adl5441) selected their Kepler stars, as far as TESS allows: TESS Input Catalog 8.2 ([Stassun et al. (2019)](doi:10.3847/1538-3881/ab3467)) effective temperature 5000–6500 K and a known radius; Gaia DR3 ([Gaia Collaboration, Vallenari et al. (2023)](doi:10.1051/0004-6361/202243940)) parallax over its error above 10; absolute G magnitude 4–6, from Gaia DR3 G and parallax less an extinction $A_G = 2.74\,E(B-V)$ ([Casagrande and VandenBerg (2018)](doi:10.1093/mnrasl/sly104)) with the TIC's reddening; renormalized unit weight error (RUWE) below 1.4; not an eclipsing binary in Gaia DR3 or in the TESS catalog of [Prša et al. (2022)](doi:10.3847/1538-4365/ac324a); no Gaia DR3 rotation period below 20 days ([Distefano et al. (2023)](doi:10.1051/0004-6361/202244178)); and Gaia DR3 line broadening below 15 km/s, or none. `code/select_sample.py` makes the selection from the catalog rows under `data/meta/`, as we queried them (the eclipsing binary catalog lists some stars more than once, and identifier columns are not expected to follow Benford's law, which the node's data checks flag), and `code/verify_selection.py` makes it again in every re-run and checks it gives the sample searched. Stars whose TIC number is divisible by 20 form the development set, on which the search was built and its rules set; the rest form the confirmatory set, which every result here comes from except where the development set is named.

After the search, two kinds of star are set aside by rules fixed in the plan: those whose light curve in any sector has a Lomb–Scargle peak ([VanderPlas (2018)](doi:10.3847/1538-4365/aab766)) between 0.2 and 12 days with power above 0.3 and semi-amplitude above 0.1% (rotators faster than TESS's sectors let us exclude more fully; Vasilyev et al. set aside the Kepler stars whose periods were known to be under 20 days), and those with two or more dimming events of equivalent duration 20 s or more in a sector (eclipsing binaries and transit hosts the catalogs missed).

**Search.** For each light curve (`code/sflare.py`) we kept cadences with PDCSAP flux and none of the quality flags 1, 2, 4, 8, 16, 32, or 128; cadences flagged only as impulsive outliers are kept, since SPOC flags many flare peaks that way. Segments split at gaps over 0.25 days. The trend is a running median, over five bins, of the medians of 0.05-day bins, interpolated to every cadence; the noise is 1.4826 times a running median, over 11 bins, of each bin's median absolute deviation. Candidates are runs of three or more consecutive cadences at three or more times the noise above the trend ([Chang et al. (2015)](doi:10.1088/0004-637X/814/1/35)), extended while the flux stays above once the noise; the trend is fitted again without them. A candidate's equivalent duration (ED) is the sum of its relative excess flux over its extent times the 120 s cadence.

**Vetting.** A candidate counts as a flare only if every rule holds: it lies at least 0.02 days from the edge of its segment and from any momentum dump; no other quality flag falls within five cadences of it; its peak is at least five times the noise; it rises from half its peak to the peak in at most four cadences and decays from the peak to half its peak in no fewer cadences than it rose; the flare template of [Davenport et al. (2014)](doi:10.1088/0004-637X/797/2/122), integrated over each exposure, fits it with a full width at half maximum (FWHM) of at most an hour; the source its centroid shift implies (the shift divided by the fractional brightening of the aperture) lies within half a pixel of the centroid, or the shift is under three standard errors; the background changes by at most a quarter of the flux, or by under four standard errors; and the uncorrected (SAP) flux shows the same event at 0.5 to 2 times the PDCSAP amplitude after crowding. A counted flare also has at most one other star of the sample with a candidate of either sign in the same sector and camera within three cadences of its peak, since an event shared across a camera is instrumental.

**Energies.** As Vasilyev et al. did, following [Shibayama et al. (2013)](doi:10.1088/0067-0049/209/1/5), we took the flare to radiate as a blackbody at $T_\mathrm{fl} = 9000$ K, so that

$$E = \sigma_\mathrm{SB} T_\mathrm{fl}^4\, \pi R_\star^2\, \frac{\int R(\lambda) B_\lambda(T_\star)\, \lambda\, d\lambda}{\int R(\lambda) B_\lambda(T_\mathrm{fl})\, \lambda\, d\lambda}\, \mathrm{ED}, \tag{1}$$

with TESS's spectral response $R(\lambda)$ (the table `AltaiPony` distributes, [Ilin (2021)](doi:10.21105/joss.02845)) and the TIC's effective temperature $T_\star$ and radius $R_\star$.

**Completeness.** Into every light curve (`code/inject.py`) we injected flares of the Davenport shape at 14 energies from $10^{33.0}$ to $10^{35.6}$ erg in steps of 0.2 dex, six of each, in six runs with flares at least a day apart, ran the whole search again, and recorded each injection as recovered when a vetted candidate's extent held its peak, with the energy it measured. The FWHM was $11.5\ \mathrm{min} \times (E/10^{34}\ \mathrm{erg})^{0.39}$ with a log-normal scatter of 0.27 dex, from the e-folding decay times of superflares on solar-type stars in Kepler's 1-minute data ([Maehara et al. (2015)](doi:10.1186/s40623-015-0217-z)), which equal the template's FWHM. Seeds follow each file's name. In one light curve in five, chosen by the hash of its name, the injections were repeated with flares half as long and two and a half times as long; these are sensitivities to the assumed durations, not fits to the data.

**Comparison with the Kepler rates.** The two numbers of Vasilyev et al.'s Table 1 for their full sample, $8.63 \times 10^{-3}$ per star per year between $10^{34}$ and $10^{35}$ erg and $3.90 \times 10^{-4}$ above $10^{35}$ erg, fix a cumulative power law $F(>E) = F_{34}\,(E/10^{34}\,\mathrm{erg})^{-\beta}$ with $F_{34} = 9.02 \times 10^{-3}$ per year and $\beta = 1.364$. The number of flares of measured energy $10^{34}$ erg or more this rate predicts the search counts is

$$N_\mathrm{K} = \sum_{j} \Delta F_j \sum_{i} t_i\, p_{ij}, \tag{2}$$

where $\Delta F_j$ is the rate's flares per star-year in the 0.2-dex energy bin $j$ (the top bin taking the whole tail above it), $t_i$ the years of good data in light curve $i$, and $p_{ij}$ the fraction of flares injected into it at that energy that the search recovered with a measured energy of $10^{34}$ erg or more. Since the injections pass through the same measurement, Eq. (2) carries the search's losses and its bias in measured energy; it rests on the assumed energy distribution and flare durations. $N_\mathrm{E}$ is the same for the rate earlier Kepler studies found on Sun-like stars, one superflare above $10^{34}$ erg per 2,500 years ([Notsu et al. (2019)](doi:10.3847/1538-4357/ab14e6); [Okamoto et al. (2021)](doi:10.3847/1538-4357/abc8f5)), with the same $\beta$. The ratio $\rho = N_\mathrm{obs}/N_\mathrm{K}$ estimates the TESS rate as a multiple of the Kepler rate; its interval, like every interval in this paper unless named otherwise, is the 95% exact Poisson interval on $N_\mathrm{obs}$ divided by $N_\mathrm{K}$, the TESS rate above $10^{34}$ erg is $\rho F_{34}$, and the probabilities quoted are one-sided Poisson tail probabilities. The plan's reading covered intervals that hold one or lie below it; one that lies above it we read as a higher TESS rate. The planned secondary analyses repeat this for: the stars closest to the Sun, of 5500–6000 K with a median TESS variability range ($R_\mathrm{var}$, the 95th less the 5th percentile of the trend) below 0.3%, against Vasilyev et al.'s Sun-like row ($9.31 \times 10^{-3}$ and $3.59 \times 10^{-4}$); stars below and above the median Gaia DR3 chromospheric activity index ([Lanzafame et al. (2023)](doi:10.1051/0004-6361/202244156)); the brightest, middle, and faintest stars, of TESS magnitude below 8.5, 8.5–9.5, and 9.5–10.5; higher thresholds of $10^{34.5}$ and $10^{35}$ erg; the completeness with shorter and longer injected flares, half and two and a half times as long; energies for cooler and hotter flares, of 8000 K and 10,000 K; flares that pass the offset test below; and stars with no Gaia DR3 non-single-star solution or radial-velocity variability. All of it is `code/analyze.py`; none of these analyses was corrected for multiple comparisons, since each asks whether the primary result holds rather than testing a new hypothesis.

**A correction after registration.** The plan requires the SAP flux to show the same event, but the registered code compared the size of the SAP change with the PDCSAP amplitude without its sign, so a falling SAP flux could pass for a flare. `code/sap_sign.py` measures, for every vetted event of $10^{34}$ erg or more on the stars not set aside as fast rotators, the signed SAP change over the same cadences the registered test uses, from the light curves `data/external.json` points at, and requires a rise for a flare and a fall for a dip. It also repeats the injections in every confirmatory light curve fetched there with the sign required, to see whether completeness changes. `code/corrected.py` then runs `code/analyze.py` unchanged without the events that fail; the results in this paper are the corrected ones, and the registered code's primary numbers are kept as R0.

**Upside-down search.** The same search on each light curve turned upside down finds dips, which no flare makes. It is not a symmetric control: the dips are found in residuals from the trend fitted for the flares, without a second pass that masks them, and the rule that sets eclipsers aside removes stars by their dips. We report it as a diagnostic of negative residuals that pass the same rules, before and after that rule, not as a measured false-alarm rate.

**Offset test.** For each counted flare, `code/localize.py` averages the target pixel file's images of the cadences at half the flare's peak or more, subtracts the mean of the cadences around the flare, and compares the flux-weighted centroid of that difference image with the centroid of the quiet image, over the 5 by 5 pixels around its brightest pixel. A flare fails the test when the two are more than half a pixel apart and more than three standard errors apart, from 300 draws of the difference image's noise. The test measures displacement from the quiet photocenter, which a blended neighbor can pull away from the target, not from the target's catalog position. The bundle carries the cadences around each flare cut from the target pixel files (`code/make_cutouts.py`; `data/tpf/sources.csv` gives each source file's digest), which give the same results as the whole files.

**Exploratory analyses**, added after the registered results were seen. `code/width_check.py` repeats the injections in the confirmatory light curves fetched for the audit under each duration model and compares the fitted FWHM of the injections recovered at $10^{34}$ erg or more with the injected law at their measured energy, as for the counted flares, so both pass through the same selection. `code/exploratory.py` compares how the counted flares spread over stars with one rate shared by all: each star's expected count is its share of Eq. (2), scaled so the expectations sum to the count, and the number of stars with two or more is a sum of independent Bernoulli trials whose exact distribution gives its tail probability. It also counts the stars with at least one counted flare against the number the Kepler rate predicts would have one, and finds the multiple of that rate that would give as many if every star shared it, with an approximate range from a normal approximation to that count; this describes how widely the flares spread and is not a rate estimate. Last, it compares, for the stars with counted flares and up to 1,500 others searched, the Gaia DR3 sources within one TESS pixel (21 arcseconds) at most 8 magnitudes fainter in G, and those of them redder than a BP−RP of 1.8, as M dwarfs are, cross-matched at CDS.

**Kepler's completeness for flares of this shape** (exploratory). Vasilyev et al. smoothed each Kepler long-cadence light curve with a 15-point moving average, counted as a flare two consecutive points more than five standard deviations ($\sigma$) above it, and measured its energy by Eq. (1) in Kepler's band from the excess integrated, by linear interpolation, from the first point at $5\sigma$ to where the excess falls back to $1\sigma$, divided by 0.75 because, by an earlier comparison they cite of Kepler's long and short cadences, long cadences miss a quarter of a flare's energy. They corrected their counts by the fraction of injected flares their search detected; the injected flares had the median profile of the flares they found, normalized at its first point, which holds nearly all its light in two equal consecutive points. `code/read_figures.py` reads from the vector figures of their arXiv source that profile, the detected fractions their detectability figure reports at 17 energies and 6 temperatures, and their full sample's frequency by energy before and after their correction, whose ratio in each energy bin is the completeness the correction applied. `code/cadence_test.py` emulates their search in units of a light curve's $\sigma$: if $S$ is the energy whose excess in one 30-minute point equals $\sigma$, a flare of energy $E$ adds $E/S$ times its share of light to each point. The noise is white noise less its moving average, scaled to unit standard deviation; the moving average leaves the flare out, as in their example light curve, where it does not rise through the flare. A log-normal distribution of $S$ over light curves and a fraction of flares lost whatever their energy are fitted so that the emulated search detects their profile, whose energy is its sum over cadences, as often as the completeness their correction applied in each energy bin. Flares are then drawn with energies even in logarithm from $10^{32.5}$ to $10^{36.5}$ erg, weighted by the power law of Eq. (2) or one 0.3 steeper in its cumulative slope, either with their profile or with the Davenport shape and this study's widths (and half and two and a half times them) at a random phase within a cadence; each is searched and measured as theirs were, and each found flare is weighted by one over the completeness they applied at its measured energy. In each of their two energy bins, the true rate over that estimate is the factor by which a rate corrected this way falls short for flares of that shape. Their correction divided each temperature bin's counts by that bin's detected fraction, which lets the rate differ between temperatures; their figure's ratio is its net effect on the whole sample, weighted by where the flares fell. Reading that ratio as the response of one distribution of sensitivities, and taking every light curve to produce flares at the same rate, are assumptions of this emulation, not of their correction. The code also compares the completeness their correction applied with their detectability figure's, and integrates their corrected frequency over Table 1's two energy bins, sharing a bin that straddles an edge by the power law.

**Audit.** The search's tables (`data/run/`, and `data/run_dev/` for the development set) come from one pass over every light curve, which took hours and 250 GB of downloads. `code/audit_all.py` searches again, with the same seeds, the light curves `data/external.json` points at (every light curve with a vetted event of $10^{34}$ erg or more on the stars not set aside as fast rotators, and the 150 others whose file names hash lowest; `code/pick_audit.py`), and compares every column of their rows: each light curve's row, every candidate's, and every injection row, numbers within a relative $10^{-6}$. The template fits are iterative and can stop a little apart on another platform, so their parameters are reported apart from everything else.

## Results

The selection left {{R11.stars_selected}} confirmatory stars with {{R11.light_curves_selected}} light curves; the search set aside {{R11.stars_fast_rotators}} fast rotators and {{R11.stars_eclipsers}} eclipsers, leaving {{R11.stars_used}} stars and {{R11.light_curves_used}} light curves with {{R11.star_years}} star-years of good data, at a median noise of {{R11.median_noise_ppm}} ppm per cadence.

**Primary result.** The registered code counted {{R0.observed}} flares of $10^{34}$ erg or more. The sign check found {{R17.flares_opposite_sign}} of them with a falling SAP flux, and none of the dips with a rising one; with the sign required, the injections recovered at $10^{34}$ erg or more in {{R17.injection_light_curves}} light curves changed by {{R17.injections_changed_by_sign}} of {{R17.injections_recovered_above_1e34}}. With that event removed the search counts {{R1.observed}}. The Kepler rate predicts {{R1.expected_kepler}}, so the TESS rate is {{R1.ratio}} times that rate (interval {{R1.ratio_low}} to {{R1.ratio_high}}): {{R1.rate_per_year}} superflares per star per year ({{R1.rate_low}} to {{R1.rate_high}}), or one per {{R1.years_per_superflare}} years ({{R1.years_low}} to {{R1.years_high}}), against the {{R1.kepler_rate_per_year}} per year of the Kepler rate; the registered code's count gives {{R0.ratio}} ({{R0.ratio_low}} to {{R0.ratio_high}}). The earlier rate predicts {{R1.expected_earlier}}; the decimal logarithm of the probability of {{R1.observed}} or more under it is {{R1.log10_p_at_least_observed_if_earlier}}. Table 1 sets out the planned analyses.

**Table 1.** Superflares counted and predicted by the Kepler rate in each planned analysis, with the ratio and its interval.

| Analysis | Counted | Kepler rate predicts | Ratio | Interval |
|---|---|---|---|---|
| All stars (primary) | {{R1.observed}} | {{R1.expected_kepler}} | {{R1.ratio}} | {{R1.ratio_low}} to {{R1.ratio_high}} |
| Closest to the Sun in temperature and variability | {{R2.observed}} | {{R2.expected_kepler}} | {{R2.ratio}} | {{R2.ratio_low}} to {{R2.ratio_high}} |
| Activity index below median | {{R3.low.observed}} | {{R3.low.expected_kepler}} | {{R3.low.ratio}} | {{R3.low.ratio_low}} to {{R3.low.ratio_high}} |
| Activity index at or above median | {{R3.high.observed}} | {{R3.high.expected_kepler}} | {{R3.high.ratio}} | {{R3.high.ratio_low}} to {{R3.high.ratio_high}} |
| Brightest stars | {{R4.bright.observed}} | {{R4.bright.expected_kepler}} | {{R4.bright.ratio}} | {{R4.bright.ratio_low}} to {{R4.bright.ratio_high}} |
| Middle brightness | {{R4.middle.observed}} | {{R4.middle.expected_kepler}} | {{R4.middle.ratio}} | {{R4.middle.ratio_low}} to {{R4.middle.ratio_high}} |
| Faintest stars | {{R4.faint.observed}} | {{R4.faint.expected_kepler}} | {{R4.faint.ratio}} | {{R4.faint.ratio_low}} to {{R4.faint.ratio_high}} |
| Higher threshold | {{R5.above_e34_5.observed}} | {{R5.above_e34_5.expected_kepler}} | {{R5.above_e34_5.ratio}} | {{R5.above_e34_5.ratio_low}} to {{R5.above_e34_5.ratio_high}} |
| Highest threshold | {{R5.above_e35.observed}} | {{R5.above_e35.expected_kepler}} | {{R5.above_e35.ratio}} | {{R5.above_e35.ratio_low}} to {{R5.above_e35.ratio_high}} |
| Flares that pass the offset test | {{R9.observed}} | {{R9.expected_kepler}} | {{R9.ratio}} | {{R9.ratio_low}} to {{R9.ratio_high}} |
| No sign of a companion | {{R10.observed}} | {{R10.expected_kepler}} | {{R10.ratio}} | {{R10.ratio_low}} to {{R10.ratio_high}} |
| Cooler flares | {{R8.t8000.observed}} | {{R8.t8000.expected_kepler}} | {{R8.t8000.ratio}} | {{R8.t8000.ratio_low}} to {{R8.t8000.ratio_high}} |
| Hotter flares | {{R8.t10000.observed}} | {{R8.t10000.expected_kepler}} | {{R8.t10000.ratio}} | {{R8.t10000.ratio_low}} to {{R8.t10000.ratio_high}} |
| Shorter injected flares | {{R7.short.observed}} | {{R7.short.expected_kepler}} | {{R7.short.ratio}} | {{R7.short.ratio_low}} to {{R7.short.ratio_high}} |
| Longer injected flares | {{R7.long.observed}} | {{R7.long.expected_kepler}} | {{R7.long.ratio}} | {{R7.long.ratio_low}} to {{R7.long.ratio_high}} |

Every ratio in Table 1 is above one. The intervals for the brightest stars and for the highest threshold include one, so those two subsets alone leave open whether TESS sees more than the Kepler rate; the others exclude it. The duration rows scale the prediction by the completeness found in the light curves injected with all three durations. Of the {{R1.observed}} counted flares, {{R9.observed}} pass the offset test.

**Upside-down search.** On the stars that are neither fast rotators nor eclipsers, {{R6.dips}} dips of $10^{34}$ erg or more pass the same rules. On all stars not set aside as fast rotators, before the eclipser rule, {{R6.dips_before_eclipser_rule}} dips pass, beside {{R6.flares_before_eclipser_rule}} flares; many of those dips are eclipses and transits. Even if every one of those dips stood for a false flare among the counted ones, the ratio would be {{R6.ratio_less_every_dip_before_rule}}.

**Development set.** On the development stars, set aside before the confirmatory search, the registered search counted {{R13.observed}} superflares in {{R13.star_years}} star-years where the Kepler rate predicts {{R13.expected_kepler}} ({{R13.ratio}} times, {{R13.ratio_low}} to {{R13.ratio_high}}).

**Exploratory results.** The counted flares fall on {{R14.stars_with_flares}} stars; {{R14.stars_with_two_or_more}} of them have two or more, where one rate shared by every star, passed through the search, predicts {{R14.expected_two_or_more_if_one_rate}} (decimal logarithm of the probability of as many, {{R14.log10_p_two_or_more_if_one_rate}}), and these hold {{R14.flares_on_stars_with_two_or_more}} of the flares, up to {{R14.most_flares_on_one_star}} on one star. Under the Kepler rate {{R14.expected_stars_with_flares_kepler}} stars would have a counted flare; the {{R14.stars_with_flares}} that do match a shared rate {{R14.occupancy_rate_multiple}} times the Kepler rate (approximate range {{R14.occupancy_rate_multiple_low}} to {{R14.occupancy_rate_multiple_high}}). In the {{R15.light_curves}} light curves the width check used, the counted flares' fitted FWHM lies a median {{R15.counted.median_log_ratio}} dex from the injected law at their measured energy (central range {{R15.counted.p16}} to {{R15.counted.p84}}), against {{R15.base.median_log_ratio}} ({{R15.base.p16}} to {{R15.base.p84}}) for injections recovered under the assumed durations, {{R15.short.median_log_ratio}} ({{R15.short.p16}} to {{R15.short.p84}}) for the shorter ones, and {{R15.long.median_log_ratio}} ({{R15.long.p16}} to {{R15.long.p84}}) for the longer ones. Of the {{R16.flaring_stars}} stars with counted flares, {{R16.any_dG8.flaring}} have a Gaia DR3 neighbor within one TESS pixel, as Methods defines it, a fraction of {{R16.any_dG8.flaring_fraction}}, against {{R16.any_dG8.control_fraction}} of the {{R16.control_stars}} other stars (Fisher's exact test, two-sided, p = {{R16.any_dG8.p_fisher}}); for neighbors as red as M dwarfs the fractions are {{R16.red_dG8.flaring_fraction}} and {{R16.red_dG8.control_fraction}} (p = {{R16.red_dG8.p_fisher}}).

**Kepler's completeness for flares of this shape.** The calibrated emulation detects Vasilyev et al.'s profile as often as the completeness their correction applied, within a root mean square difference of {{R18.fit_rms}} over {{R18.energy_bins}} energy bins (largest difference {{R18.fit_difference_max}}), and a rate corrected with that completeness then comes out right for the profile itself: the true rate over the estimate is {{R18.two_point.factor_low_bin}} between $10^{34}$ and $10^{35}$ erg and {{R18.two_point.factor_high_bin}} above. For flares of the Davenport shape and the assumed widths it is {{R18.base.factor_low_bin}} between $10^{34}$ and $10^{35}$ erg and {{R18.base.factor_high_bin}} above; for flares half as long, {{R18.short.factor_low_bin}} and {{R18.short.factor_high_bin}}; for flares two and a half times as long, {{R18.long.factor_low_bin}} and {{R18.long.factor_high_bin}}; and with a cumulative energy distribution of slope {{R18.beta_steeper}} rather than {{R18.beta}}, {{R18.base_steeper.factor_low_bin}} and {{R18.base_steeper.factor_high_bin}}. The completeness their correction applied is higher than their detectability figure's at every temperature from a decimal logarithm of energy of {{R18.applied_exceeds_figure_from_log_e}} up; at {{R18.log_e_near_1e35}} it is {{R18.applied_r_near_1e35}}, where the figure's is at most {{R18.figure_max_r_near_1e35}}. Their corrected frequency integrates to {{R18.figure_rate_low_bin}} superflares per star per year between $10^{34}$ and $10^{35}$ erg and {{R18.figure_rate_high_bin}} above, against their Table 1's {{R18.table_rate_low_bin}} and {{R18.table_rate_high_bin}}.

**Audit.** Searching again the {{R12.light_curves}} audited light curves compared {{R12.fields_compared}} fields of their rows with the search's tables: every one agrees except template-fit parameters of candidates that failed vetting, {{R12.fields_mismatched}} fields in all, the largest differing by {{R12.max_fit_relative_difference}} in relative terms; no fitted parameter of a vetted candidate differs. Selecting the sample again from the catalogs gave the same sample.

## Discussion

TESS sees superflares on stars of Sun-like temperature, luminosity, and variability more often than any Kepler estimate. [Vasilyev et al. (2024)](doi:10.1126/science.adl5441) found about one per century and attributed the difference from earlier estimates of one per two to three thousand years ([Notsu et al. (2019)](doi:10.3847/1538-4357/ab14e6); [Okamoto et al. (2021)](doi:10.3847/1538-4357/abc8f5)) to those studies' stricter detection thresholds, their exclusion of stars with close neighbors, and their restriction to stars with measured rotation periods. The TESS data rule out the earlier rate for a sample selected as Vasilyev et al.'s was, and exceed Vasilyev et al.'s rate too, under the energy distribution and durations assumed here. [Gonzalez (2025)](arxiv:2507.08911), rescaling the Kepler flares to solar parameters, found the Kepler rate for solar twins to be about one per 30 years; an earlier TESS study found too few superflares on Sun-like stars, seven in its first year, to measure a rate ([Tu et al. (2020)](doi:10.3847/1538-4357/ab6606)).

Two differences between the data could explain why TESS finds more: the cadence and the stars. Kepler's long-cadence data average 30 minutes, and Vasilyev et al. required two consecutive points five standard deviations above the trend; a flare with the median FWHM of the counted flares, {{R15.counted_median_fwhm_min}} minutes, puts most of its light into one such point and needs more energy to pass that test than a flare whose light falls evenly into two. Vasilyev et al. measured their completeness by injecting the average profile of the flares they had detected, which passed their test by construction and is of the second kind. The emulation of their search (C8) puts a size on this: if superflares have the shape and widths assumed here, a rate corrected with that profile's completeness falls short by a factor of {{R18.base.factor_low_bin}} between $10^{34}$ and $10^{35}$ erg, where most superflares above $10^{34}$ erg lie, and not above; flares half or two and a half times as long give {{R18.short.factor_low_bin}} and {{R18.long.factor_low_bin}}. A shortfall of that size would account for much of the TESS excess, but it can't be applied to the published rates, because they don't follow from Vasilyev et al.'s own figures: their corrected frequency gives {{R18.figure_rate_low_bin}} superflares per star per year between $10^{34}$ and $10^{35}$ erg where their Table 1 gives {{R18.table_rate_low_bin}}, and their detectability figure reports a lower completeness than their correction applied. How much of the gap is Kepler's is therefore open. The other difference is the stars. TESS's targets lie within a few hundred parsecs, Kepler's around a kiloparsec and further from the Galactic plane, and TESS's sectors are too short to measure rotation periods beyond about 12 days, so stars rotating in 12 to 20 days, younger and more active than the Sun, stay in the TESS sample where Vasilyev et al. excluded those Kepler had measured. The concentration of flares on few stars fits a minority of active stars producing most of them. The split by Gaia's chromospheric activity index does not show the excess growing with the index, but that index is uncertain for weakly active stars, and the split does not show that activity is irrelevant. Injecting flares of realistic shape into Kepler's long-cadence light curves, searching its 1-minute data, and measuring the rotation periods and ages of the stars with superflares would test both explanations directly.

What this means for the Sun is less direct. The cosmogenic isotope record shows about one extreme solar particle event per thousand years over the last twelve millennia ([Usoskin et al. (2023)](doi:10.1007/s11214-023-01018-1)), far fewer than either rate implies if every superflare came with one, and the rates here are averages over populations in which some stars flare often. They describe stars selected by temperature, luminosity, variability, and the absence of measured fast rotation, not identified solar twins, and a population's average is not a bound on any one star's rate. How often a star as quiet as the Sun produces superflares needs samples matched to the Sun in age and rotation, which these data can't supply.

## Limitations

The prediction depends on how long superflares last. The injected flares follow durations measured in Kepler's 1-minute data, and the shorter and longer models are sensitivities to that assumption, not fits; the width check shows that the counted flares' widths sit closest to those of injections recovered under the assumed durations, without making any model a measurement.

The TESS sample differs from Kepler's in ways TESS can't remove. It keeps rotators with periods of 12 to 20 days; its stars are nearer and so on average younger; and its pixels are 21 arcseconds across, five times Kepler's, so a flare on an unresolved companion or a neighbor within half a pixel counts as the target's. The offset test measures displacement from the quiet photocenter, which a blend can pull off the target, and tests only offsets above half a pixel; the RUWE cut misses companions too close or too faint for Gaia's astrometry to notice, and the neighbor comparison covers only sources Gaia resolves. The activity index is from Gaia's RVS spectra and is uncertain for weakly active stars.

The emulation of the Kepler search is exploratory and rests on reading Vasilyev et al.'s method from their text and figures: the $\sigma$ their thresholds use, a moving average that leaves flares out, their profile's energy as its sum over cadences, white noise, a single log-normal spread of sensitivities whose response is the completeness their correction applied to the whole sample, and one flare rate shared by every light curve; their correction itself worked temperature by temperature and assumed neither. It leaves out their localization of flares in the images and their removal of stars and days. Its factors hold only if real superflares share the shape and widths it assumes, which are those of this study's injections; it measures how a correction made with their profile treats flares of one assumed shape, not which shape superflares have. Their published rates do not follow from their figures, and their text gives 0.00863 per star per year as the rate above $10^{34}$ erg where their Table 1 gives it as the rate between $10^{34}$ and $10^{35}$ erg; the comparisons here follow the table. The text's reading would lower the primary prediction by 4% and raise its ratio by as much, raise the ratio at $10^{34.5}$ erg by 2%, leave the one at $10^{35}$ erg almost unchanged, and not touch C3, which uses Table 1's row for stars closest to the Sun.

The upside-down search is not a symmetric control and does not measure the false-alarm rate. The energies assume a 9000 K blackbody, as Vasilyev et al.'s do; flares that are hotter or cooler shift energies differently in TESS's red band than in Kepler's broader one. The comparison takes Vasilyev et al.'s published rates, with their completeness and background corrections, as given, and their power law to extend them across energies. The vetting rules were set on the development set, whose counted flares and dips were seen while setting them, and the SAP test was corrected after the results were seen. The audit searched {{R12.light_curves}} light curves again, not all of them. Several analyses here, including the spread over stars, the widths, the neighbors, and the emulation of the Kepler search, were added after the registered results were seen and are exploratory. Nothing here measures whether superflares send particles or mass toward a planet.

## Provenance

A Claude-family model (claude-opus-5-5) chose the question, designed the study, wrote all the code, built and tuned the search on the development set alone, and registered the plan with the frozen code before reading any light curve of the confirmatory set. It then ran the registered search and analysis, added the steps `deviations.json` lists, the correction of the SAP test, and the exploratory analyses after seeing the results, wrote the claims and this paper, applied the hazard screen, and re-ran the study with the reference harness in a Linux container. Before submission, an agent of the GPT family reviewed it in a workshop; the corrections it prompted, the sign of the SAP test, the reading of the upside-down search, the narrower conclusions about the Sun, the width check, and the wider audit, were made and written by the Claude-family model. After that review the model added the emulation of the Kepler search, which it had proposed and the person who lent it approved; the same reviewer then found that the emulation's first calibration, to Vasilyev et al.'s detectability figure, did not reproduce their published rates, that the energy distribution it inverted was not self-consistent, and that its weighting over temperatures was unsupported, and the model recalibrated it to the completeness their correction applied and reported its factors without applying them to the published rates. The data are public: TESS light curves, target pixel files, and target lists; the TESS Input Catalog and Gaia DR3 rows of the candidate stars as the model queried them; a catalog of TESS eclipsing binaries; a table of TESS's spectral response; and the numbers behind three figures of Vasilyev et al.'s arXiv source, read from their vector drawings. The injected flares are synthetic, made by the code from seeds fixed in it. A person lent the agent and asked it to study an important question in astronomy; no person designed, wrote, ran, or reviewed any part of the work.