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

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), and 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- to M-type main-sequence (MS) stars (Yang & Liu 2019).

Incapable of resolving the 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 into 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 Transiting Exoplanet Survey Satellite (TESS) missions (Koch et al. 2010; Howell et al. 2014; Ricker et al. 2015), collecting high-quality and long consecutive photometry, have advanced the field of flare studies, exposing various properties and correlations of flares to their hosting stellar types across the Hertzsprung–Russell (H-R) 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 (WDs), remains very limited. Hot subdwarfs are evolved stars with B–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 stars (Heber 2009). WDs are remnants of ∼97% of the 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 convective envelopes, preventing 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, WDs are found with magnetism whose strength is in a large range from a few kilogauss 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 have been claimed to have a detectable magnetic field (Vos et al. 2021; Pelisoli et al. 2022). Therefore, whether these magnetic compact stars can present flare events is an open and 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 that flare hunting requires intensive photometric monitoring. The 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 variations 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 a pioneering survey with the aim of finding flare events in sdB/sdO stars and WDs.

In line with this study, the light curves of compact stars can present complex variability, bringing 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 sporadic outbursts in cool hydrogen-atmosphere pulsating WDs (DAVs; Bell et al. 2015; Hermes et al. 2015; Bell et al. 2016), which can be explained by limit cycles due to sufficiently resonant three-mode couplings (Luan & Goldreich 2018). Scaringi et al. (2022) also reported localized thermonuclear bursts in three accreting WDs. 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 the light curves of compact stars.

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 of 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 flare identification and validation techniques.

This paper is structured as follows. In Section 2, we describe the TESS photometry collected on our compact stars of interest. Section 3 provides details of our methods for preliminary identification of flare candidates, including a new approach designed to address short-term periodic variations, which enhances our ability to identify flares in compact stars. We also describe several indispensable steps to exclude contamination from cataclysmic variables (CVs) and solar system objects (SSOs). Section 4 then explains our process of validating the preliminary candidates using a random forest (RF) 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 (Ricker et al. 2015) was launched on 2018 April 18, and began observing in 2018 July. After completing its primary mission in 2020 July and the first extended mission (EM1) in 2022 September, TESS is currently conducting its second extended mission (EM2), which will last approximately 3 yr.

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 Cycles 1 and 3 pointing to the southern ecliptic hemisphere and Cycle 2 to 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'' pixel⁻¹ over an observational period of ∼27 days. Due to the overlaps among the sectors in each cycle, a portion of the TESS targets were observed multiple times during the survey. Targets in the vicinity of the ecliptic poles were observed continuously for almost 1 yr.

TESS is equipped with four 10.5 cm optical telescopes and four identical cameras with a red bandpass covering the wavelength range of 600–1000 nm. Each camera has a field of view of 24° × 24° and consists of four 2k × 2k back-illuminated CCDs, which continuously read out at 2 s intervals to guide spacecraft. These 2 s data are then stacked, compressed, and stored on the spacecraft until TESS reaches its perigee after every 13.7 days of 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) ¹² 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 a 2 minute cadence that have been collected for these compact stars. We retrieved their light curves from the Mikulski Archive for Space Telescopes (MAST), ¹³ 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 two-thirds 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 flux, where common instrumental systematics have been removed using the co-trending basis vectors. In each light curve, all epochs with quality bit flags equal to 1, 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. ¹⁴

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 CVs, due to their complexity and irregularity. We also inspect whether each candidate is caused by an SSO encounter event.

3.1. 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 constituting small segments of the light curves. Meanwhile, the Savitzky–Golay filter easily overfits the transit signals, which reduces their depth and wraps the detrended light curve around the transit. This will introduce extra outliers and some false alarm detections of flares (see Figures 1 and 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.

Zoom In Zoom Out Reset image size

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. 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 day. We thus set the window length to 0.3 day in the application of the biweight filter, which is considerably long enough for the duration of most flares. Such window length can also effectively filter out other transient events with typical timescales over 0.3 day in compact stars, such as pulsation outbursts in DAVs (Bell et al. 2017).

3.2. Detrending Procedure

Flares typically last from minutes to hours, and rarely longer than 1 day (Ilin & Poppenhaeger 2022). In some cases, the light curves of compact stars exhibit periodic variations, on a timescale of order ≲1 day, 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 a 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 short-term periodic variations.

We first 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–2 days 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 represent the true period of variation in brightness. 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 day. Both the large and small flare profiles are preserved and can be easily detected through this procedure (Figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.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

$\begin{eqnarray}&&{\sigma }_{\mathrm{MAD}}\approx 1.4826\,\mathrm{MAD},\end{eqnarray} \tag{ 1 }$

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).

3.4. 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

$\begin{eqnarray}&&{F}_{{\rm{s}}}=\displaystyle \frac{{F}_{{\rm{d}}}-1}{{\sigma }_{\mathrm{MR}}},\end{eqnarray} \tag{ 2 }$

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 data sets 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₁ consecutive measurements above N₂ × σ from the median, i.e., greater than N₂ 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₁ = 2 and N₂ = 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.

3.5. Rejection of CVs and SSO Encounters

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 timescales from seconds to millennia (Bruch 2022). The majority of outbursts are not triggered by flare events (e.g., superhumps), leading to severe contamination in our identification. We thus have to disregard those light curves by querying the type of each compact star in the SIMBAD astronomical database ¹⁵ (Wenger et al. 2000).

When a small foreground 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 SkyBoT ¹⁶ service (Berthier et al. 2006) to identify the false signals caused by SSO encounters, which have 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 28, 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 system ¹⁷ (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.

Zoom In Zoom Out Reset image size

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 nonflare 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 nonflare events from the preliminary flare candidates.

Owing to the substantial number of preliminary candidates, we adopt RF (Breiman 2001), a supervised machine-learning algorithm, to automatically evaluate the confidence level of each candidate. We generate batches of flare and nonflare 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 of it being an intrinsic flare event. Candidates with a probability that exceeds a certain confidential threshold are accepted as validated candidates.

4.1. RF

4.2. Training Set

As a supervised learning method, RF requires a training set that is composed of objects with known labels. We chose 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 wider variety of flare events, which might be underrepresented in previous survey data due to incomplete samples of low-energy flares (Feinstein et al. 2020).

4.2.1. Flare Model

Davenport et al. (2014) generated an empirical flare template using classical (single peak) flares discovered in 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 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 used 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 time of the full width at half-maximum (FWHM) of the flux (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.

4.2.2. Parameter Fitting of Simulated Flares

To create simulated flares that accurately represent real ones as closely 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 the profiles of flares, we concentrate entirely 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 only provides 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.

Zoom In Zoom Out Reset image size

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 the 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 autocorrelation 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.

4.2.3. 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 nonflare 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 nonflare events. These discrepancies in amplitude and duration concur with the intrinsic divergences between flares, which typically demonstrate substantial enhancement in brightnesses over prolonged durations, and nonflares, such as noise and artifacts, which display more subtle amplitude variations over shorter timespans. As with the simulated flare events, the simulated nonflare 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 performed 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 nonflare 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 nonflare events in the training set.

We generate a total of 5000 simulated events for flares and nonflares 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.

Zoom In Zoom Out Reset image size

4.3. Feature Selection and Extraction

4.4. Hyperparameter Tuning and Model Training

4.5. Validating Preliminary Flare Candidates

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 to 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).

Table 4. Catalog of All 193 Flaring Compact Stars

TIC	Spectral Type	Object Type	LTESS	${N}_{\mathrm{flares}}$	Flare Freq.	$\mathrm{log}({L}_{\mathrm{flare}}/{L}_{\mathrm{TESS}})$	Contam.
			(erg s−1)		(sector−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
6997163	sdB+F/G/K	HotSubdwarf	9.37(1.06)1032	1	0.5	−4.54	0
20656977	DA+M	WhiteDwarf	3.69(0.04)1030	1	1.0	−3.68	0
21860382	DA	WhiteDwarf	9.54(0.07)1029	1	0.333	−3.60	2
23226265	cand-DAV	WhiteDwarf	9.86(0.06)1029	6	3.0	−3.58	1
23385704	WD	WhiteDwarf_Candidate	1.23(0.01)1030	1	0.333	−4.20	2
23992223	DA	WhiteDwarf	6.20(0.05)1030	1	0.5	−4.78	1

Note. The table presents the details of the 193 flaring compact stars. Column (1): TIC ID. Column (2): spectral type given by TASC WG8 target list. Column (3): object type queried from the SIMBAD database. Column (4): stellar luminosity in the TESS bandpass. Column (5): number of observed flares. Column (6): flare occurrence frequency per TESS sector. Column (7): logarithmic fractional flare luminosity. Column (8): level of contamination (Contam.) of the target (see Section 5.2.1). The description of the object types is presented in SIMBAD (https://simbad.cds.unistra.fr/guide/otypes.htx). Uncertainties are given in parentheses.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

5.1. Flare Properties and Fractional Flare Luminosity

For each identified flare event, we measure properties including duration, amplitude, equivalent duration (ED), and energy in the TESS bandpass. We also compute the fractional flare luminosity for each of the flaring compact stars to measure their strength of flare activity quantitatively. 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 level of contamination of 0 (see Section 5.2.1) are shown in Figure 7(b).

Zoom In Zoom Out Reset image size

5.1.1. Flare Amplitude and ED

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

$\begin{eqnarray}&&A=\displaystyle \frac{{F}_{\mathrm{peak}}-{F}_{\mathrm{quiescent}}}{{F}_{\mathrm{quiescent}}},\end{eqnarray} \tag{ 3 }$

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 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 level of quiescent brightness to equal the excess energy released by the flare. It is calculated as follows:

$\begin{eqnarray}&&\mathrm{ED}=\int \displaystyle \frac{{F}_{\mathrm{flare}}\,(t)-{F}_{\mathrm{quiescent}}}{{F}_{\mathrm{quiescent}}}{\rm{d}}t,\end{eqnarray} \tag{ 4 }$

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).

5.1.2. 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 crossmatching of the flaring compact stars with Gaia Data Release 3 (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 on 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 5 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:

$\begin{eqnarray}&&{L}_{\mathrm{TESS}}=4\pi {d}^{2}\times {F}_{\mathrm{TESS}},\end{eqnarray} \tag{ 5 }$

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:

$\begin{eqnarray}&&{F}_{\mathrm{TESS}}={10}^{-0.4T}\times {F}_{0},\end{eqnarray} \tag{ 6 }$

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₀ = 4.03 × 10⁻⁶ erg s⁻¹ cm⁻², 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

$\begin{eqnarray}&&{E}_{\mathrm{TESS}}=\mathrm{ED}\times {L}_{\mathrm{TESS}},\end{eqnarray} \tag{ 7 }$

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%.

5.1.3. Fractional Flare Luminosity

For each flaring compact star, we calculate the fractional flare luminosity, represented as ${L}_{\mathrm{flare}}/{L}_{\mathrm{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 denoted as L_fl/L_Kp in Lurie et al. (2015) to characterize the level of flare activity 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, L_Kp has been substituted with the bolometric luminosity L_bol (see, e.g., Yang & Liu 2019; Ilin et al. 2021).

We calculate this ratio for a star using the following equation:

$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{flare}}}{{L}_{\mathrm{TESS}}}=\displaystyle \frac{\sum \,{\mathrm{ED}}_{i}}{{t}_{\mathrm{obs}}},\end{eqnarray} \tag{ 8 }$

where ED_i represents the ED of each flare of the star, and t_obs is the total observation time of the star by TESS. A higher value of this ratio signifies enhanced flare activity, indicative of an increased level of stellar magnetic activity.

We crossmatched the flaring compact stars with Gaia DR3 (Gaia Collaboration et al. 2016, 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 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.

Zoom In Zoom Out Reset image size

5.2. 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.

5.2.1. Evaluating the Level of Contamination

Due to the relatively low angular resolution of TESS (21'' pixel⁻¹), 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 tpfplotter ¹⁸ (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). ¹⁹ This enhancement allows for a more convenient assessment of the level of contamination 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 a greater number of communities, for instance, the exoplanet and variable star research communities. Our script is publicly available on Github. ²⁰

Zoom In Zoom Out Reset image size

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 contamination level (CL) for each star, as shown in Table 4. A CL = 0 indicates that the target is the only object within the aperture, without any significant contaminating photometry from nearby stars. A CL = 1 denotes that the target is the brightest object within the aperture, but there are other dim stars present. Lastly, a CL = 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 CL = 0, indicating an uncontaminated target star. Figure 9(b) displays minor contamination, warranting a CL = 1. Cases of severe contamination, as demonstrated by Figures 9(c) and (d), are assigned with a CL = 2.

In our categorization, 44 stars have a CL = 0, 62 stars fall into the CL = 1 category, and 87 stars are classified as CL = 2, which corresponds to 279, 256, and 481 flares, respectively. These levels of contamination 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 level of contamination does not necessarily disqualify a star as a flare candidate. However, stars with a high level of contamination warrant additional caution and follow-up investigation to confirm the source of the flare events.

5.2.2. Detecting Companion Stars

To address potential contamination from companion stars, we first exclude flaring stars that are labeled as binary in the TASC WG8 target list. We then locate the rest of the flaring stars on the Gaia H-R diagram. If a target falls on the MS, it indicates that the compact star has a brighter MS companion dominating the observed brightness, which is especially relevant for intrinsically faint WDs. Therefore, we exclude the flaring compact stars positioned within the MS of the Gaia H-R diagram (G_BP − G_RP > 1 and G_abs < 12, see Figure 8). We also exclude 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 ²¹ (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 speedyfit ²² to fit the SEDs of the sdB/sdO stars. Initially, we used spectral models exclusively from the Tübingen NLTE Model-Atmosphere Package (Werner & Dreizler 1999; Rauch & Deetjen 2003; Werner et al. 2003). However, the fitting residuals indicated a clear infrared 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–4500 K, consistent with K/M dwarfs. A more detailed description of the fitting process of speedyfit is given in Vos et al. (2017, 2018).

5.2.3. 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 dwarfs. Therefore, for WDs, there is a high probability that the flares originate from their red companions, especially given the comparable levels of luminosity 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 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 that a 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 a CL = 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 as hot subdwarfs or hot subdwarf candidates, but labeled as a hot subdwarf in the 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.

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

TIC	Gaia DR3 Source ID	Name	SpClass	T	GBP − GRP	Gabs	${N}_{\mathrm{flares}}$	Flare Freq.	$\mathrm{log}({L}_{\mathrm{flare}}/{L}_{\mathrm{TESS}})$
				(mag)	(mag)	(mag)		(sector−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
6997163	650908345319118720	GALEXJ08259+1307	sdB	14.203	0.0751	4.5317	1	0.5	−4.54
52078744	4704482467645621376	GALEXJ01077-6707	sdB	13.861	−0.2549	4.6422	1	0.167	−5.02
118327563	5000760581717433088	CD-38222	sdB	10.515	−0.319	4.1129	1	0.333	−5.84
150265769	781164326766404736	Feige34	sdO	11.5152	−0.4893	4.305	1	0.5	−5.47
167976324	5465148904077059968	GALEXJ10078-2924	sdO	12.779	−0.3275	4.0243	1	0.25	−5.48
202507151	1615596448547786496	PG1524+611	sdB+G0V	12.39	0.1957	3.5597	1	0.143	−6.27
206688085	6602969028791558016	GALEXJ22568-3308	sdB+F	13.269	0.4541	2.8338	1	0.25	−5.67
219974863	6468929937072637312	EC20217-5704	sdOB	12.4	−0.2958	3.9193	1	0.333	−4.84
231712886	6499440216512468096	EC23257-5443	sdB+F	13.568	1.1109	6.0068	3	1.0	−4.02
234281664	6404530097924835968	BPSCS22956-94	He-sdB	12.901	0.1036	4.5531	6	1.2	−3.76
262846506	2560685069816099968	PHL1079	sdB+G7V	12.901	0.1101	4.6389	1	0.333	−5.69
368628965	23330399790876032	PG0232+095	sdB+G1V	11.77	0.7273	3.2486	4	2.0	−3.95
443619867	3476266612927121536	EC12219-2618	sdB+MS	14.235	0.3135	3.9936	1	1.0	−4.12

Note. This table lists the 13 sdB/sdO stars in our refined sample. Column (1): TIC ID. Column (2): Gaia DR3 source ID. Column (3): target name. Column (4): spectral classification given by Culpan et al. (2022). Column (5): TESS magnitude from TIC v8.2 (Stassun et al. 2019; Paegert et al. 2021). Column (6): BP-RP color index from Gaia DR3. Column (7): absolute Gaia G magnitude from Gaia DR3. Column (8): number of observed flares. Column (9): frequency of flare occurrence per TESS sector. Column (10): logarithmic fractional flare luminosity.

Download table as:  ASCIITypeset image

5.3. 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. This can be expressed as

$\begin{eqnarray}&&{dN}\,(E)={{kE}}^{-\alpha }{dE},\end{eqnarray} \tag{ 9 }$

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 power-law 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

$\begin{eqnarray}&&\mathrm{log}\,(\nu )=\beta +(1-\alpha )\mathrm{log}\,(E),\end{eqnarray} \tag{ 10 }$

where $\beta =\mathrm{log}\,(\tfrac{k}{1-\alpha })$, and ν denotes the frequency of flares, i.e., the number of flares per unit of 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:

$\begin{eqnarray}&&\nu =\sum \,\displaystyle \frac{{N}_{i}\,(\gt E)}{{t}_{i}},\end{eqnarray} \tag{ 11 }$

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. Figure 10 displays the cumulative FFDs with corresponding best-fit lines, the fitted α values and uncertainties, and the number of flares included in the fitting.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Based on Cycles 1–5 of TESS photometry, we comprehensively investigated the flaring activity observed in sdB/sdO stars and WDs and identified 1016 flare events from 193 compact stars (Tables 3 and 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 filter out false positives effectively. 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 open-source 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 an analysis using Gaia DR3 data to identify binary systems with a bright MS companion, 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 of the light curves (only ∼25% flaring compact stars have a CL = 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 (approximately 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 that the observed flares may not be from the cool MS companions. We aim to explore this scenario further by 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 the lowest level of contamination from nearby objects (CL = 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 contamination 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 contamination by nearby active stars, as energy calculations use stellar luminosity (see Equation (7)), which is greater in sdB/sdO stars. The diminishment of this discrepancy for the CL = 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 a higher level of flare activity. Again, this potential contamination effect diminishes when analyzing the CL = 0 subset. Therefore, while definitive compact star flare attribution remains complex presently, comparing WD and sdB/sdO flare distributions provides initial insights into 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 of α ∼ 2 when fitted with the entire flare sample, which is polluted by flares originating from other late-type MS stars (F- to M-type stars). These MS stars, with FFDs having an index of α ∼ 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–M. They found that the FFD from A-type stars has an index of 1 < α < 1.5 when using Kepler photometry, a finding that markedly deviates from the α ∼ 2 observed in cooler F- to 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 the 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 late-type MS stars.

This observation prompted 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 show magnetic activity. 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. It is important to note that some previous work has found alternative explanations for putative flares on A-type stars (Pedersen et al. 2017).

We thus hypothesize that some of the stellar flares from the refined sdB/sdO sample may originate solely from the 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 in 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 that stellar flares on sdB/sdO stars may operate through mechanisms similar 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 that differential rotation might suffice to generate local magnetic fields in B/A-type 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 interest 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 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'' pixel⁻¹, 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 well integrated with Kepler/K2 photometry, which is useful for visualizing target stars and their surrounding environments. Beyond the research on stellar flares, the identification charts generated by tpfi could also be invaluable for other studies, for instance, exoplanet detection and variable star research, underscoring its broader astronomical applicability.

We acknowledge the support from the National Natural Science Foundation of China (NSFC) through grant Nos. 12273002, 12090040, 12090042, 11988101, and 11933004. This work is supported by the International Centre of Supernovae at Yunnan Key Laboratory (Nos. 202302AN360001 and 202302AN36000102) and the science research grants from the China Manned Space Project. R.S. acknowledges support from the INAF mini-grant on "Hot subdwarfs and white dwarfs: Pulsations, Binaries, and Planetary Systems." S.C. has financial support from the Centre National d'Etudes Spatiales (CNES, France). T.C. is supported by the LAMOST fellowship as a Youth Researcher, which is supported by the Special Funding for Advanced Users, budgeted and administrated by the Center for Astronomical Mega-Science, Chinese Academy of Sciences (CAMS), and acknowledges funding from the China Postdoctoral Science Foundation (2023M730297). T.W. acknowledges the support from the National Key Research and Development Program of China (grant No. 2021YFA1600402), the B-type Strategic Priority Program of the Chinese Academy of Sciences (grant No. XDB41000000), the Youth Innovation Promotion Association of the Chinese Academy of Sciences, and the Yunnan Ten Thousand Talents Plan Young & Elite Talents Project. All the TESS data used in this paper can be found in MAST (MAST Team 2021a, 2021b). The authors gratefully acknowledge the TESS team and all who have contributed to making this mission possible. Funding for the TESS mission is provided by the NASA Explorer Program.

Facility: TESS -

Software Astropy (Astropy Collaboration et al. 2013,2018, 2022), astroquery (Ginsburg et al. 2019), lightkurve (Lightkurve Collaboration et al. 1812), matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020).

We present all 23 flare events in the refined sample of 13 sdB/sdO stars in Figure 11. The flare event from TIC 234281664 in TESS Sector 28 is a complex flare event composed of two classical (single-peak) flares.