Flare Hunting in Hot Subdwarf and White Dwarf Stars from Cycles 1–5 of TESS Photometry

Stellar flares are critical phenomena on stellar surfaces, which are closely tied to stellar magnetism. While extensively studied in main-sequence (MS) stars, their occurrence in evolved compact stars, specifically hot subdwarfs and white dwarfs (WDs), remains scarcely explored. Based on Cycles 1–5 of TESS photometry, we conducted a pioneering survey of flare events in ∼12,000 compact stars, corresponding to ∼38,000 light curves with a 2 minute cadence. Through dedicated techniques for detrending light curves, identifying preliminary flare candidates, and validating them via machine learning, we established a catalog of 1016 flares from 193 compact stars, including 182 from 58 sdB/sdO stars and 834 from 135 WDs, respectively. However, all flaring compact stars showed signs of contamination from nearby objects or companion stars, preventing sole attribution of the detected flares. For WDs, it is highly probable that the flares originated from their cool MS companions. In contrast, the higher luminosities of sdB/sdO stars diminish companion contributions, suggesting that detected flares originated from sdB/sdO stars themselves or through close magnetic interactions with companions. Focusing on a refined sample of 23 flares from 13 sdB/sdO stars, we found their flare frequency distributions were slightly divergent from those of cool MS stars; instead, they resemble those of hot B/A-type MS stars having radiative envelopes. This similarity implies that the flares on sdB/sdO stars, if these flares did originate from them, may share underlying mechanisms with hot MS stars, which warrants further investigation.

Stellar flares are abrupt and intense phenomena that manifest as a rapid increase in luminosity across a broad wavelength coverage (Benz & Güdel 2010), posing significant impacts on the habitability of orbiting exoplanets (e.g., Vida et al. 2017).The underlying mechanisms that generate stellar flares are believed to be analogous to those of solar flares (Davenport 2016), closely related to magnetic activity on stellar surfaces.They are triggered by the sudden release of magnetic energy through the reconnection of twisted magnetic field lines within the stellar atmosphere.In general, stars with deeper convective layers are capable of generating more vigorous dynamo mechanisms (Charbonneau 2010), leading to stronger surface magnetic fields and higher frequencies of flares.This relationship manifests as an observed trend where the incidence of flare stars increases progressively from F-type to M-type main-sequence stars (Yang & Liu 2019).
Incapable of resolving details of stellar surfaces, studies of stellar flares rely on time-resolved photometric or spectroscopic observations (Walkowicz et al. 2011).Before the era of space-borne photometry, various types of observations of flares, for instance, radio waves (Bastian et al. 1988), X-rays (Schmitt et al. 1993), and optical spectroscopy (Pettersen & Hawley 1989), offer different insights to their properties but lacking in continuous data coverage for their temporal evolution.While time-resolved spectroscopy enables detailed plasma diagnostics and detection of mass ejections of single flare events (e.g.Namekata et al. 2021), continuous photometric monitoring provides a broader view, allowing for statistical analyses of flare frequencies and energies over extended periods.The Kepler, K2, and TESS missions (Koch et al. 2010;Howell et al. 2014;Ricker et al. 2015), collecting high-quality and long consecutive photometry, have fertilized the field of flare studies, exposing various properties and correlations of flares to their hosting stellar types across the HR diagram (e.g., Davenport 2016;Günther et al. 2020;Yang et al. 2023).
Despite these advances, the knowledge of evolved compact stars, specifically hot subdwarfs and white dwarfs, remains very limited.Hot subdwarfs are evolved stars with B to O spectral types (sdB/sdO stars), burning helium in the core with a thin hydrogen envelope.They are located at the blue end of the horizontal branch, thus also known as extreme horizontal branch (EHB) stars (Heber 2009).White dwarfs (WDs) are remnants of ∼ 97% stars in the Milky Way when they cease nuclear fusion, resulting in a degenerate core typically composed of carbon and oxygen and cooling down for the rest of their lives (Saumon et al. 2022).Contrary to cool MS stars, hot compact stars do not present deep convec-tive envelopes, which prevents the generation of surface magnetic fields through dynamo processes.
Nevertheless, strong magnetism has been claimed in WDs and sdB/sdO stars (Bagnulo & Landstreet 2021;Vos et al. 2021;Pelisoli et al. 2022).Now that WDs are found with magnetism whose strength is in a large range from a few kG up to 1000 MG (e.g., Landstreet & Bagnulo 2019).Their magnetism is explained by fossil mechanisms or the conservation of magnetic flux.Recent observations have shown the presence of spots on an increasing number of WDs (e.g., Kilic et al. 2015;Hoard et al. 2018;Hermes et al. 2021), implying that flares may be able to occur on these stars, as both phenomena are closely tied to magnetic fields.Moreover, spot modulation is empirically not restricted to purely convective or strongly magnetic WDs (Hermes et al. 2017), thereby opening the possibility for flare occurrences on WDs having a radiative envelope or relatively weak magnetic field.On the other hand, magnetism in sdB/sdO stars is very rare (Landstreet et al. 2012) and only a few of them were claimed to have detectable magnetic field (Vos et al. 2021;Pelisoli et al. 2022).Therefore, whether these magnetic compact stars can present flare events is an open interesting question to understand the mechanism between flare and magnetism.
However, it is a challenge to search for flare events in sdB/sdO stars and WDs since their sample is relatively small due to their low brightness and flare hunting needs intensive photometric monitoring.TESS mission provides the first opportunity for such research as it collected extensive photometry for a large fraction of sdB/sdO stars and WDs brighter than G < 16 (Ricker et al. 2015;Charpinet et al. 2019).TESS enables the investigation of various types of brightness variation in those compact stars, for instance, brightness modulation from different binary effects (see, e.g., Schaffenroth et al. 2022) and rapid variation from gravity and pressure pulsations (see, e.g., Romero et al. 2022; Baran et al. 2023).
Here we initiate the pioneering survey with an aim that tries to find flare events in sdB/sdO stars and WDs.
In line with this study, the light curves of compact stars can present complex variability, which brings difficulty in flare identification.For example, Pietras et al. (2022) excluded stars with spectral types earlier than F1 when detecting stellar flares in TESS photometry, due to their diverse and rapid brightness variations that can easily cause false positives.Moreover, previous literature has reported that sporadic outbursts in cool hydrogen-atmosphere pulsating white dwarfs (DAVs) (Bell et al. 2015(Bell et al. , 2016;;Hermes et al. 2015), which can be explained by limit cycles due to sufficiently resonant 3mode couplings (Luan & Goldreich 2018).Scaringi et al. (2022) also reported localized thermonuclear bursts in three accreting white dwarfs.Besides, self-lensing flares of WDs with compact companions can produce periodic brightening events if viewed close to edge on (Nir & Bloom 2023).These various outbursts have different characteristic profiles and durations compared to stellar flares.Careful inspection is therefore necessary to conclusively identify stellar flares amidst the complex variability in compact star light curves.
Although there are several existing methods to detect flare events automatically (e.g.Davenport 2016;Van Doorsselaere et al. 2017;Günther et al. 2020;Ilin et al. 2021), they generally rely on simple outlier detection in flattened light curves without considering event morphology.As a result, they struggle to distinguish flares from artifacts or other rapid brightening events.Yang & Liu (2019) found serious contamination of previous flare catalogs by various pulsators, rapid rotators, and transits.Other methods for detecting flares rely on machine learning algorithms (e.g., Feinstein et al. 2020), but these may perform poorly when analyzing light curves with fast and complex brightness modulations.Therefore, detrending the light curves of compact variable stars correctly and distinguishing intrinsic flares from other transient events requires dedicated techniques for flare identification and validation.
The paper is structured as follows.In Section 2, we describe the TESS photometry collected on our concerning compact stars.Section 3 provides details of our methods for preliminary identification of flare candidates, including a new approach designed to address short-term periodic variations that enhances our ability to identify flares in compact stars.We also describe several indispensable steps to exclude contamination from cataclysmic variables and solar system objects.Section 4 then explains our process to validate the preliminary candidates using a Random Forest classifier trained on a series of simulated data.We present the resulting flare catalog and analyze the properties of these events in Section 5, also considering potential contamination and establishing a refined sample of compact stars.Finally, Section 6 provides a summary of our results and a discussion of their implications regarding flare production mechanisms in compact stars.

TESS PHOTOMETRY
The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) was launched on 18 April 2018, and began observing in July 2018.After completing the primary mission in July 2020 and the first extended mission (EM1) in September 2022, TESS is currently conduct-ing its second extended mission (EM2), which will last approximately three years.
Both the primary mission and EM1 of TESS consist of two observation cycles each, while EM2 comprises three cycles (Cycles 1-2, 3-4, and 5-7, respectively).Cycles 1-3 each included 13 sectors, with Cycle 1 and 3 pointing on the southern ecliptic hemisphere and Cycle 2 on the northern hemisphere.Cycle 4, totaling 16 sectors, continued observations in the northern hemisphere and added several sectors along the ecliptic plane.The recently completed Cycle 5, consisting of 14 sectors, covered both the northern and southern ecliptic hemispheres.Each sector covers a 24 • × 96 • strip of the sky with an angular resolution of 21 arcsec pixel −1 over an observational period of ∼27 d.Due to the overlaps among the sectors in each cycle, a portion of TESS targets were observed multiple times during the survey.Targets in the vicinity of the ecliptic poles were observed continuously for almost one year.
TESS is equipped with four 10.5 cm optical telescopes and four identical cameras having a red bandpass covering the wavelength range from 600 nm to 1000 nm.Each camera has a field of view of 24 • × 24 • and consists of four 2 k×2 k back-illuminated CCDs, which continuously read out at 2-second intervals for spacecraft guiding.These 2-second data are then stacked, compressed, and stored on the spacecraft until TESS reaches its perigee after every 13.7-day orbital period, at which point TESS downlinks the collected data to Earth.
During Cycles 1-5 (Sector 1-69), TESS observed 11 618 compact stars, including 7505 WDs and 4113 sdB/sdO stars, proposed by the TESS Asteroseismic Science Consortium (TASC) 1 Working Group 8 (WG8), which focuses on variability in evolved compact stars (see, e.g., Charpinet et al. 2019;Bognár et al. 2020).There are 38 102 light curves with 2-minute cadence that have been collected for these compact stars.We retrieve their light curves from the Mikulski Archive for Space Telescopes (MAST)2 , which were processed through the standard pipeline of the Science Processing Operations Center (SPOC) (Jenkins et al. 2016).Table 1 lists the detailed observations of these compact stars: nearly twothirds of the targets were observed in one or two sectors, while about one-tenth were visited by more than five sectors.We analyze all these 38 102 light curves and characterize flare events using the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) flux, where common instrumental systematics has been  2, 3, 4, 5, 6, 8, 10, 13, and 15 are masked out, which are the recommended quality flags in the TESS Data Product documentation 3 .

FLARE CANDIDATES DETECTION
Flares are shown as consecutive positive outliers in light curves.To detect these outliers, we first detrend the light curves to remove astrophysical variability.Then, we identify groups of outliers in the detrended light curves that meet certain criteria, which are regarded as preliminary flare candidates since their profiles have not yet been considered.We exclude the light curves of known cataclysmic variables (CVs) due to their complexity and irregularity.We also inspect whether each candidate is caused by a solar system object (SSO) encounter event.

Detrending Algorithm and Window Length
For many studies on flare detection, the Savitzky-Golay filter (Savitzky & Golay 1964) has been adopted to detrend the light curve by various groups (e.g.Ilin et al. 2021).However, this filter is cadence-based and cannot function as designed on time series with observational gaps.Hippke et al. (2019) compared the relative performance of various detrending algorithms for transit discovery.They concluded that generally it is optimal to use a time-windowed slider with an iterative robust location estimator based on Tukey's biweight (Mosteller & Tukey 1977).This result can also be applied to flare detection since transits and flares are both signals con-3 https://outerspace.stsci.edu/display/TESS/2.0+-+Data+ Product+Overview stituting small segments of the light curves.Meanwhile, the Savitzky-Golay filter easily overfits the transit signals, which reduces their depth and warps up the detrended light curve around the transit.This will introduce extra outliers and bring some false alarm detection of flares (see Figure 1 & 11 in Hippke et al. 2019).Therefore, instead of using the Savitzky-Golay filter, we adopt the biweight filter implemented in Wōtan (Hippke et al. 2019), an open-source Python package for time-series smoothing, to detrend the light curves.
Next, we need to correctly determine the window length for the biweight filter before detrending.A narrow window can remove the stellar variability effectively but may overfit the large flares, which introduces underestimations of their duration and amplitude, as shown in Figure 1(b).To avoid overfitting flares, the window length should be several times longer than their durations.A recent work by Howard & MacGregor (2022, hereafter H22) identified 3792 flares from 226 M-dwarfs based on 20-s cadence photometry from TESS Cycle 3, concluding that only ∼5 % of the detected flares have a duration longer than 0.1 d.We thus set the window length to 0.3 d in the application of the biweight filter, which is considerably long enough for the durations of most flares.Such window length can also effectively filter out other transient events with typical timescales over 0.3 d in compact stars, such as pulsation outbursts in DAVs (Bell et al. 2017).

Detrending Procedure
Flares typically last from minutes to hours, rarely longer than 1 d (Ilin & Poppenhaeger 2022).In some cases, the light curves of compact stars exhibit periodic variations, on a timescale order ≲1 d, induced by factors such as pulsations or binary effects.Such variations cannot be easily removed if a wide window is chosen to preserve large flares during the detrending process, as shown in Figure 1(a).In addition, it can result in the failure to detect flares with amplitudes comparable to those periodic variations.We thus have to implement an automated procedure for detrending light curves, with a particular focus on identifying and removing these shortterm periodic variations.
We firstly normalize the light curves by dividing them by the total median and compute the Lomb-Scargle periodogram (Lomb 1976;Scargle 1982).If the period of the highest peak, P , is in the range of 0.05 d to 2 d and with significant confidence (we arbitrarily adopt 4× the median noise level), then this light curve will be considered to have short-period variability, otherwise it will be directly detrended.In some special cases, instead of the period of the maximum power signal, its harmonics rep- resent the true period of brightness variation.Therefore, the light curve will be folded with P , 2P , and 4P , respectively.We then apply a median filter with a sliding window of P fold /50 to each folded light curve to generate models of the short-period variability, where P fold is the folded period of the light curve (Figure 2).After evaluating the standard deviations of the residuals of all models, the optimum model will be selected and subtracted from the light curve.Lastly, the differential light curves will be detrended using the biweight filter with a window length of 0.3 d.Both large and small flare profiles are preserved and can be easily detected through this procedure (Figure 3).

Removing Detached Eclipses
The aforementioned procedure can effectively detrend most light curves, but it fails, in particular, for a few with detached eclipses.The folded models cannot correctly recover the periodic variations because detached eclipsing signals require many harmonics to be represented in Fourier transformation.To avoid this case, we therefore remove the eclipses before applying the biweight filter.The eclipses are masked in the light curve where more than three consecutive epochs are below −3σ MAD from the median before subtracting the folded model.Here σ MAD is a robust standard deviation calculated using the median absolute deviation (MAD), given by σ MAD ≈ 1.4826 MAD, assuming that the data is normally distributed (Huber 1981).We then extend each masked segment on both sides until an epoch is above the median of the light curve, ensuring the entire eclipse is covered (Figure 2).

Searching for Flare Candidates
After detrending the light curve, we slide a 2-day window along it and calculate its rolling σ MAD to estimate the local noise level.We then standardize the light curve by where F d is the normalized detrended light curve, σ MR is the rolling σ MAD and F s is the standardized light curve.
Standardization is a technique that transforms data to have a mean of 0 and a standard deviation of 1, commonly used when datasets have varying scales.In this case, we use standardized light curves, which represent the deviation of the flux measurements from the median in units of σ MR , to simplify identifying and validating preliminary flare candidates in subsequent analyses.
To identify groups of outliers as preliminary flare candidates, we search for at least N 1 consecutive measurements above N 2 × σ from the median, i.e., greater than   N 2 in the standardized light curve.Lower values of N 1,2 recover more low-amplitude flares but increase the number of false positives due to artifacts, and vice versa.We empirically select a threshold of N 1 = 2 and N 2 = 3, compromising between the efficiency of recovering flare events and minimizing false positives.
Although the above criteria are sufficient for identifying flare candidates, the start of their rising phases and the end of their decaying phases may be cut off if they lie below the 3σ threshold.We therefore extend the profiles of flare candidates, similar to the extension process in Section 3.3, to provide a more accurate estimation of their parameters.Considering the sharp rise and gradual decay of flare events, we extend the left and right sides of each candidate event until one and two consecutive epochs, respectively, fall below σ MAD from the median, i.e., less than 1 in the standardized light curve.

Rejection of CVs and SSO Encounters
Cataclysmic variables (CVs) have also been proposed in the TASC WG8 target list as they consist of a WD primary and a mass transferring secondary.However, it is difficult to identify flare events in the light curves of CVs since they exhibit rich variable properties across different time scales from seconds to millennia (Bruch 2022).The majority of outbursts are not triggered by flare events (e.g.superhumps), leading to severe pollution in our identification.We thus have to disregard those light curves by querying the type of each compact star from the SIMBAD astronomical database4 (Wenger et al. 2000).
When a small, foreground solar system object (SSO) (e.g.asteroid, comet) moves across the aperture mask of a target, it tends to cause a symmetrical profile of increasing flux on the light curve as shown in Figure 4.This phenomenon has been extensively discussed by Pál et al. (2018).The SSO encounter events can be easily misidentified as flare candidates and hence need to be excluded.We use the SkyBoT5 service (Berthier et al. 2006) to identify the false signals caused by SSO encounters, which has been implemented in the Kepler/K2 images (Berthier et al. 2016).For each flare candidate, we perform a cone search in SkyBoT with a radius of 8 TESS pixels, corresponding to 2.8 ′ , at the peak time, which returns a list of known SSOs located in the vicinity of the target.However, a very faint SSO will have little impact on the light curve even if it is identified with flare events (see e.g.Günther et al. 2020).We thus also query the visual magnitude of the SSOs from the JPL Horizons system6 (Giorgini et al. 2001) and only remove those with V mag < 19, which will retain the real flares occurring on the targets passing by faint SSOs.

FLARE CANDIDATES VALIDATION
While the criteria described above are effective in identifying preliminary flare candidates, they may also capture other types of transient events or artifacts, resulting in a considerable number of false positives.Despite our efforts to minimize them through the aforementioned steps, some non-flare events may still exist due to various reasons such as poor data quality, uncategorized CVs in SIMBAD, and targets affected by unknown SSOs.As a result, a more rigorous validation process is required to eliminate non-flare events from the preliminary flare candidates.
Owing to the substantial number of preliminary candidates, we adopt Random Forest (RF; Breiman 2001), a supervised machine learning algorithm, to automatically evaluate the confidence level of each candidate.We generate batches of flare and non-flare events through a series of simulations, which serve as the training inputs for the RF classifier.The classifier then evaluates each candidate and returns a confidence probability for it being an intrinsic flare event.Candidates with a probability that exceeds a certain confidential threshold are accepted as validated candidates.

Random Forest
Random Forest (RF) is an ensemble machine learning method that can be used for both regression and classification problems.It is constructed by combining multiple decision trees, with each tree being trained on a different subsample of the training set generated by bootstrap aggregation, also known as bagging (Breiman 1996).For classification tasks, the RF classifier outputs the class probabilities of an input sample by averaging the predicted class probabilities from each tree in the forest.Unlike a single decision tree, which can be unstable and sensitive to noise, RF combines the results of numerous weakly correlated trees, providing more accurate and robust predictions.
Due to its high accuracy, low variance, and ease of application, RF has been widely used for various classification tasks (e.g.Brink et al. 2013;Wyrzykowski et al. 2015Wyrzykowski et al. , 2016;;Godines et al. 2019).Previous literature has also reported that RF outperforms other machine learning methods when classifying variable stars (see, e.g., Richards et al. 2011;Brink et al. 2013;Pashchenko et al. 2018).Given these successful applications and outstanding performance in comparisons, we choose to apply RF to validate the preliminary flare candidates.For this research, we use scikit-learn (Pedregosa et al. 2011), an open source Python package for machine learning, to build the RF classifier.

Training Set
As a supervised learning method, RF requires a training set that is composed of objects with known labels.We choose to construct the training set from simulations rather than relying on data from previous surveys.Using simulations enables the creation of a comprehensive and accurate training set with minimal label contamination because it allows for precise control of input parameters and the generation of a broader variety of flare events, which might be underrepresented in previous survey data due to incomplete samples of low-energy flares (Feinstein et al. 2020).The M22 flare template used the convolution of a Gaussian and a double exponential to model the morphology of the flares, which can be parameterized with three variables: amplitude, the full time width at half the maximum of the flux (FWHM, also known as t 1/2 in Kowalski et al. 2013) and center time (similar to t peak in Davenport et al. 2014, the moment the flare peaks).The flare template was normalized to a relative flux scale, ranging from 0 (before and after the flare occurs) to A, the amplitude (at the flare peak).By changing the parameters from specific distributions, we can generate various simulated flare events using this flare template.

Parameters Fitting for Simulated Flares
To create simulated flares that accurately represent real ones as close as possible, we use observed data to fit the distributions of flare parameters.By drawing parameters from these fitted distributions, we can generate simulated flare events that are statistically representative of real ones.Since the center time does not affect profiles of flare, we completely concentrate on determining the distributions of its amplitude and FWHM.
We specifically choose to use the parameters of the 3792 M-dwarf flares listed in H22, as they provide the photometric signal-to-noise ratio (S/N), σ peak , of the flare peaks, which corresponds to the amplitude A in the M22 flare template.By using σ peak to represent the amplitude, we ensure that the simulating flux is directly comparable to the standardized flux of the preliminary flare candidates, which enables us to use the standardized flux as input for the RF classifier without any additional transformation.Additionally, since H22 provides only the rise and decay times of the flares without their FWHM, we apply an approximate method to obtain them.Based on the M22 flare template, the rise time is roughly equal to 0.6×FWHM.We therefore divide the rise time of each flare by 0.6 to estimate the corresponding FWHM for each flare.
The flare amplitude and FWHM generally exhibit a strong correlation, with higher amplitude flares tending to have longer FWHMs.In H22, the histograms of these parameters follow skewed distributions with long tails toward larger values, characteristic of log-normal distributions (Figure 5).Since taking the natural logarithm of a log-normally distributed random variable yields a normal distribution, we transform the amplitude and FWHM values by taking their natural logarithms.This allows us to simplify the fitting by using a bivariate normal distribution for the transformed data.A bivariate normal distribution is a two-dimensional generalization of the normal distribution, parameterized by the means, standard deviations, and correlation coefficient of the two variables.Specifically, we fit for the means of the natural logarithms of amplitude (µ A ) and FWHM (µ FWHM ), the standard deviations of the natural logarithms of amplitude (σ A ) and FWHM (σ FWHM ), and the correlation coefficient (ρ) between the natural logarithms of amplitude and FWHM.We use the Python package emcee (Foreman-Mackey et al. 2013) to perform a Markov Chain Monte Carlo (MCMC) analysis to fit the model.The walkers are initialized by slightly perturbing the median and standard deviation of the transformed data shown in Figure 5 for the means and standard deviations of the natural logarithm of amplitude and FWHM.The initial correlation coefficient ρ for each walker is randomly drawn from a uniform distribution over [0, 1].We run emcee using 128 walkers for 10 000 steps, discarding the first 1000 steps as burn-in, which is sufficient for the chains to converge.The acceptance fraction is 0.551 and the mean auto-correlation time is 51.5 steps.The resulting estimate for ρ is 0.093 ± 0.016, while the estimates for the other parameters are shown in Figure 5.

Generating Simulated Events
Using the best-fit parameters obtained from the MCMC analysis, we define a bivariate normal distribution to generate paired samples of the natural logarithm of amplitude and FWHM values, which are then input into the M22 flare model to produce a series of simulated flare events.As the flux of these events is standardized, it is essential to introduce noise modeled by a standard normal distribution to create a more realistic representation of the observed data, ensuring that the simulation closely mimics the properties observed in intrinsic flare data.
In addition to simulating flare events, we also generate simulations of non-flare events to establish a contrasting group in the training set.To distinguish these events from flares, which exhibit an asymmetric profile characterized by a rapid rise and gradual decay, we employ a symmetric Gaussian profile parameterized by sigma and amplitude.The sigma values are derived from a log-normal distribution with parameters µ = 3 and σ = 1, while the amplitude values are sampled from a half-normal distribution with parameter σ = 1.Through this parameterization approach, the simulated flare events have a broader range of amplitude values that extend to considerably higher levels, and generally exhibit increased duration compared to non-flare events.These discrepancies in amplitude and duration concur with the intrinsic divergences between flares, which typically demonstrate substantial brightness enhancements over prolonged durations, and non-flares such as noise and artifacts, which display more subtle amplitude variations over shorter timespans.As with the simulated flare events, the simulated non-flare events are incorporated into noise modeled by a standard normal distribution.Both types of simulated events have a 2-minute cadence, which is consistent with the TESS photometry we used here.
Before incorporating the simulated events into the training set, we perform additional steps to ensure their close resemblance to the preliminary flare candidates.In accordance with the criteria used for identifying preliminary flares in Section 3.4, we require the generated events to have a minimum of two consecutive epochs with flux values greater than 3. Subsequently, we extend these epochs using the previously described algorithm, which involves expanding the left and right sides of each candidate event until one and two consecutive epochs, respectively, have flux values less than 1, and then truncating the left epochs.Moreover, for non-flare events with five or more epochs, we require that the peak epoch must be located on the right side of the event, thereby creating a clearer distinction between flare and non-flare events in the training set.
We generate a total of 5,000 simulated events for flares and non-flares, separately.These events are subsequently partitioned into training, validation, and test sets following a 3:1:1 ratio.A selection of events from the training set is illustrated in Figure 6.

Feature Selection and Extraction
Feature extraction is essential in machine learning applications, as it converts time-series data with varying durations into structured vectors suitable for input into the classifier.Moreover, the careful selection of an appropriate and informative feature set is also crucial to ensure the critical information embedded within the data.The chosen features should effectively encapsulate the intrinsic characteristics of the events, allowing the classifier to distinguish more accurately between distinct categories and yield more reliable results.
In this study, we employ the Python package tsfresh (Christ et al. 2018) to facilitate feature selection and perform feature extraction.This package is designed for machine learning applications and is capable of automatically identifying and extracting relevant features from time-series data using hypothesis testing.We initially use tsfresh to select a comprehensive list of features, from which we subsequently handpick a set of valuable features that efficiently distinguish between flare and non-flare events.These selected features are presented in Table 2.
For the training of the RF classifier, we transform each simulated event into a single 9 × 1 vector, corresponding to the selected features.The feature extraction process not only vectorizes the events but also emphasizes the distinguishing characteristics that differentiate flare and non-flare events, enabling the classifier to focus on the most significant aspects of the data.

Hyperparameter Tuning and Model Training
Following feature extraction, we proceed with hyperparameter tuning, a crucial step in optimizing the classifier's performance.This process involves determining the optimal combination of hyperparameters, which are selected based on their capacity to maximize the score of the RF model on the validation set.
We employ the Python package optuna (Akiba et al. 2019) to tune the following hyperparameters of the RF classifier within the specified ranges, considering a balance between model complexity and computational efficiency: • n estimators: The number of trees in the forest.
Range: (100, 1000, step=100) • max depth: The maximum depth of the trees.Range: (3, 17, step=2) • max features: The number of features to consider when seeking the best split.Range: (2, 5, step=1) The GridSampler in optuna is used to perform a grid search, training the RF classifier with all possible combinations of hyperparameters and evaluating their accuracy scores on the validation set.The following hyperparameters produce the highest accuracy score (0.993) on the validation set: • n estimators: 100 • max depth: 7 • max features: 2 We select the corresponding model as the final model.Its performance is assessed on the test set, yielding an F1 score of 0.994, which is calculated as the harmonic mean of precision and recall.Moreover, we can extract the feature importance after the training process, which assigns an importance score to each feature based on its contribution to the decisionmaking process.The significance of our selected features is listed in Table 2, which represents the relative contribution of each feature to the classification process, offering insights into their relevance in the context of flare identification.The first location of maximum feature, capturing the relative location of peak flux, emerges as the most important, which is reasonable given that flares exhibit a rapid rise and gradual decay.

Validating Preliminary Flare Candidates
After training the RF classifier, we apply it to the preliminary flare candidates.As outlined in Section 4.3, we extract the corresponding set of features for these candidates.In cases where the candidates contain missing values, we employ a linear fitting method to interpolate them prior to feature extraction.The resulting feature vectors are then fed as input into the RF classifier.
The RF classifier computes a probability score for each candidate, indicative of its likelihood of being an intrinsic flare.We filter the preliminary flare candidates by applying a probability threshold of 0.5.Candidates with a probability score above this threshold are classified as validated flare events, while others are considered false positives.This filtering effectively refines the flare candidate list, retaining only the events with a high probability of being intrinsic flares for further investigation.
To further improve the reliability and accuracy of the flare validation process, we perform an additional manual verification of the filtered flare candidates through visual inspection.During this process, we identify several common sources of false positives, including glitches in the light curves, unclassified CVs in the SIMBAD database, SSO encounters not captured in the SkyBoT query (potentially due to unknown SSOs or bright SSOs passing by without falling within the queried cone area), and low amplitude candidates exhibiting flare-like profiles but suffering from poor light curve quality.

RESULTS
From a total of 38 102 available 2-minute cadence light curves of 11 618 compact stars, we identified 7584 events as preliminary flare candidates.After excluding 578 events that are attributed by SSOs, we validated the remaining 7006 events with the help of our RF classifier, which flagged 3558 of these events as false positives.Finally, we confirmed 1016 flares to be real events with further visual inspection of the remaining 3448 candidates.These confirmed flare events originated from 193 compact targets, which include 182 events from 58 sdB/sdO stars and 834 events from 135 WDs.Table 3 provides the observed properties of all confirmed flare events, such as peak time, amplitude, and energy in the TESS bandpass (see Section 5.1.2).Table 4 presents the key attributes of the flaring compact stars, including their classification, flare occurrence frequency, and logarithmic fractional flare luminosity (see Section 5.1.3).

Flare Properties and Fractional Flare Luminosity
For each identified flare event, we measure properties including duration, amplitude, equivalent duration, and energy in the TESS bandpass.We also compute the fractional flare luminosity for each of the flaring compact stars to quantitatively measure their strength of flare activity.Figure 7(a) shows histograms comparing the distributions of flare properties and fractional flare luminosity between WDs and sdB/sdO stars, while comparisons exclusively between compact stars with a pollution level of 0 (see Section 5.2.1) are shown in Figure 7(b).

Flare Amplitude and Equivalent Duration
The flare amplitude, defined as the maximum increase in flux during the flare event relative to the flux of the star in its quiescent state, is calculated as: where A represents the flare amplitude, F peak is the maximum flux observed during the flare event and F quiescent is the quiescent flux.
The equivalent duration (ED; Gershberg 1972) provides an estimate of the energy output of a flare relative to the quiescent emission of the star.Expressed in units of time (e.g., seconds), the ED represents the duration for which the star would need to emit at its quiescent brightness level to equal the excess energy released by the flare.It is calculated as follows, where the integral is computed by a trapezoidal sum of the light curve between the start and stop times of the flare.The uncertainty of ED is calculated following Davenport (2016).

Flare Energy in the TESS Bandpass
Typically, the total energy emitted by a stellar flare (i.e., the bolometric flare energy) can be estimated if the effective temperature and radius of the star are known (see, e.g., Shibayama et al. 2013;Günther et al. 2020).
However, the result of the cross-matching of the flaring compact stars with Gaia DR3 shows that only a small subset (∼10%) of our samples have these stellar parameters available, which restricts our ability to calculate the bolometric flare energy across the entire sample.
Given these limitations, we refocused towards estimating the flare energy in the TESS bandpass (E TESS ), a measure that provides insight into the observable energy released during the flare event.It is noteworthy that the bolometric flare energy can be significantly greater, by a factor of approximately five or more, than the flare energy estimated in the TESS bandpass (see H22).This discrepancy arises because flares, often approximated by a 9000 K blackbody (Jackman et al. 2023), primarily emit energy in the high-energy range at short wavelengths, which fall outside the TESS bandpass.
We first calculate the luminosity of the star in the TESS bandpass using the following formula: where d is the distance to the star derived from its parallax listed in Gaia DR3.The flux of the star in the TESS bandpass F TESS is calculated via the following relation: where T is the TESS magnitude derived from the TESS Input Catalog version 8.2 (TIC v8.2; Stassun et al. 2019;Paegert et al. 2021) and F 0 = 4.03 × 10 −6 erg s −1 cm −2 , which is the flux corresponding to T = 0 (Sullivan et al. 2015).The uncertainty in the calculated luminosity L TESS accounts for uncertainties in both the TESS magnitude and the parallax.Lastly, we compute the flare energy in the TESS bandpass by: where ED is the equivalent duration of the flare and L TESS is the calculated luminosity of the star in the TESS bandpass.
For seven of the flaring compact stars, which correspond to 121 flare events, the parallax parameters were not available in Gaia DR3.As a result, we were only able to calculate the energies for the remaining 923 flares.The median value of uncertainties associated with E TESS is 22.1%.

Fractional Flare Luminosity
For each flaring compact star, we calculate the fractional flare luminosity, represented as L flare /L TESS , which is the total luminosity emitted by flares relative to that emitted by the star through the TESS bandpass.Serving as an intuitive metric, it quantifies the intensity of stellar flare activity.The metric was originally We calculate this ratio for a star using the following equation, where ED i represents the ED of each flare for the star, and t obs is the total observation time for the star by TESS.A higher value of this ratio signifies enhanced flare activity, indicative of an increased level of stellar magnetic activity.We cross-matched the flaring compact stars with Gaia Data Release 3 (DR3; Gaia Collaboration et al. 2016Collaboration et al. , 2023) ) to acquire their BP-RP color index (G BP − G RP ) and absolute Gaia G magnitude (G abs ), and plot them on the Gaia Hertzsprung-Russell (H-R) diagram (Figure 8).We distinguish between sdB/sdO stars and WDs in the diagram and illustrate the fractional flare luminosity across the flaring compact stars.

Contamination Check
Prior to analyzing the flare activities, it is vital to refine our sample to more accurately characterize the flare activity inherent to sdB/sdO stars and WDs.This refinement ensures the exclusion of targets where flare activities might potentially originate from contaminating sources other than the compact stars themselves, such as nearby objects or companion stars of the compact target.

Evaluating Pollution Level
Due to the relatively low angular resolution of TESS (21 arcsec pixel −1 ), there is a considerable risk that the photometry may be contaminated by nearby objects.To address this concern, we developed an open-source script, tpfi, which generates identification charts for TESS target pixel files, as illustrated in Figure 9.Our script, based on tpfplotter7 (Aller et al. 2020), introduces a notable feature: the ability to visualize the identical sky coverage of the target pixel file as provided by the Digital Sky Survey (DSS)8 .This enhancement allows for a more convenient assessment of the contamination level of a target.In addition to its applicability to TESS, we extended the use of tpfi to Kepler and K2 missions.We anticipate that this enhancement will make the script valuable to broader communities, for instance, the exoplanet and variable star research communities.Our script is publicly available on Github9 .
Following these considerations, we employed tpfi to create identification charts for all compact stars in our sample.We thus can conduct a visual examination of each of the 193 compact stars in our sample to assess the potential for contamination from nearby sources.In order to systematically quantify the contamination, we defined the pollution level (PL) for each star, as shown in Table 4.A PL=0 indicates that the target is the only object within the aperture, without any significant polluting photometry from nearby stars.A PL=1 denotes that the target is the brightest object within the aperture, but there are other dim stars present.Lastly, a PL=2 refers to the cases where the target is not the brightest object within the aperture or near a much brighter star.To illustrate, Figure 9(a) represents a scenario with a PL=0, indicating an unpolluted target star.Figure 9(b) displays minor contamination, warranting a PL=1.Cases of severe contamination, as demonstrated by Figures 9(c) and (d), are assigned with a PL=2.
In our categorization, 44 stars have PL=0, 62 stars fall into the PL=1 category, and 87 stars are classified with PL=2, which corresponds to 279, 256, and 481 flares, respectively.These pollution levels offer a crude estimate of potential contamination.They are particularly useful in downstream analyses, where the potential impact of such contamination on our flare detection results must be considered.We also note that a high PL does not necessarily disqualify a star as a flare candidate.However, stars with a high PL warrant additional caution and follow-up investigation to confirm the source of the flare events.

Detecting companion stars
To address potential contamination from companion stars, we firstly exclude flaring stars that are labelled as binary in TASC WG8 target list.We then locate the rest flaring stars on the Gaia H-R diagram.If a target falls on the main sequence, it indicates the compact star has a brighter MS companion dominating the observed brightness, which is especially relevant for the intrinsically faint WDs.Therefore, we exclude the flaring compact stars positioned within the main sequence of the Gaia H-R diagram (G BP − G RP > 1 and G abs < 12, see Figure 8).We also exclude the stars that are not provided with G BP − G RP or G abs in Gaia DR3, as their positions on the H-R diagram are unavailable.After the filtering, 13 flaring compact stars remained, comprising seven WDs and six sdB/sdO stars.
To further examine the likelihood of companion stars, we constructed the spectral energy distributions (SEDs) for these stars using the VO Sed Analyzer (VOSA)10 ( Bayo et al. 2008).Our examination revealed that all seven WDs display a significant red/near-infrared flux excess, suggesting the presence of cool companions.In contrast, the SEDs of the six sdB/sdO stars displayed less noticeable red/near-infrared excess.This may be attributed to the inherent high brightness of sdB/sdO stars, which can dominate the SED when paired with a cool MS star.We then employed speedyfit11 to fit the SEDs of the sdB/sdO stars.Initially, we used spectral models exclusively from the Tübingen NLTE Model-Atmosphere Package (TMAP; Werner & Dreizler 1999;Werner et al. 2003;Rauch & Deetjen 2003).However, the fitting residuals indicated a clear IR excess beyond what the model predicted for all six sdB/sdO stars.This discrepancy was resolved when we incorporated a second component in the fitting, modeled using ATLAS9 spectra (Castelli & Kurucz 2003).The results of the MCMC fitting suggest that the effective temperatures of the cool companions are within the range of 3500 K to 4500 K, consistent with K/M-dwarfs.A more detailed description of the fitting process of speedyfit is given in Vos et al. (2017Vos et al. ( , 2018)).

Selecting the Refined Sample
Our previous analysis found that all flaring compact stars show evidence of contamination from nearby objects or the presence of an unresolved companion.As a result, we have been unable to conclusively attribute any specific stellar flares solely to an individual compact star.However, the origin of the stellar flares in binary systems containing compact stars is still unclear and demands careful examination.
For these flaring compact stars with cool companions, the origin of the flares warrants careful consideration, since flares are common in cool MS stars, particularly M-Xing et al. dwarfs.Therefore, for WDs, there is a high probability that the flares originate from their red companions, especially given the comparable luminosity levels of WDs and cool MS stars.This makes it reasonable for flares from the cool companions to be detectable in the combined light curve of the WD binary systems.However, the situation is less clear-cut for sdB/sdO systems.The comparative high brightness of sdB/sdO stars over cool MS stars suggests the flare contribution from a cooler companion would be less apparent.If the stellar flares indeed originate from the cool companion, they would need to be extraordinarily energetic to be noticeable in the combined light curve.Such high-energy flares from a cool MS companion are relatively rare, implying that the flares observed in sdB/sdO systems may have a different origin.Another possible explanation is that the companions to compact objects are more magnetically active compared to their isolated counterparts, causing more frequent high-energy flares.
To conduct a more detailed investigation into the origin of flares in sdB/sdO systems, a refined sample of sdB/sdO stars is necessary for further in-depth analysis of flare activity.To ensure the purity of the refined sample, we cross-match all 16 flaring sdB/sdO stars with PL=0 with the known hot subdwarf catalog in Culpan et al. (2022).There are three stars not present in the known hot subdwarf catalog, which are misclassified hot subdwarfs or hot subdwarf candidates, but labelled as hot subdwarf in TASC WG8 target list.This left us with a refined sample of 13 sdB/sdO stars, corresponding to 23 flare events.The detailed parameters of the 13 selected sdB/sdO stars are presented in Table 5, while the 23 flare events are shown in Fig. 11.Xing et al.

Flare Frequency Distribution
We here investigate the Flare Frequency Distributions (FFDs) of different subsets of the flaring compact stars.The FFD, typically described by a power law (Lacy et al. 1976), represents the frequency of flare occurrences as a function of their energy.It can be expressed as, where N represents the number of flares occurring within a specific observation duration, E denotes the flare energy, k is a proportionality constant, and α is the power-law index (Jackman et al. 2021).This powerlaw index, α, plays a critical role in investigating flare production mechanisms, and its derivation is essential to our study of flare activities in compact stars.Due to differences in observation durations for various flaring compact stars, a direct computation of their FFD is not feasible.We hence employ the cumulative FFD to estimate α instead.The cumulative FFD is derived by integrating Equation ( 9), leading to where β = log( k 1−α ), and ν denotes the frequency of flares, i.e. the number of flares per unit time with energy exceeding a specific threshold.
We compute three distinct cumulative FFDs, each representing flare events from all observed sdB/sdO stars, WDs, and the refined sample of sdB/sdO stars (see Section 5.2.3).For every flare event with energy E, we calculate the corresponding ν as follows: where N i (> E) is the number of flares with energies greater than E and t i is the observation time for each flaring compact star.Subsequently, we use the MCMC method to fit a linear regression model to each cumulative FFD.This approach effectively mitigates binning effects, providing an advantage over the conventional method of histogram generation and straight-line fitting (Maschberger & Kroupa 2009).We only fit the portion of the cumulative FFDs where we consider the sample complete.events from 193 compact stars (Table 3 & 4).We pioneered a new method for flare detection, specifically designed to address short-term periodic variations in light curves.This method enhances flare detection capabilities when applied to the light curves, despite the complexities introduced by common phenomena in compact stars such as pulsation and binary effect.We also considered potential interferences from CVs and SSOs and eliminated them.In addition, we implemented a validation step using machine learning to effectively filter out false positives.These rigorous detection and validation processes assured the reliability of the detected flares, enabling us to establish the first flare catalog of compact stars, which contains 58 sdB/sdO stars and 135 WDs with 182 and 834 flare events, respectively.
We then thoroughly examined potential contamination caused by nearby objects and companion stars of the compact stars.This included developing an opensource script, tpfi, to generate identification charts for TESS target pixel files, thereby enabling a detailed visual inspection for potential contamination sources near the target.We also conducted analysis using Gaia DR3 data to identify binary systems with a bright MS com-panion, as well as SED fitting to reveal infrared excess from cool companion stars.Through these extensive analyses, we found evidence that all flaring compact stars showed signs of contamination from nearby objects or unresolved companion stars.As a result, we cannot conclusively attribute any specific stellar flares solely to an individual compact star at this stage.
For WDs, there is a high probability the observed flares originate from the cool MS companion.With comparable luminosities and vigorous magnetic activity generating frequent energetic flares, cool companions (e.g.K/M dwarfs) can overwhelm the emission of the WD.Considering the large number of 7505 WDs we searched for stellar flares, we still did not confidently detect flares intrinsic to any individual WD.This may partly be attributed to the low angular resolution of TESS causing severe pollution in the light curves (only ∼ 25% flaring compact stars have a PL=0).However, our detection provides some indication that the likelihood of observable flares occurring on WDs may be low.Alternatively, if stellar flares do occur on WDs, they may operate via a different mechanism with shorter timescales on the order of the dynamical timescale of WDs (∼seconds), which cannot be captured by TESS cadence.
Nevertheless, the origin of stellar flares in sdB/sdO+MS binaries still warrants investigation.With significantly greater luminosity, a flare from the cool MS companion would need to be extraordinarily energetic to be detectable in the combined light curve from the system, which is very rare.This suggests the observed flares may be not from the cool MS companions.We aim to explore this scenario further through analyzing our refined sample of sdB/sdO stars.By focusing on high-purity sdB/sdO stars with negligible contamination, we can conduct detailed analyses to understand the origins of these flare events occurring in sdB/sdO binaries.We finally selected a refined sample consisting of 13 sdB/sdO stars, corresponding to 23 flare events (Table 5).
Distributions comparing various flare properties and fractional flare luminosity between WDs and sdB/sdO stars are shown in Figure 7. Figure 7(a) displays the overall sample, while Figure 7(b) shows the subset with lowest pollution level from nearby objects (PL=0).We note that all compact stars in our sample show evidence of contamination from nearby objects or companion stars, substantially impacting the results.Despite pollution concerns, no obvious distinctions emerge in duration, amplitude or ED distributions when contrasting WD and sdB/sdO flares, although these distributions are more concentrated in Figure 7(b).Higher flare energies in sdB/sdO stars (Figure 7(a)) may link to contam-ination by nearby active stars, as energy calculations use stellar luminosity (see Eqn. 7), which is greater in sdB/sdO stars.The diminishment of this discrepancy for the PL=0 subset in Figure 7(b) supports this contamination explanation.Furthermore, higher fractional flare luminosities for WDs could connect to their lower luminosity compared to sdB/sdO stars, allowing the flares from MS stars to more readily override and become detectable, and thus show higher level of flare activity.Again, this potential contamination effect diminishes when analyzing the PL=0 subset.Therefore, while definitive compact star flare attribution remains complex presently, comparing WD and sdB/sdO flare distributions provides initial insights on their distinct origins.
When we proceeded to analyze the FFDs, a notable observation emerged related to the power-law index α (see Figure 10).We found that the FFDs from both sdB/sdO stars and WDs have an index α ∼ 2 when fitted with the entire flare sample, which is polluted by flares originated from other late-type MS stars (F-M type stars).These MS stars, with FFDs having an index α ∼ 2 (e.g., Althukair & Tsiklauri 2023), significantly skew the index α when fitted with the entire sample.However, when we move to the refined sample of sdB/sdO stars, we obtain an FFD that is less steep (α = 1.70), which indicates a higher proportion of high energy flare events.
Before delving into potential explanations for this phenomenon, it is pertinent to note that several studies have investigated stellar flares across spectral types A to M. They found that the FFD from A-type stars has an index 1 < α < 1.5 when using Kepler photometry, a finding that markedly deviates from the α ∼ 2 observed in cooler F-M type stars ( Švanda & Karlický 2016;Yang & Liu 2019;Bai & Esamdin 2020;Althukair & Tsiklauri 2023).Yang et al. (2023) also reported α = 1.76 ± 0.19 for A-type stars based on TESS data.Although this value is higher than that obtained via Kepler photometry, which is possibly due to increased contamination from lower angular resolution of TESS, it remains lower than indices for cooler MS stars.Such deviations in α have been interpreted as indicative of differing flare mechanisms between early-type (B/A-type) and latetype MS stars.
This observation prompts us to compare sdB/sdO stars with B/A-type MS stars, especially given the observed deviation in α between the FFDs derived from our entire sample and the refined sample of sdB/sdO stars.Notably, both sdB/sdO stars and B/A-type MS stars possess a radiative envelope, distinguishing them from cooler MS stars with convective envelopes.Despite the general expectation that these stars lack strong magnetic fields and are unlikely to produce flares through known dynamo mechanisms, as a result of the absence of a deep convective envelope (Charbonneau 2010), previous literature has reported that they might do show magnetic activities.For instance, Balona (2021) reported flare events on B/A-type MS stars and argued that observed rotational modulation in flaring B/A-stars suggests strong surface magnetic fields.
We thus hypothesize that some of the stellar flares from the refined sdB/sdO sample may originate solely on sdB/sdO itself, or through magnetic reconnection involving the sdB/sdO and its close companion.This hypothesis stems from the low incidence of superflares in late-type MS stars, the similar structure between hot MS stars and sdB/sdO stars, and most importantly, the decrease of the power-law index α from the FFD of the refined sdB/sdO sample compared to other samples.Analogous conclusions were previously proposed for B/A-type MS stars (Balona 2021;Maryeva et al. 2023).We also propose stellar flares on sdB/sdO stars may operate through similar mechanisms to those detected in hot MS stars, if the flares in the refined sdB/sdO sample do arise exclusively from the sdB/sdO star.However, the nature of magnetic fields responsible for spots and flares in hot MS stars remains an open question.Švanda & Karlický (2016) assumed these fields could be amplified by dynamo processes in the convective cores of A-type stars, subsequently becoming unstable and rising as magnetic ropes through the radiative envelope.Balona (2019) suggested differential rotation might suffice to generate local magnetic fields in B/Atype stars, presenting an exciting direction for future research to deepen our understanding of magnetic activities in stars with radiative envelopes.It is crucial to emphasize, however, that our hypothesis is preliminary and warrants further in-depth investigation.
We recall that, to our knowledge, there has been no dedicated survey searching for stellar flares in compact stars previously, although some searches have been conducted focusing on the outbursts in WDs (see e.g., Bell et al. 2016).The lack of extensive surveys is largely due to challenges like complex light curve detrending (Pietras et al. 2022).Our catalog could be the first step to such research, definitely, triggering new interests in stellar activity in highly evolved compact stars.The 13 sdB/sdO stars in the refined sample could be good candidates for future inspection.If confirmed, these candidates would be the first compact stars to exhibit a flare event.Moreover, while we cannot yet conclusively attribute any flares solely to an individual compact star, characterizing flare events in compact binary systems merits deeper investigation.For instance, Morgan et al. (2016) demonstrated enhanced magnetic activity in M dwarfs with close WD companions compared to their isolated counterparts by analyzing their flare rates.Our results can help to examine such magnetic interactions across various compact binary systems.
We finally propose some prospects for our discoveries of flare events in sdB/sdO stars and WDs.Our methods are readily adaptable for similar analyses of Kepler and K2 photometry, which boast a higher angular resolution of roughly 4 arcsec pixel −1 , reducing contamination from nearby objects significantly.We anticipate the ongoing photometry from the second extension mission of TESS, which will further enable continuous monitoring of flare events in the entire sample.In addition, our method and pipeline can be helpful for flare hunting in other types of stars, or detecting other types of transient events in compact stars, with high confidence and feasibility.Furthermore, our tool, tpfi, is wellintegrated with Kepler/K2 photometry, which enhances contamination reduction.Beyond the research on stellar flares, the identification charts generated by tpfi can also be invaluable for other studies, for instance, exoplanet detection and variable star research, underscoring its broader astronomical applicability.

Figure 1 .
Figure 1.Detrending the light curve of TIC 219244444 (RR Cae, a WD eclipsed by an M-dwarf every ∼0.3 d; Maxted et al. 2007) in TESS Sector 3 using the biweight filter with two different window lengths: (a) long of 0.3 d, (b) short of 0.05 d.Top: raw flux (black points) and trend flux (red solid line).Bottom: detrended flux (black points), 3σMAD from the median (red dashed line) and zoom windows of one small flare (left) and one large flare (right).The large flare is overfitted with the short window, while the small flare cannot be detected with the long window.Only a 6-day segment of the light curve is shown in this figure.

Figure 2 .
Figure 2.An example of generating the model of the short-period variability in the light curve of TIC 219244444.Top: The light curve of TIC 219244444 in TESS Sector 3 (left) and its Lomb-Scargle periodogram (right).The red arrow indicates the period used to fold the light curve (P fold = 0.304 d), which is twice the period at the max power.Bottom: The folded light curve and its model (red solid line) obtained by applying the median filter (left) and the residual of the model (right).The epochs masked by the processing in Section 3.3 are marked with grey points.The model fits the short-period variability well except the ingress and egress of the eclipse.

Figure 3 .
Figure3.The same as Figure1, but using our detrending procedure instead of directly using the biweight filter.The trend flux (red solid line) is the combination of the unfolded model of the short-period variability and the trend flux obtained using the biweight filter afterwards.The grey points in the top panel are removed (see Section 3.3) in the detrended light curve, which causes the missing data (grey dashed line) in the right zoom window of the bottom panel.

Figure 4 .
Figure 4.An SSO encounter event found in the light curve of TIC 275308213 in TESS Sector 8.The light curve shown in the left panel is the photometry result using the aperture shown as orange region in the right panel, where presents four frames of the target pixel file during the process that an SSO, marked by red star symbol, encounters the target star.The corresponding fluxes of these moments are also marked in the light curve.all 11 months of 1-minute cadence data for GJ 1243, an active M4 star, available from Kepler Data Release 23.This flare template has been widely used for constructing light curves of flare stars (e.g.Günther et al. 2020).Recently, Tovar Mendoza et al. (2022, hereafter M22) reanalyzed the same data for GJ 1243 using the Kepler Data Release 25, where the light curve processing was improved.They generated an updated analytic and continuous flare template, addressing the limitations of the flare template in Davenport et al. (2014).We thus adopt the flare template in M22 to generate the simulated flare events here.The M22 flare template used the convolution of a Gaussian and a double exponential to model the morphology of the flares, which can be parameterized with three variables: amplitude, the full time width at half the maximum of the flux (FWHM, also known as t 1/2 inKowalski et al. 2013) and center time (similar to t peak inDavenport et al. 2014, the moment the flare peaks).The flare template was normalized to a relative flux scale, ranging from 0 (before and after the flare occurs) to A, the amplitude (at the flare peak).By changing the parameters from specific distributions, we can generate various simulated flare events using this flare template.

Figure 5 .
Figure 5. Histograms of the flare amplitude A and FWHM values, along with their natural logarithms, for the flare samples from H22. Left panels: The distributions of the natural logarithms of amplitude and FWHM (black step lines) with the fitted normal distributions (red solid lines).Right panels: The distributions of the original amplitude and FWHM values (black step lines) with the fitted log-normal distributions (red solid lines).The median values (vertical dashed lines) and standard deviations of the natural logarithms of amplitude and FWHM, provided in the upper right of the left panels, are used to initialize the walkers for the MCMC analysis.The resulting estimates from the MCMC analysis are shown in the upper right of the right panels.

Figure 6 .
Figure 6.Examples of simulated events in the training set: (a) Flares, (b) Non-flares.The dashed lines and dotted lines indicate standardized flux values of 3 and 1, respectively.

Figure 7 .
Figure 7. Histograms comparing flare properties and fractional flare luminosity between WDs (blue lines) and sdB/sdO stars (red lines).The top panels (a) show distributions for all flare events and flaring compact stars, while the bottom panels (b) display only flare events and flaring compact stars with a PL of 0. The left four panels in both rows show flare property distributions: duration, logarithmic amplitude, logarithmic ED, and logarithmic energy in the TESS bandpass, either for (a) 834 WD and 182 sdB/sdO flares or (b) 193 WD and 86 sdB/sdO flares.The right panels illustrate distribution of the logarithmic fractional flare luminosity value, either for (a) 135 flaring WDs and 58 flaring sdB/sdO stars or (b) 28 flaring WDs and 16 flaring sdB/sdO stars.

Figure 8 .
Figure 8. Gaia H-R diagram of the flaring compact stars included in our study.The flaring WDs and sdB/sdO stars are signified by diamond and circle markers, respectively, with other compact stars from TASC WG8 target list denoted with black dots.The colors of the markers correspond to the log(L flare /LTESS) value for each flaring star.The stars surrounded by blue squares are those with a pollution level of 0 (see Section 5.2.1).The red dash-dotted lines indicate the thresholds used for filtering out stars with companions (see Section 5.2.2).

Figure 9 .
Figure 9. Examples of the identification charts provided by tpfi: (a): an ideal situation that the aperture only hosts the target with no other bright stars nearby, (b): the aperture hosting the target star and another dimmer star, (c): the aperture hosting numerous stars with comparable brightness alongside the target, (d): only the target resides within the aperture, but in proximity to a bright star.In each chart, the right panel overlays the Gaia DR3 catalog onto the target pixel file, with the target denoted by a white cross symbol.The size of the circle represents the relative brightness of the stars, as indicated by the Gaia G magnitude.The red region indicates the default aperture mask used by SPOC for photometry extraction.The left panel shows the same sky coverage, using DSS2 red images, with the same orientation.
table lists the 13 sdB/sdO stars in our refined sample.Columns are TIC ID, Gaia DR3 source ID, target name, spectral classification given byCulpan et al. (2022), TESS magnitude from TIC v8.2(Stassun et al. 2019;Paegert et al. 2021), BP-RP color index from Gaia DR3, absolute Gaia G magnitude from Gaia DR3, number of observed flares, flare occurrence frequency per TESS sector, and logarithmic fractional flare luminosity.
Figure 10.Cumulative FFDs and corresponding power law fits.The blue and red open circle markers represent flares from all WDs and sdB/sdO stars, respectively.The red filled circle markers signify flares from the refined sdB/sdO sample.The dashed and dash-dot black lines represent the bestfit power laws to the cumulative FFDs of the full and refined samples, respectively.The fitted α values and uncertainties, along with the number of flares used for fitting N fit , are shown in the bottom left.

Figure 11 .
Figure 11.All 23 flare events (red lines) in the light curves of the 13 sdB/sdO stars in our refined sample.The grey dashed lines are 3σMAD from the median.TIC IDs and Sector numbers are indicated in the top right corner of each panel.The flare event in the light curve of TIC 234281664 in TESS Sector 28 is a complex flare event composed of two classical (single peak) flares.

Table 1 .
Number of evolved compact stars and light curves observed in TESS Cycles 1-5 with data available from N sectors (N from 1 to ⩾ 6).

Table 2 .
(Christ et al. 2018)atures Selected for Classification.Note-These features are computed using the Python package tsfresh(Christ et al. 2018).For details on these features and algorithms, see https://tsfresh.readthedocs.io/en/latest/text/list of features.html.

Table 3 .
Catalog of All 1016 Flares Observed across 38 102 Light Curves of 11 618 Compact Stars at 120 s Cadence during TESS Cycles 1-5.The table presents the details of the 1016 flare events.Columns are TIC ID, Sector, start time, peak time, stop time, signal-to-noise ratio, amplitude, equivalent duration, and energy in the TESS bandpass.The start, peak, and stop time are in Barycentric TESS Julian Date (BTJD).The signal-to-noise ratio (SNR) is the maximum value in the standardized light curve during the flare event.Uncertainties are given in parentheses.(This table is available in its entirety in machine-readable form.)

Table 4 .
Yang & Liu 2019;Ilin et al. 20212019)ars.Note-The table presents the details of the 193 flaring compact stars.Columns are TIC ID, spectral type given by TASC WG8 target list, object type queried from the SIMBAD database, stellar luminosity in the TESS bandpass, number of observed flares, flare occurrence frequency per TESS sector, logarithmic fractional flare luminosity, and pollution level (Poll.) of the target (see Section 5.2.1).The description of the object types is presented in https: //simbad.cds.unistra.fr/guide/otypes.htx.Uncertainties are given in parentheses.(Thistable is available in its entirety in machine-readable form.)denotedasLfl/LKpinLurie et al. (2015)for characterizing the flare activity level of the flaring stars with Kepler photometry.This metric is now widely used in various flare research (see, e.g.,Davenport 2016;Davenport et al. 2019).In certain cases, the L Kp has been substituted with the bolometric luminosity L bol (see, e.g.,Yang & Liu 2019;Ilin et al. 2021).

Table 5 .
Catalog of Flaring Compact Stars in the Refined Sample.