Correcting Stellar Flare Frequency Distributions Detected by TESS and Kepler

The habitability of planets is closely connected with the stellar activity, mainly the frequency of flares and the distribution of flare energy. Kepler and TESS find many flaring stars are detected via precise time-domain photometric data, and the frequency and energy distribution of stellar flares on different types of stars are studied statistically. However, the completeness and observational bias of detected flare events from different missions (e.g. Kepler and TESS) vary a lot. We use a unified data processing and detection method for flares events based on the light curve from Kepler and TESS. Then we perform injection and recovery tests in the original light curve of each star for each flare event to correct the completeness and energy of flares. Three samples of flaring stars are selected from Kepler and TESS, with rotating periods from 1 to $\sim$ 5 days. Adopting a hot-blackbody assumption, our results show that the cumulative flare frequency distributions (FFDs) of the same stars in Kepler and TESS bands tend to be consistent after correction, revealing a more natural flaring frequency and energy distribution. Our results also extend the low-energy limit in cumulative FFD fitting to $10^{31.5-33}$ erg on different types of stars. For solar-type stars, the average power-law index of cumulative FFD ($\alpha_{\rm cum}$) is $-0.84$, which indicates that low-energy flares contribute less to the total flare energy. With a piecewise correlation between $\alpha_{\rm cum}$ and $T_{\rm eff}$, $\alpha_{\rm cum}$ first rises with $T_{\rm eff}$ from M2 to K1 stars, then slightly decreases for stars hotter than K1.

Stellar flares are closely related to the stellar activity. The explosion of the stellar flares can release enormous energy, mainly caused by the magnetic reconnection process in the coronal region (Pettersen 1989;Benz & Güdel 2010;Shibata & Takasao 2016). Thus, detecting stellar flares give clues to understanding the stellar magnetic activity, which reflects the theory of stellar dynamo and internal structure. Stars of different spectral types have different convective depths. For example, low-mass stars (M * < 0.3 M ), i.e. stellar effective temperature cooler than M4, are fully convective, while more massive stars (0.3M < M * < 1.2M ) have a radiation layer inside the convective layer. According to Walkowicz et al. 2011 etc., flares on low-mass stars are usually more frequent.
Stellar flares usually release intense ultraviolet radiation and are crucial to understanding the habitability of exoplanets (Buccino et al. 2006). The far-ultraviolet radiation produced by flares can increase atmospheric erosion. Strong ultraviolet radiation reduces the habitability of rocky planets. Although Proxima Centauri b (Anglada-Escudé et al. 2016) is in the liquid-water habitable zone of its host star, the intensity of ultraviolet radiation on its surface is about 400 times that on the Earth, and there are many superflares (∼ 8/year, E bol >10 33 erg; Güdel 2004;Welsh et al. 2007; Abrevaya et al. 2020). Howard et al. (2018) even reported the first naked-eye superflare detected from Proxima Centauri. Thus, the habitability of Proxima Centauri b is still debatable (Ribas et al. 2016). Research on ultraviolet energy and flare frequency is essential for studying habitability, especially for the rocky planets around M dwarfs. However, a moderate near-ultraviolet flux may be necessary for prebiotic chemistry on the surface of rocky planets around low-mass stars (Buccino et al. 2007), e.g. TRAPPIST-1 planets (Gillon et al. 2017) and Proxima Centauri b .
To characterize the flare energy and frequency, the flare frequency distribution (FFD) is adopted and usually described by a power-law relation, i.e. dN/dE ∝ E α (Dennis 1985). The index α is utilized to constrain the magnetic activity for different stars. Yang & Liu (2019) presented a flare catalog of the Kepler mission using the long-cadence data in Data Release 25. They found that the FFDs from F-type stars to M-type stars obey this power-law relation with α ∼ -2, indicating that they have a similar mechanism on generating flares. But the deviation of FFDs on A-type stars with α ∼ -1.1 implies a different mechanism. However, other works do not support similar α for stars of different type. The α of solar flares is about -1.75 (e.g., Crosby et al. 1993;Shimizu 1995;Aschwanden et al. 2000). For latertype main-sequence stars, measured α values range from -2.0 to -1.4 (e.g., Paudel et al. 2018). The α value of -2.0 is critical. If α < -2, the total flare energy distribution is dominated by low-energy flare events (e.g. Güdel et al. 2003). In this case, the flares with low energy become important for redistributing the energy in the stellar atmosphere. For example, nanoflares with energies of ∼ 10 24 erg have been suggested to heat the quiescent solar corona (Parker 1988;Hudson 1991).
For the detection of flares, the current popular definition follows Chang et al. (2015) (e.g. Davenport 2016;Yang & Liu 2019, etc.). To eliminate false-positive events, the flares require three consecutive points fulfilling the same criteria. Although the selected flares ensure the accuracy of the study of the morphological characteristics of a flare (e.g., Davenport et al. 2014;Yun et al. 2017;Yan et al. 2021), many real flare events with few sampling points are eliminated due to light curve cadence. Therefore, the number of low-energy or short-duration flares is underestimated. In addition, due to the detrending of a light curve to remove variable shape caused by the stellar rotation, and the uncertainty in reconstructing the flare shapes, the uncertainties in calculation of flare energy are ∼ 60% − 65% (e.g., Yang & Liu 2019).
Since both the detection completeness of flare and the accuracy of energy calculation are important while fitting the FFDs. In this paper, we try to improve the completeness and the precision of energy calculation in detecting flares based on Kepler and TESS data, which can be used to correct the FFDs of different types of stars. Section 2 describes the data and samples chosen in this work. Section 3 describes the methods of data processing and FFD correction. Section 4 demonstrates our results based on our stellar samples. Finally, the summary and further discussion are presented in Section 5.

DATA AND SAMPLE SELECTION
In this paper, we use the photometry data from Kepler and TESS to study stellar flares because of their high photometric precision and continuity. Both Kepler and TESS provide large samples of stellar light curves, which is the basis for studying the stellar flares statistically.
Kepler was launched in 2009 with the primary scientific goal of exploring the structure and diversity of planetary systems by surveying a large sample of stars in a region of the Cygnus and Lyra constellations of our Galaxy. It is a Schmidt telescope design with a 0.95 m aperture and a spectral bandpass from 400 nm to 900 nm (Van Cleve & Caldwell 2016). The photometric precision was expected to be 20 ppm for a 12th magnitude G2V star with a 6.5 hr integration, though solar-type stars themselves turned out to be noisier than expected (Gilliland et al. 2011(Gilliland et al. , 2015. It has monitored the targets continuously and simultaneously, with a sampling interval of up to 1 minute for short-cadence (SC) data and 29.4 minutes for long-cadence (LC) data. Only 5000 targets have SC data, and their average time of observation is about two months. More than 200,000 targets have LC data, covering the whole Kepler mission for 48 months (from Q1 to Q17).
TESS was launched in 2018 with the primary goal of searching for transiting Earth-sized planets around nearby and bright stars. Its four 10 cm optical cameras simultaneously observe an entire field of 24°×96°, with a red-optical bandpass covering the wavelength range from about 600 to 1000 nm (Tenenbaum & Jenkins 2018). Note the TESS band is redder than the Kepler band. TESS measured light curves of over 200,000 preselected stars in its first two-year primary mission with a cadence of 2 minutes. Each sector of TESS covers about 27 days.
In 2019 July, TESS moved its field of view toward the northern ecliptic hemisphere for its second-year mission and observed the sky monitored by Kepler. Therefore, it allows us to statistically analyze the same flaring stars from both Kepler and TESS (e.g., Davenport et al. 2020). Note that the photometric precision of different missions and the contamination of nearby stars caused by the large pixel scale of TESS (see Section 3.2 and Section 4.1) mean that we need to deal with TESS data carefully.
In previous works, Yang & Liu (2019) (hereafter Yang19) presented a flare catalog of the Kepler mission using the LC data of Data Release 25, comprising 3420 flaring stars and 162,262 flare events. Günther et al. (2020) (hereafter Gün20) performed a study of stellar flares for the 24,809 stars observed with a cadence of 2 minutes during the first two months of the TESS mission. Tu et al. (2020) (hereafter Tu20) studied superflares on solar-type stars from the first year's Observations from TESS. We collected the flare events on solar-type stars from their catalog and plotted the distribution of flare energy in Figure 1(a). In a relatively complete detection range (flare energy between 10 34.5 and 10 36 erg), the power-law index (α) of the FFDs from Yang19, Tu20 and Gün20 is -2.45±0.06, -1.64±0.20, and -2.17±0.11 respectively (see Figure 1 (b) Figure 1. Different flare catalogs of solar-type stars. (a) The histogram of 15,056 flares on 507 stars from Yang19 (solid blue), 989 superflares on 307 stars from Tu20 (solid green), and 502 flares on 107 stars from Gün20 (solid red). The blue, green, and red dashed lines are the cumulative distribution functions from Yang19, Tu20 and Gün20, respectively. N f is the number of flares. (b) The blue, green, and red dots are the FFDs from Yang19, Tu20 and Gün20, respectively. The blue, green, and red solid lines are the fitted power laws of the respective FFDs in the energy range of 10 34.5 − 10 36.0 erg. The FFDs of different catalogs using are fitted the equation log10(dN (E)/(dE · dt)) = α · log10(E) + C (see Section 3.5). The text in the figure indicates the fitted parameter values and the errors of α.
When using a complete method of detecting flare in light curves from different missions, the white-light FFDs of the same spectral type stars obtained from different missions should be similar and closer to the actual FFDs, especially in the flare energy region with better detection completeness for each mission. Therefore, we selected our samples from these three catalogs of flaring stars. First, we choose the samples observed by both Kepler and TESS to correct the difference due to bandpass, data cadence, photometric precision etc.. As can be seen in Section 2.1, most targets are M-type stars. Then we add stars of other spectral types to extend the sample to investigate the flare statistical characteristics of flares for different stars (see Section 2.2).

Kepler Flaring Star with TESS Light Curve
Assuming that the actual FFD of a star is time-independent, we should obtain the same FFD on each star based on Kepler and TESS data. However, due to bandpass, photometric precision, and incompleteness of detection in different missions, the calculated number of flares and energy will lead to different FFDs via Kepler and TESS. We hope to find a method to correct the bias of FFDs in different missions. With the goal of correcting the FFDs to be closer to the actual values, we first choose the flaring stars observed by both Kepler and TESS missions.
We cross-matched Yang & Liu (2019)'s catalog with the TESS Target List 2 and got 193 targets within a distance of less than 1 arcsec. We sorted the target list according to the value of N f (the number of flares in the Kepler mission) and selected stars with N flare greater than 300 or stars with both N f between 100 and 300 and Tmag between 12.13 and 14.13 which have relatively better precision in light curves of TESS. Finally, 13 stars are selected-one G-type star, two K-type stars, and 10 M-type stars. The parameters of each star are listed in Table 1. These Kepler flaring stars with TESS light curves (hereafter denoted as Sample-1) have Tmag from 12.13 to 14.46 (see Table 1). These stars are faint for the TESS mission, and their photometric light curves from TESS are not as good as those from Kepler for detecting flares. For example, for the faintest star, TIC 27685142, a solar-type star with Tmag=14.46, we can only detect 10 flare events using the TESS data (see Section 3.2). We cannot fit the cumulative FFD of the TESS mission for stars with few detected flares. Therefore, we verify the criterion while detecting flares in Section 3.2, and cut the sample more strictly in Section 4.1 to exclude large uncertainties of the cumulative FFDs.

Select of Stars for the Completeness of Spectral Type
In Section 2.1, we have selected 13 stars with different spectral types. As the most interesting case, the statistics of flare events for the solar-type stars are important for understanding the quiescent Sun, but Sample-1 only contains one solar-type star. Note we adopt stars with effective temperature between 5100 K and 6000 K as solar-type stars. To involve more solar-type stars, we select solar-type stars with a rotation period between 1 and 5 days and randomly choose 84 targets in Tu et al. (2020) (hereafter denoted as Sample-2). The cut of rotation period is because the burst frequency of flares is related to the star's rotation period (Davenport 2016;Yang & Liu 2019;Tu et al. 2020, etc.). Sample 2 will be used to investigate the flares characteristics of solar-type stars.
Furthermore, to investigate the variation of a star's FFD with spectral type, we need to add more stars with different spectral types. We select the stars from the catalogs in the papers of Günther et al. (2020), Tu et al. (2020) and Yang & Liu (2019) (see Table 1). The stars are classified into K4-K2, K8-K5, M2-M0, and M5-M3 according to their effective temperature via Pecaut & Mamajek (2013). For each stellar type, we select the top five stars with the highest flare frequency in Günther et al. (2020)'s catalog, while the star with the highest flare frequency in Yang & Liu (2019)'s catalog is selected. In total, we select 24 stars as a sample including different spectral types (hereafter denoted as Sample-3).
In Sample-2 and Sample-3, the Tmag of TESS flaring stars varies from 7.78 to 12.14, while the Kmag of Kepler flaring stars in Sample-3 varies from 10.96 to 16.48. The light curves of these flaring stars are precise enough to detect flares (see Section 3.2).
To more accurately estimate flare energy based on stellar parameters, and to exclude contamination from binaries (e.g., Notsu et al. 2019), we cross-matched these stars to GAIA DR3 (Fouesneau et al. 2022). We also listed the updated stellar parameters and renormalised unit weight error (RUWE) from GAIA DR3 in Table 1 and applied them in Section 4.

METHOD
To correct the completeness of flare detection based on the data from Kepler and TESS, we reprocess the light curve uniformly and correct the detection efficiency by the following steps for each star: (i) detrend the light curve, (ii) find flare events, (iii) perform injection and recovery tests for each flare event, and (iv) correct the number and energy of flare events based on the detection probability and energy recovery ratio obtained from step (iii), respectively. Finally, (v) we fit the cumulative FFD to obtain the power-law parameters. Figure 2 demonstrates the procedures we adopted in this section, which will be introduced in detail one by one as follows.

Detrending the Light Curve
Flare detection and energy estimation are sensitive to the detrending method. In this subsection, we aim to correct the bias of the detrending methods in Kepler and TESS light curves, via a homogeneous detrending method.
As usual, we use the publicly available version of the Kepler or TESS photometry data, i.e. the pre-search data conditioning simple aperture photometry (PDCSAP) flux (Smith et al. 2012). The PDCSAP fluxes are obtained from a systematic reprocessing of the entire database using a Bayesian approach to remove systematic noises in both the short-and long-term light curves .
First, we smooth out the secular modulation signal. Note, the sampling stars are from the catalog of Kepler and TESS flaring stars (Yang & Liu 2019;Tu et al. 2020;Günther et al. 2020), which has excluded out pulsating stars and binaries. Thus, the secular modulation signal is mainly caused by the spots rotating with star.
We detrend the light curve using the Savitzky-Golay filter (Savitzky & Golay 1964) (see function scipy.signal.savgol filter function in SciPy package (Virtanen et al. 2020). The Savitzky-Golay filter is based on local polynomial least-squares fitting in the time series and can filter out noise while ensuring that the shape and width of the signal remain unchanged.
There are two crucial parameters-window length and polyorder-when we use the Savitzky-Golay filter to detrend light curves. If the window length is too small, it will smooth out the characteristics of light variation of flares, especially long-lived ones. Conversely, if the window length is too large, it will not be able to smooth the rotation modulation signal for fast rotating stars. The optimal window length depends on the photometric error, sampling rate, light-curve trend, and filter order (Sadeghi & Behnia 2018). Thus, calculating the optimal window length requires a prior knowledge of the light-curve trend caused by stellar rotation modulation.
Yang & Liu (2019) used the LombScargle periodogram (Lomb 1976;Scargle 1982) to calculate the most significant period of each block, then smooth each block with the filter width of one-fourth of the most significant period. Their lower limit on filter width is 6 hr and the upper limit is 24 hr. We also set the window length as 1/4 of the stellar rotation period, while the stellar rotation period is from the previous catalog calculated through the LombScargle periodogram. Since all the stellar rotation periods in our samples are between 1 to 5 days, the lower and upper limits on window length are 6 hr and 30 hr, respectively. We set the polynomial order as 3 to filter the PDCSAP light curve, as in the previous work (e.g. Davenport 2016).
The detrending method adopted in this paper differs from previous works of finding flaring stars. Günther et al. (2020) computed the Lomb-Scargle periodogram and removed the periodic signal from light curve using a sine wave fit. Tu et al. (2020) detected superflares without detrending the light curves, while a quadratic function is adopted to remove the long-term stellar variability to calculate the energy of superflares.

Finding Flares and Energy Calculation
After detrending the light curve, we need to choose criteria to identify flare events and get the parameters of each flare. In this section, we will introduce our criteria for detecting flares for Kepler and TESS data. Furthermore, since the flare energy is crucial and is model-dependent, we introduce our method to calculate the flare energy based on the blackbody assumption.

Flare Detection
We identify flares based on the three parameters N 1,2,3 suggested by Chang et al. (2015) (the definitions of N 1,2,3 can be found in Appendix A). Chang et al. (2015) recommend that N 1,2,3 should be at least larger than 3, 1, and 3, respectively. We focus on the completeness of flare detection in this work rather than a detailed study of flare morphology. Therefore, we do not verify false positives visually for detected flare events.
To determine the parameters N 1,2,3 , which are related to the photometric precision, we estimate the combined differential photometric precision (CDPP) using the Savitzky-Golay method (see the estimate cdpp function in the Lightkurve package; (Lightkurve Collaboration et al. 2018)). For these 13 Kepler flaring stars with TESS light curve with Tmag from 12.13 to 14.46, the CDPPs of the Kepler light curve are around 200 ppm, and then the N 1,2,3 values are set as 3, 2, and 3, respectively. The CDPPs of the TESS mission for these targets are 1470-9480 ppm. To identify more flares with such photometric accuracy of TESS, the N 1,2 are set as 2, and N 3 is set as 2 to detect flares with shorter duration, when we detect flares from TESS light curves. Davenport (2016) also found that the flare recovery in Kepler is improved via long-cadence data and that adjusting the scale factor N 3 from 3 to 2 did not enhance the recovered flares obviously for short-cadence data. We also visually inspect the choice of parameter N 1,2,3 , especially for the light curve of TESS, and believe that the flare events we found are confirmed.
After we set the detection criteria of N 1,2,3 , some flares can be detected at the beginning or end of each segment. These flares are suspected when such events are only partially observed or correlated with the detrending precision at the end or beginning of available light-curve data (Chang et al. 2015). Therefore, we discard flare events that overlapped with the 4 data points (0.005 days) and six data points (0.125 days), at both beginnings and ends of each light-curve segment, for the TESS and Kepler data, respectively. Additionally, there are impulsive outliers in the data marked as '512' (Tenenbaum & Jenkins 2018). We also abandon the flares that overlap with these outliers, which are usually doubted. Consequently, the available observation time is shortened in the FFD calculation. Take TESS data in Sector 1 as an instance: after excluding both the beginning and end of each segment in Sector 1, as well as the impulsive outliers, the available observational time is reduced from 27.88 to 25.36 days. We use the reduced observational time to calculate the FFD of each star.
We compare our flare events detected with those in previous works to validate the choice of N 1,2,3 adopted in this paper. We choose TIC 206592394 because both Günther et al. (2020) and Tu et al. (2020)

Flare Energy Calculation
According to previous work, the bolometric energy of a flare can hard be calculate precisely, because of the complexity of the flare generation mechanism. The flares have different time-dependent spectra, even if flares on the same star. The flare energies may depend on the area and magnetic structure of starspots (Okamoto et al. 2021) during different flare events. At present, a common method for calculating the flare energy is to assume that the hot-blackbody radiation with a constant temperature of 9000K or 10,000K in the active area of the stellar surface, and that only the active area changes during the flares. The assumption is not correct, since the temperature of the stellar flares will change significantly during the injection and release of energy. Furthermore, the hot-blackbody assumption cannot well represent the Balmer jump and underestimates the energy radiated in the near-UV region by factors of 2-3 (Kowalski et al. 2016(Kowalski et al. , 2018(Kowalski et al. , 2019. Accurate estimation of flare energy requires more observations of multiple wave bands (e.g., France et al. 2016;Richey-Yowell et al. 2022) and observation of time series for each event.
For comparison with previous work, we simply adopt the hot-blackbody assumption in this paper. The FFDs from F to M stars all obey a power-law relation with α ∼ 2, indicating that they have the similar burst mechanism of magnetic reconnection (e.g., Yang & Liu 2019;Günther et al. 2020;Jackman et al. 2021). Then, we assume that the same spectral shape among flares on different stars with different flare energies, i.e. temperatures of flares are all set as 9000 K (Shibayama et al. 2013;Davenport 2016;Yang et al. 2018;Yang & Liu 2019;Günther et al. 2020). There is another method to calculate the energy of flares, e.g. Wu et al. (2015) and Tu et al. (2020) integrate the energy based on the enhancement in the light curve, assuming the enhancement is constant in all wavelengths. Refer to Appendix B for exact energy calculation formulae and differences between the two methods.
After finding flare events, the pipeline records the amplitude (Amp.), equivalent duration (ED) and duration (Dur.) of flares (Davenport 2016;Ilin et al. 2021). These parameters can be used to calculate the flare energy in the observational band. Since the energy calculation depends on the parameters ED, the uncertainty is sensitive to photometric precision, i.e. the recorded error of ED is from the uncertainty on the flare flux values (see Davenport 2016 Eqn. (2) for details). For example, TIC 312590891 in Sample-2 with Tmag=9.45 (CDPP=255 ppm) has a median error δ ED = 0.076 s for all its flares found in TESS data, while a fainter star, TIC 27454242 in Sample-1 with Tmag=13.15 (CDPP=2554 ppm), has a larger median error δ ED = 0.264 s due to the 10 times larger photometric uncertainty. Additionally, TIC 27454242 is also observed by Kepler (KIC 12314646) with a much better photometric precision of 159 ppm. Thus, the δ ED from Kepler data is reduced to 0.049 seconds. Obviously, for the flare events with the same ED, the worse photometric precision causes larger δ ED , or the larger energy uncertainty.
Besides the photometric error, the uncertainties of the stellar parameters (e.g. T eff , R * ) also contribute to the energy uncertainty. Thus the calculated energy may be not precise for individual flare events. However, we focus on the statistical results of FFDs in this paper, especially the slope of FFDs, which is not sensitive to a consistent deviation due to uncertainties of stellar properties. Therefore we do not consider the uncertainty of calculated flare energy .

Validation of Flare Detection and Energy Calculation via Our Method
To verify the completeness of flare detection, especially for fast rotating stars with various periodic light-curve amplitudes (e.g., Okamoto et al. 2021), we analyze the distribution of flare parameters with stellar parameters. The detrending light curve procedure is related to the stellar rotation period and the periodic amplitude of the light curve (LCAmp). The Tmag of the star indicates the photometric accuracy, which is directly related to the procedure for finding flares. Stellar rotation periods are taken from Günther et al. (2020) for TESS flaring stars, and from Reinhold et al. (2013) and McQuillan et al. (2014) for Kepler flaring stars. We obtain the LCAmp following McQuillan et al. (2013). Panel (a-c) of Figure. 4 shows the distribution of the normalized flare amplitude with the stellar parameters (Tmag, P rot , LCAmp) in Sample-2. Obviously, the normalized amplitudes of detected flares are positively correlated with the T-mag, but they have no obvious correlation with P rot and LCAmp. The results indicate that the detrending method remove the variation module in light curves well.  In order to validate the flare detection and energy calculation via our method, we compare the detected flare events in Sample-2 and Sample-3 with previous work. For this comparison, we use the same stellar parameters as in the previous work. In Section 4, we used the latest stellar parameters from the GAIA DR3 catalog.
We detect 792 flare events on Sample-2 stars, including ∼95% of superflares in the catalog of Tu et al. (2020). i.e. 169 flare events are cross-matched with the previous catalog of superflares, which includes 178 events. Figure 5(a) shows the histogram and cumulative distribution function (CDF) of flares in Tu et al. (2020) and all flares we detected. More flares with energy below 10 34 erg are detected in our methods.
Note that we calculate the energy of flare using the assumption of the hot-blackbody radiation, while Tu et al. (2020) integrate the energy based on the enhancement in the light curve (see Appendix B for detail about the two ways to calculate flare energy). In order to compare the similarity of flares between this work and Tu et al. (2020), we correct the energy of flare events we detected via a factor r in Eqn. B9 to correct the difference due to the two calculation methods. Figure 5(b) shows the histogram in Tu et al. (2020) (green line) and the cross-matched flares we detected (red line). We also use a KolmogorovSmirnov (K-S) test to test the similarity of the two distributions, and this returns a P-value of 0.77, proving that our detection of superflares has nearly the same completeness and energy distribution as Tu et al. (2020).   (Figure 6(b)).
To confirm the reliability of our flare detection and energy calculation, we also compare the flare distribution as a function of energy detected by our pipeline with that detected by Günther et al. (2020), as shown in Figure 6 ours returns a P-value of 0.63, which shows that our detection completeness and energy calculation are consistent with Günther et al. (2020).
We conclude that the flare events we detected with higher energy are relatively consistent with previous works in both Sample-2 and Sample-3, which validates the detection method adopted by this work. Furthermore, we detect more flares with lower energies and can improve the completeness of flares to constrain the slope of the distribution with energy (α, see Section 3.5).

Inject and recovery test
As seen in section 3.2, flares with higher energy are detected more completely in different works. However, flares with lower energies are not completely detected. Therefore, we perform injection and recovery tests to correct the completeness via the detection probability and improve the energy estimation via the energy recovery ratio.
We synthesize flare shape with exactly the same amplitude and duration of each flare we detected in all samples, using the empirical model of Davenport et al. (2014). Then, we inject synthesized flare events at random epochs in the original light curve of the corresponding star. To avoid the overlapping of flare events, we limit the injection frequency to 0.5 per day, and make sure each segment has at least one injection. The flares are set individually and do not overlap with each other. We generate ∼2000 injections for each flare from TESS, and ∼4000 to 5000 injections for each Kepler flare. Then, we detrend the light curve with injected flares and find flares using the methods introduced in Section 3.1 and Section 3.2. Based on the simulation of injection and recovery, the detection probability of each flare is defined as the fraction of recovered flare events, while the energy recovery ratio of each flare is the average of the ratio of recovered energy to the energy we injected initially. Finally, we obtain the detection probability and energy recovery ratio for each flare event based on a statistical analysis of the synthetic injection and recovery data (see the characterize flares, using AltaiPony package (Davenport 2016;Ilin et al. 2021).
We analyze the detection probability and energy recovery ratio of 84 targets in Tu et al. (2020) (Sample-2). The total 792 flare events we detected on these 84 TESS flaring stars returned a detection probability and energy recovery ratio for our injection and recovery test. Figure 7 (a) and (b) show the histogram and CDF of detection probability and energy recovery ratio. 50% of flare events have a detection probability of 0.8 or more, while 30% of flares have a detection probability of less than 0.2. Thus the number of flares detected with low detection probability is incomplete and needs to be corrected to estimate the real frequency. Furthermore, 60% of flare events have an energy recovery ratio less than 0.625, i.e. their energies are underestimated by at least 37.5%. Note also that 9% of flare events that have an energy recovery ratio of more than 1.0, which means their total energies are overestimated. Therefore, the recovered energy probably deviates from the real value probably and needs to be corrected according to the energy recovery ratio.
Usually, lower-energy flares lead to lower detection probabilities and larger energy uncertainties due to lower signalto-noise ratio. As seen in Figure 8

Correction for Each Flare Event
We denote the detection probability of a flare event as P det and the energy recovery ratio as R rec . Then the actual number of flare events should be enlarged by a factor of 1/P det . Assuming a flare event with an energy of E recovery , the actual energy should be corrected as E correct = E recovery /R rec . These corrections of numbers of flares and energies are adopted hereafter.
Note that if the detection probability or the energy recovery ratio is too small, the correction factor will return a large and unreliable value. Thus to avoid these extreme cases, we set a lower limit of 0.01. I.e. if the calculated probability or recovery ratio of a flare is less than 0.01, we skip over the correction to keep its energy and number as the original values and leave them as incomplete flares. However, only five flares (0.6%) in Sample-2 have such small P det and the minimum R rec is 0.3.
Finally, we make corrections according to each flare event's recovery probability and energy recovery ratio, flare by flare, star by star.

Fitting the flare frequency distribution
Flares with a power-law distribution in energy (Lacy et al. 1976) are described by: where N is the number of flares that occur in a given observational period dt, E is the flare energy, dt is the time spanning of observation, k is a constant of proportionality, and α is the power-law index.
One can obtain the standard cumulative FFD by integrating the number of flares in unit time with energy range greater than a certain value. i.e. the frequency of flares (denoted by ν) with energy ≥ E per day should be log 10 (ν) = β + α cum log 10 (E), where β = log 10 ( k 1+α ) and α cum =α+1 Jackman et al. 2021). The power-law index α cum represents the slope of the cumulative FFD. If a star has α cum < -1 ( i.e. α< -2), the low-energy flares contribute the majority of total energy emitted by flares. Conversely, on the stars with α cum > -1 (i.e. α> -2), the high-energy flares dominate the total energy (Paudel et al. 2018;Jackman et al. 2021).
Note that the detected FFD is not exactly a power-low distribution, and the fitted power-law indices do not sarisfy the equation of α cum =α+1 (Maschberger & Kroupa 2009). Thus, we cannot directly transfer α cum to α. In Section 4, we will fit either α cum or α consistently with previous works for comparison.
Since the fitted power-law index depends on the energy range we choose, the choice of energy range should be complete for flare detection. As described in Section 3.4, we correct the completeness of flares via the detection probability P det . For a given flare star, we can plot P det of each flare (e.g. Figure 8 (a)) and obtain a smoothed correlation curve via a Savitzky-Golay filter. We set two lower limits of P det , i.e. 0.65 and 0.15, to preserve the flare samples with high completeness and with moderate completeness. Then the number and energy of flares in three samples are corrected to fit the FFD.
We take KIC 7131515 as an instance. As shown in Figure 8 (c), the number of flares we detected after correction is obviously more than that before correction, and the energy extends as low as 10 32.5 erg. For flares with energy ≥ 10 34 erg, the histograms of corrected and original flares are similar, while those with energy 10 32.5−34.0 erg are very different.
The power-law index of cumulative FFDs of each star before and after correction is listed in Table 2 (Sample-1) and Table 3 (Sample-2 and Sample-3). The electronic version of Table 1-3 can be found on the website 3 . We will demonstrate the results in the next section. Note-"RangeAll" means the fitting within the whole energy range; "Range65" means the fitting within the energy range of P det ≥ 65%; "Range15" means the fitting within the energy range of P det ≥ 15%; a Because of the requirement of sufficient energy bins, we only use the whole energy range to fit the cumulative FFDs obtained from the TESS data for Sample-1 stars. The fitting of FFD fails for "Range65" and "Range15" due to the limited number of energy bins. For faint stars in Sample-1, considering their worse photometric precision from TESS, we correct the flares with P det ≥ 0.001 to enlarge the numbers of flares. (see Section 3.4) (This table is available in its entirety in machine-readable form.)

Kepler Flaring Stars with TESS Light Curves
For the Kepler flaring stars cross-matched with TESS Target List, we selected 13 targets (Sample-1) to conduct injection and recovery tests and correct every event using the detection probability and energy recovery ratio. We perform a more refined source screening to compare α cum obtained from Kepler and TESS flares, respectively.
Many fewer flare events are detected by the TESS light curve than by Kepler due to precision and the length of observation. Thus some stars only include a few flares with limited energy ranges. TIC 27685142 (KIC 11872364) does not have enough flare events to fit a cumulative FFD. Therefore, we excluded it from the subsequent analysis.
To avoid contamination by blended flaring stars, we also exclude stars with nearby sources. For the Kepler mission, about 10% of flares are probably polluted by nearby stars because some sources are very close to each other on CCDs (Shibayama et al. 2013;Gao et al. 2016). TESS has a larger pixel scale and causes more contamination by blended flaring stars. We cross-matched the Yang & Liu (2019)'s catalog with the TESS Target List using a distance of 21" (the TESS pixel scale) to check for contamination by nearby stars. In these 13 targets we selected, we found that three have nearby stars within 21", i.e. TIC 120629872 (KIC 2692704), TIC 406952118 (KIC 12646841) and TIC 120637232 (KIC 2692708). In the analysis of cumulative FFDs, we also exclude these three stars. Thus, there are six stars left for further analysis. To extend the energy range of FFDs, we use the whole energy range to fit α cum,T ESS in both original and corrected cases, while the energy range to fit α cum,Kepler is chosen as P det > 15%.
As a typical case, the cumulative FFDs of TIC 7454242 (KIC 12314646) from the TESS and Kepler mission are shown in Figure 9. The corrected FFDs are much closer than the original FFDs, which indicates that the corrections adopted in this paper are useful for correcting the detecting bias from Kepler and TESS data. Original α cum =-1.031±0.044 (d) For the remaining nine stars, Figure 10 shows a comparison of the α cum values of TESS and Kepler missions obtained before and after correction. Nearly all original α cum obtained from TESS data are larger than those from Kepler, which indicates the obvious incomplete detection of flares with low energies via TESS. Additionally, TESS tends to detect flares with larger energy due to its redder observational band than Kepler's (Doyle et al. 2019;Tu et al. 2021). Before correcting for flare events, we get an average difference of -0.33 between the two missions (∆α cum = α cum,Keplerα cum,TESS ), while this decreases to 0.04 after correction. This indicates that the slope of cumulative FFDs obtained from TESS and Kepler missions is consistent.
After the correction, the α cum of some targets obtained by TESS and Kepler data are still different, e.g. TIC 27010191, TIC 137759349 and TIC 138296722. We speculate that the FFDs have changed during different observational epochs from Kepler and TESS. The open green square is the original cumulative FFD detected from TESS data by our pipeline, and the green dot is the cumulative FFD after correction for each flare event using the detection probability and energy recovery ratio. The dashed-dotted and solid green lines indicate the fitted lines before and after correction, respectively. Similarly, for flare events detected from Kepler data, we use blue to present FFDs before and after correction.

Solar-type Stars
For the 84 solar-type flaring stars from Tu et al. (2020)'s catalog (see Section 2.2), 32, 53, and 59 stars have enough flare events to fit cumulative FFDs individually when the fitting energy range is set to P det ≥ 65%, ≥ 15% and the whole energy range, respectively. We fit these α cum in three different energy ranges, as shown in Figure 11, and denote the three setting groups with different P det as A, B, and C, respectively.
The average α cum in group A and B are nearly the same with a tiny difference of 0.024 (see the top panel of Figure  11), which indicates that the completeness is improved in the energy range P det > 15%, via the correction method in this paper. The average α cum in group C is ∼ 0.2 larger than that in B or A, which means the flares with lower energy are still incomplete and only partially corrected by our method. The open olive squares and solid red dots indicate differences of six stars, before and after correction, respectively. The red and olive horizontal lines indicate the average values of the difference before and after correction, which are -0.33 and 0.04, respectively. Note that because the number of flare events detected by the TESS light curve is small, we fit the cumulative FFD in the whole energy range for Sample-1 stars obtained from the TESS light curve, while the original α cum,Kepler come from the fitting of the original cumulative FFD with the energy range of P det ≥ 65%, and the corrected α cum,Kepler come from the fitting of the corrected cumulative FFD with the energy range of P det ≥ 15%.
Based on the updated stellar parameters from the GAIA DR3 catalog, we removed the stars with effective temperatures not in the range of 5100−6000K, and possible binaries (RUWE>1.4) in the analysis of solar-type stars. We collect the flare events on the remaining 53 stars to obtain a typical α cum of a solar-type star. α cum before and after correction is -1.24±0.10 and -1.08±0.04(P det > 15%), respectively.
In previous works, both Yang & Liu (2019) and Tu et al. (2020) fit the FFD in the energy range above 1034 erg, although Tu et al. (2020) only focus on superflares. In this paper, the original detected number of flares peaks at ∼ 10 33.5 erg. The lower energy limit with P det = 15% is extended to 10 33.0 erg. Thus our results of α cum are available in lower energy ranges.
To compare with the result for the Sun and other previous works, we also fit the α of FFDs of all flare events on the remaining 53 stars in Sample-2. I.e. α = −1.93 ± 0.12 with an energy range P det ≥ 65% before correction, and α = −2.00 ± 0.06 with an energy range P det ≥ 15% after correction.
Fitting α due to all flares collected on a group of stars assumes different weights for each star based on flare counts. It may be dominated by stars with many more flares. Even if all stars have the same α, the stacked FFD for all flares on these stars can have a different slope，because of the different counts and detectable energy range of FFDs for different stars. Thus, we also calculate the averaged α = −1.55 ± 0.29 (α cum = −0.84 ± 0.24), assuming equal weights for each star. The results is larger than the fitted value when combining all flares, but is consistent within 2 σ uncertainty.

Stars of Different Spectral Type
Unlike solar-type stars, stars with lower temperatures, especially M dwarfs, are usually more active. In this section, according to the updated stellar parameters from GAIA DR3, we collect 22 stars in Sample-3 and eight stars in Sample-1 (Section 4.1), as well as 45 stars in Sample-2. Note that the total of 75 flaring stars from TESS and Kepler have stellar rotating period in the range 1 to 5 days, with only two exceptions-TIC 121785583 (KIC 4355503) and TIC 120257064 (KIC 3934090) from Sample-1 (with the stellar rotation periods of 5.45 and 5.82 days).
To compare with different previous works, we calculate the average value of α cum from Günther et al. (2020)'s catalog for each stellar type and take the standard deviation of each type as the uncertainty. Jackman et al. (2021) gives the FFD fits of main-sequence stars with the Next Generation Transit Survey (NGTS, Wheatley et al. 2018), providing a positive correlation between α cum and the stellar types later than K5. The values from previous works are listed in Table 4 and shown in Figure 12. Note-Nstar means the number of flaring stars; "-" indicates no available value. a We calculate the average αcum in each subset and take the standard deviation as its error. b Jackman et al. (2021) fitted an average α combining flares on all stars, and adopted αcum value as α + 1, based on the assumption that all stars taken together were equivalent to a single, average, star (see the reference for detail). K8-K5 K4-K2 Solar-type Fitting in this work, M3-G0 cum from Gun20, M5-G0 corrected cum in this work average cum in this work, M5-G0 average cum from Gun20, M5-G0 average cum from Gun20, < M4 average cum from Jackman21, M5-K2 3) fitting the distribution relationship of αcum with T eff from M2 to G1 stars. The Pearson correlation between αcum and T eff is R=0.63 and R=-0.27 for stars from M2 to K1 and from K0 to G1, respectively. The coefficient of determination for fitting linear formulae are 0.39 and 0.10 in these two segments, respectively.
The average α cum of the targets from Günther et al. (2020) based on TESS is independent of spectral type with Pearson correlation R = −0.074. Their average α cum on K8-K5 or M2-M0 stars is consistent with those from Jackman et al. (2021). Note that both previous works do not correct the flare events via injection and recovery. After correction, our work provides α cum with a more obvious correlation to the stellar type for stars cooler than K1. The Pearson correlation between α cum and T eff is R=0.63, if we fit the positive trend as shown in Figure 12.
For solar-type stars, there is also a negative trend of α cum , which is also revealed in Table 4 of Okamoto et al. (2021)) for fast rotating solar-type stars. Therefore, we use a piecewise linear formula to fit the trends of α cum , in the range 3400-5900 K. The fitted linear formula is as follows: α cum = 0.44(±0.01) * T eff 1000 K − 2.84(±0.05) for 3400 K ≤ T eff < 5150 K −0.53(±0.03) * T eff 1000 K + 2.18(±0.18) for 5150 K ≤ T eff ≤ 5900 K Note that the coefficients of determination are 0.39 and 0.10 in the different segments, respectively. Thus, more data from solar-types stars with well-determined parameters are crucial for validating the negative trend of α cum for solar-type stars.
According to the fitted correlation, α cum increases from -1.34 for an M2 star (T eff =3400K) to -0.57 for a K1 star (T eff =5150 K). This predicts a higher frequency of lower-energy flares in cooler stars with thicker convective envelopes. Our result is consistent with the α − Ω dynamo theory, in which the magnetic field depends on the depth of the convective envelope (e.g., Moffatt 1978;Parker 1979). Conversely, for stars hotter than K1 (T eff ∼5150 K), α cum starts to decrease with T eff instead. For a solar-type star with T eff = 5550K, Equation 3 provides α cum = −0.79, which is consistent with the average value of α cum = −0.84 ± 0.24 in Sample-2. For stars with T eff in the range between 4170 and 5940 K, the calculated α cum >-1.0 according to Eqn. 3, thus the total energy is mostly contributed by flares with higher energy.
According to the updated stellar parameters from GAIA DR3, there are only two fully convective stars (cooler than M4 type) in our samples, with average α cum =-0.95. Eight M3 stars in our sample also show a similar α cum with an average value of -0.90. There are 26 fully convective stars with a stellar period of 1.0-5.0 days in Günther et al. (2020)'s catalog, and the average α cum =-0.80 is represented by the open gray square to the left of the vertical dashed-dotted line in Figure 12. It is a little larger than the average α cum of M0-M2 type stars from both Günther et al. (2020), Jackman et al. (2021) and this work. This indicates that the flare burst mechanism of fully convective stars may differ from other stars. The geometry of their magnetic fields is different and may have a switch in the dynamo (e.g., Stassun et al. 2011;Shulyak et al. 2017). That is why we exclude the 10 stars cooler than M2 in our sample to obtain the Equation 3.

SUMMARY AND DISCUSSION
Based on Kepler and TESS data, we consider the flare detection bias due to different missions and focus on FFD correction for flaring stars. We choose three samples with stellar rotation period ranging from 1 to ∼5 days in Section 2 according to the flaring catalogs from Günther et al. (2020), Tu et al. (2020), and Yang & Liu (2019).
We built an uniform process to detect flares and to improve the parameters of flare events in Section 3. Specifically, we select a suitable window width related to the stellar rotation period and detrend the light curve before finding flares in Section 3.1. Parametric thresholds are adopted due to photometric precision, and we calculate the flare energy by assuming a blackbody spectrum, and compare the completeness with previous works in Section 3.2. Each flare events are corrected by the detection probability P det and energy recovery ratio R rec in Section 3.4, which are obtained from an injection and recovery test in Section 3.3. Finally, we can fit the power-law index via the corrected cumulative FFD in Section 3.5, down to a lower energy limit with P det = 15%.
In Section 4.1, we analyze 13 Kepler flaring stars with TESS light curves in Sample-1, and select nine stars with reliable FFDs to fit α cum from TESS and Kepler data, respectively. Table 2 and Figure 10 show that original α cum from TESS are larger than those from Kepler because of the less complete detection of low-energy flares for TESS. The corrected α cum for most stars from Kepler and TESS become closer, as shown in Figure 10. We conclude that the slope of FFD or cumulative FFD (i.e., α or α cum ) after correction reveals more realistic values, and flares from different missions can be analyzed via the correction method to remove the detection bias.
We extend a large number of low-energy flares through the detection and correction for 84 solar-type flaring stars from TESS in Section 4.2. Thus we can reach a low-energy limit of 10 33 erg for cumulative FFD fitting on solar-type stars. Since the fitted α cum in the energy range P det ≥ 65%, P det ≥ 15% and the whole energy range are similar, we conclude that the flares are nearly complete in the energy range P det ≥ 15% as shown in Figure 11. Collecting all flares on 53 solar-type stars updated with stellar parameters from GAIA DR3, we obtain α cum ∼ -1.24±0.10 and -1.08±0.04 before and after correction, respectively, while the average α cum of all stars in Sample-2 is −0.84 ± 0.24. The later value is preferred in this paper (see Section 4.2), thus low-energy flares contribute less to the total flare energy on fast rotating solar-type stars.
In Section 4.3, based on stars in three samples, we obtain the α cum of stars from M5 to G1 through the correction of cumulative FFDs. Compared with previous works, we found that the α cum positively correlates to the stellar effective temperature for stars between M2 and K1. For stars hotter than K1, α cum may start to decrease with T eff slightly. A piecewise linear empirical formula is obtained to monitor the correlations as shown in Figure 12 and Eqn. 3. Assuming the spectrum is the same for all flares, for stars hotter than K5 stars and cooler than G0 stars, the total energy is mostly contributed by flares with higher energy. Fully convective stars cooler than M3 show a larger α cum than other M stars, which might indicate different dynamo mechanisms.
α cum is an excellent index to describe the FFD of a star. Coincidentally, one star (TIC 158552258 or KIC 7350067) of our Sample-1 has a confirmed super-Earth: Kepler-1646 b (Morton et al. 2016). The flaring host star (M dwarf) needs to be investigated more carefully to study the influence of flares on the super-Earth, and thus help us to know the location of the "abiogenesis zone" proposed by Rimmer et al. 2018;Günther et al. 2020.
In this work, we only study the fast rotating stars, because the rotation period is cut to 1 to ∼5 days. The rotation period depends on stellar age, which is also a crucial factor in determining stellar activity. As Wright et al. (2011) and Yang & Liu (2019) etc. show, flare activity of stars has a broken power law at the transitional Rossby number around 0.13. The Rossby number is the ratio of the stellar rotation period to the convective turnover time. Thus this is the reason why there is a large scatter of α cum for stars of different types. To obtain better control samples, flaring stars need to be selected more carefully due to various parameters, if a much larger population of flaring stars is obtained in the near future.
Additionally, the flare energy is calculated via simple models, i.e., the assumption of black-body radiation or a flat spectrum. Since Kepler and TESS have similar observational bands, thus we can correct the difference in FFD very well. However, it does not always work, especially for missions in quite different bands. Previous studies only discuss spectra from several early/mid M-dwarf flare stars (e.g., Kowalski et al. 2016Kowalski et al. , 2018Kowalski et al. , 2019. So far, we do not have enough knowledge of the spectra of G-dwarf flares. Time series spectral observation and better physical models of flares will benefit us in understanding the precise flare energy distribution at different wavelengths. And we can combine the time series data to model the flare evolution more precisely, and achieve more accurate energy estimation for these flares.
2. The positive excess plus the photometric error w i at epoch i, must be greater than N 2 σ L , i.e.
Approximately, the photometric error (w i ) and the local standard deviation (σ L ) can be considerated to be the same. The Equation A2 can be written as Therefore, when we set N 2 ≥ N 1 + 1, the criterion on N 1 is satisfied automatically. 3. To verify the confidence, the number of data points satisfying the above criteria must be at least N 3 continuously.

B. TWO METHODS TO CALCULATE FLARE ENERGY
There are two ways to calculate the flare energy. In this paper, we use the assumption of hot-blackbody radiation to calculate the flare energy (e.g., Shibayama et al. 2013;Davenport 2016;Yang & Liu 2019;Günther et al. 2020, etc.).
The luminosity of the star in the observing bandpass is The luminosity of flare in the observing bandpass is where R * is the stellar radius and R λ is the response function in the observational band (TESS or Kepler band). A flare (t) is the flare area while assuming that the flare temperature is constant throughout the duration. B λ (T flare ) and B λ (T eff ) are the Planck functions evaluated for the effective temperature of the flare and star, respectively.
The enhancement of the normalized light curve during flare events is which we can obtain from the observed light curves. We denote R λ B λ (T eff ) R λ B λ (T flare ) as η, which can be calculated through the response function of the observing bandpass and the effective temperatures of both star and flare.
Finally, we arrive at the expression for the bolometric energy of the flare, given as where σ B is the Stefan−Boltzmann constant.
Another way to calculate flare energy is to integrate the energy based on the variations in the normalized light curve, assuming the variations in the observed band are the same as that in any other band (e.g., Wu et al. 2015). Thus, the bolometric energy of the flare is given as Note that the calculated energy is usually larger than the blackbody assumption if we choose the flare temperature as 9000 K.
The calculated ratio of flare energies between the second method and the first is as follows: Taking a solar-like star (T eff =5800 K) and T flare =9000 K as examples, we estimate r = 2.93 in the TESS bandpass (600-1000 nm).