Statistical Analysis of Stellar Flares from the First Three Years of TESS Observations

In this paper, we study stellar light curves from the TESS satellite (Transiting Exoplanet Survey Satellite) for the presence of stellar flares. The main aim is to detect stellar flares using two-minutes cadence data and to perform statistical analysis. To find and analyze stellar flares we prepared automatic software WARPFINDER. We implemented three methods described in this paper: trend, difference, and profile fitting. Automated search for flares was accompanied by visual inspection. Using our software we analyzed two-minute cadence light curves of 330,000 stars located in the first 39 sectors of TESS observations. As a result, we detected over 25,000 stars showing flare activity with the total number of more than 140,000 flares. This means that about 7.7% of all the analyzed objects are flaring stars. The estimated flare energies range between $10^{31}$ and $10^{36}$ erg. We prepared a preliminary preview of the statistical distribution of parameters such as a flare duration, amplitudes and energy, and compared it with previous results. The relationship between stellar activity and its spectral type, temperature and mass was also statistically analyzed. Based on the scaling laws, we estimated the average values of the magnetic field strength and length of the flare loops. In our work, we used both single (about 60%), and double (about 40%) flare profiles to fit the observational data. The components of the double profile are supposed to be related to the direct heating of the photosphere by non-thermal electrons and back warming processes.


Introduction
Many stars show flares similar to solar flares, especially the main-sequence stars of spectral types G, K, and M. This is most likely due to a similar flare mechanism-magnetic reconnection. These stars also have a convective zone, which for the coolest objects can even span the entire star (Reid & Hawley 2005). The strength of the star's magnetic field is related to the motion of plasma due to convection (Benz & Güdel 2010). Stellar flares are highly energetic, rapid events that occur during magnetic reconnection in the stellar coronae. Stored magnetic energy is impulsively released, converted into other forms of energy and consequently radiated across the entire electromagnetic spectrum from radio to gamma-rays (Kowalski et al. 2019;Jain et al. 2011). During the impulsive phase of solar flares and their analogs occurring on late spectral-type stars, beams of nonthermal electrons, accelerated somewhere in the solar or stellar coronae, stream along magnetic lines toward the chromosphere, where they heat plasma by colliding with the dense matter near the feet of the loops. The interactions of the nonthermal electrons with dense matter also cause a strong and temporally variable emission of hard X-rays (HXR) in the feet of the loops (via the so-called bremsstrahlung process ;Brown 1971;Emslie 1978). The chromospheric matter, heated up to coronal temperatures, expands, filling the magnetic ropes (such a process is called chromospheric evaporation) and emits soft X-rays (SXR; Antonucci et al. 1984;Fisher et al. 1985). At the same time, a large amount of energy is radiated away from the loop feet at visible and ultraviolet wavelengths, significantly influencing the energy balance of the whole flare. The significant emission in the visual domain can be observed locally in solar flares (with high spatial and time resolutions), or it can be observed as a temporary variation of the integrated emission of the whole star.
The optical emission of a solar flare consists of a continuous spectrum (continuum) and chromospheric lines. During the majority of solar flares, the intensity of the continuum does not increase noticeably, while the intensity of spectral lines, formed mainly in the solar chromosphere and the transition region (TR), may increase significantly (Kowalski et al. 2015). The continuum optical emission is formed mainly in the solar/stellar photosphere, where the temperature and density are almost undisturbed during flares, and therefore the intensity of the continuum is not affected by the flares. On the contrary, hydrogen, calcium, magnesium, and other chromospheric lines are formed in the middle or higher chromosphere (Mihalas 1978), where the plasma parameters are strongly modified by different flare heating mechanisms. Therefore, the emission of these lines is strongly enhanced during flares. The main plasma parameters, which may influence the intensity of the spectral lines, are the temperature and density of the plasma. During solar flares, the values of both parameters increase due to nonthermal plasma heating and evaporation of the chromospheric plasma. The amount of such increases depends on the energy budget of all heating and cooling mechanisms.
The above-described situation is valid for most average-size solar flares. However, in some big flares the flux and the energy of the accelerated nonthermal particles are as height that heating of the low-lying solar photosphere or the lower chromosphere is possible. Then, the continuum optical emission may arise from some flare areas. Such flares are called "white-light" flares and they were observed in the past (e.g., Machado & Rust 1974;Hiei 1982;Neidig & Wiborg 1984).
Assuming the photospheric origin of white-light flare emission, the most probable mechanism that may be responsible for this continuum is H − emission (Ding et al. 1994). It is possible that the photosphere is not heated directly by nonthermal particles but rather during the process of "radiative back-warming", where the nonthermal electrons deposit their energy in the chromosphere and then the hotter chromosphere radially heats the lower-lying photosphere (Fang & Ding 1995). More recent studies show that the flare continuum emission can be formed in the thin "layer" of the chromosphere in the process of hydrogen recombination (Potts et al. 2010).
Research and observations of white-light flares from nearby stars have been rapidly developing in recent years thanks to large optical photometric surveys. For example, Davenport et al. (2014), Walkowicz et al. (2011), andHawley et al. (2014) used data from the Kepler mission, for which the time resolution was 1 minute (SC-short cadence) or 30 minutes (LC-long cadence). Other examples are Sloan Digital Sky Survey (SDSS; York et al. 2000), EVRYSCOPE (Law et al. 2014), All-Sky Automated Survey for Supernovae (ASASSN; Bersier 2016), and Next Generation Transit Survey (NGTS; Wheatley et al. 2018).
The duration times of white-light flares usually vary from minutes to hours, with a fast rise and a slower exponential decay profiles. Typical flare energies range from 10 26 to 10 32 erg. So-called super-flares with energies between 10 32 and 10 36 erg are also observed (Schaefer et al. 2000;Shibayama et al. 2013). Flare activity is thought to be correlated with fast rotation and stars of late spectral type. About 40% of M-type dwarfs are flaring stars (Yang et al. 2017;Günther et al. 2020;Howard et al. 2019). Moreover, detecting flares on late-type stars is easier because of the high flare contrast caused by a low surface temperature. The TESS mission provides us with the opportunity to study a large number and a variety of stellar flares.
In this paper, we examined light curves from the Transiting Exoplanet Survey Satellite (TESS) for the presence of stellar flares. We analyzed stellar flares from the first 39 sectors of the TESS mission. Our study investigates the preliminary statistics of stellar flares. It also includes relationships between stellar flares and such stellar properties as mass, effective temperature, and spectral type. In Section 2 we describe the TESS observations. Our software and methodology used to find flare candidates are detailed in Section 3. The results are presented and compared with previous findings in Section 4, and a discussion is provided in Section 5.

Observations
TESS (Ricker et al. 2014) is a space-based telescope that launched in April 2018. Its orbit is highly elliptical, with a period of 13.7 days. The main goal of the mission is to search for exoplanets using the transit method. TESS's instruments consist of four 10.5 cm telescopes, each covering a field of view (FoV) of 24°× 24°. In a single pointing TESS cameras cover a FoV of 24°× 96°. A given sector is observed through two satellite orbits, that is, for slightly less than a month. The first part of the mission, completed in 2020 July, consisted of observations of 26 sectors in both hemispheres covering about 85% of the sky. During the ongoing extended part of the mission, observations of the next 29 sectors (27-55) are carried out. The sectors partly overlap, which results in the presence of a continuous viewing zone in the vicinity of both ecliptic poles. During the primary part of the mission, two cadences were realized, 2 minute cadence for a sample preselected objects and 30 minute cadence for the whole covered FoV. The latter cadence is offered as Full Frame Images (FFIs) that can be used to extract photometry for any object in the FoV. During the extended part of the mission, two short-time cadences, 2 minute and 20 s, were chosen for selected objects, while FFIs are offered with a 10 minute cadence. The FFIs have an angular resolution of 21 arcsec per pixel.
Almost 330,000 stars observed with TESS in 2 minute cadence are in our sample, therefore we can study a huge number of stellar flares with a high signal-to-noise ratio (SNR). We selected calibrated, short cadence Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) light curves downloaded from the Mikulski Archive for Space Telescopes (MAST). Moreover, we required the quality flag to be equal to 0 and we normalized the light curves by dividing the flux by the mean flux of the star. If available we took the mass, radius, spectral type, and other stellar parameters from the MAST and SIMBAD databases when possible. All stars with spectral types earlier than F1 were rejected from our analysis. The effective temperatures of the flaring stars in our sample are limited to stars cooler than about 8000 K with masses smaller than about 1.7 M e . This limitation is enforced by numerous rapid changes in the light curves of hot stars, which are similar in shape to stellar flares. These changes were often misclassified by our algorithm as flares. The stellar variability of RR Lyrae, δ Scuti, γ Doradus, and SX Phoenicis is particularly problematic. After analyzing stars from the first 39 sectors of TESS observations, we found more than 25,000 flaring stars and 140,000 individual events.

Software
The software WARPFINDER (Wroclaw AlgoRithm Prepared For detectINg anD analyzing stEllar flaRes) prepared by us is written in the IDL (Interactive Data Language) and Python programming languages, which are mainly used for stellar flare identification and analysis. The software has also been adapted to download information automatically about stars from MAST and SIMBAD (Set of Identifications, Measurements and Bibliography for Astronomical Data), as well as data from the TESS satellite. Our automatic procedure also rejects most false detections such as asteroids, transits, or stellar pulsations. The appearance of an asteroid is usually observed in a stellar light curve as a strong, symmetrical increase of flux. Our software checks the asymmetry of each detection and rejects overly symmetrical phenomena. In addition, this shape is not well described by the flare profile we assumed. Many false detections occur during transits when the flux increases rapidly. To eliminate them, we added to the software the conditions regarding the relationship between the value of the signal at the beginning, at the maximum, and at the end of the flare. The most problematic, however, were the numerous false detections associated with short-period stellar pulsations. For this reason, we further limited our analysis to only stars with spectral types later than F0. The remaining false detections were rejected during partial visual inspection. The tables and diagrams containing data about the detected stellar flares and the flaring stars are the final results of using this software.

Methods
We prepared an automated, unsupervised three-step method for finding flare events from the TESS data. To search for potential flaring stars, we analyzed the 307,281 available targets with TESS PDCSAP light curves from the first 39 sectors. In our research we need data without systematic trends, so the PDCSAP flux was more appropriate than SAP (Simple Aperture Photometery) flux. In this work, we used light curves with a time resolution of 2 minutes.

Trend Method
One of the methods for finding stellar flares is the trend method. The first step is based on consecutive detrending of light curves using smoothing with several window lengths. This method was inspired by automated procedures described by Davenport et al. (2014) and tested on Kepler data. We smoothed the observed flux with the running average (smooth) function with a window length of 175 points (350 min). For this calculated trend we determined the standard deviation, and then calculated the new trend, rejecting all points protruding more than 3σ above the trend. We repeated this step recursively five times. For further analysis, we selected all points 1σ above the trend. We defined all the series of these points where at least four points protrude above the 2.5σ level as flares. They do not have to be consecutive points. There can be individual points between them where the deviation is less than 2.5σ, but greater than 1σ. The probability that a single point occurring randomly above the 2.5σ threshold is 0.0062. The likelihood of randomly getting 4 out of 7 such consecutive points drops to about 5 × 10 −12 . This gives us a chance for one false detection of a flare in one sector of properly detrended observations around 1 × 10 −7 .
We repeat above operation for the next five smoothing windows with lengths of 31, 75, 121, 221, and 311 points, respectively. The next flares determined in this way are added to the previous ones if the times of their maxima do not lie within the limits of the flares determined earlier. We have applied a few different values of smoothing windows because many flaring stars have more than one period of variability. These windows can be chosen as free parameters in our code and the values used here were selected experimentally on the data from TESS Sectors 1, 2, and 3. We checked these values for another set of thousands of flares and they seem to work pretty well for data with a 2 minute cadence.
There are many false detections caused by inadequate (usually too long) trend windows. Before the visual inspection, we verify these bad signals using other search methods in the next two steps of our automatic procedure.

Difference Method
The second method of finding stellar flares is the difference method whose idea was taken from the paper by Shibayama et al. (2013). The method is based on checking the flux difference between two consecutive points. In our method, contrary to the authors of the cited paper, we do not check the percentage distribution of differences, but the standard deviation instead. We assumed the study of points whose sigma is greater than 3σ as a default value. The spread of points is analyzed only for the positive values. All points that exceed 3σ are marked as potential flare detections and the times are compared with the hits detected by the trend method. It should be noted that only the mutual detection of probable flares in a designated time is treated by the software as a true detection and passed on to the third level of verification.

Flare Profile Method
The prepared list of potential stellar flares, which were detected by both previous methods, is then verified by a detailed analysis of their light curve profiles. In the first step of the profile method, we try to determine properly the beginning and the end time of a stellar flare. From the trend method, we obtain an estimated range of the flare duration. We use it to find the pre-flare and post-flare background approximated by a linear function. A linear background is sufficient in most cases. The standard duration of the flares usually does not exceed several dozen minutes. The percentage of stars rotating fast enough to affect the background determination is negligible in such a large sample of detected flares. Then we iteratively determine a new time of the start and the end of the flare and repeat the whole procedure 10 times. This value was selected on the basis of tests performed by us. In each iteration, we try to fit the profile for each test start and end point and choose the best fit based on the χ 2 value. Further, if the best solution is correct, we detected it as a flare. The final linear trend is fitted to the observational points and outliers are treated as a stellar flare.
In the next step, we fit the assumed flare profiles to the observational data. The first profile (Equation (1)) of a flare is given by equation (Gryciuk et al. 2017): where A, B, C, and D are parameters, t is the time, and erf is the error function. This profile is a convolution of a Gaussian profile (heating function) and exponential decay (function describing processes of energy loss). The A parameter is related to the amplitude of the flare. Parameter B determines the maximum energy release. Parameter C informs about the timescale of the energy release. The inverse of the D parameter is related to the flare decay time. We noticed that some of the flares do not fit well with only one profile. For this reason, we also performed a fit created from two similar profiles. In one case, we used profiles characterized by a small shift in the time of the maximum of both components (Equation (2)). These components have the same value as the B parameter.
In the second case, we allowed both components to have large shifts between the time of their maxima, and thus different values of the B parameter (Equation (3)).
We checked the quality of each of the flare profile fits. We used the χ 2 statistic for this purpose. The fit with the lowest χ 2 value was selected for a further analysis. For good fits, this value usually does not exceed 3.0 of the reduced χ 2 (taking into account the number of data points and the number of degrees of freedom). To distinguish a stellar flare from data noise we also used the probability density function of the F-distribution. We calculate the reduced χ 2 for the observational data with the fitted flare profile and the estimated trend. Next, we computed the cumulative distribution function for an F-distribution with the defined degrees of freedom. We accepted the detection as a flare if the probability is less than the cutoff value. After carrying out the earlier tests, we assumed that this value cannot exceed 0.1. Additional methods to reject false detections are: the calculation of the flare profile, skewness, and the bisectors of the profile. For the purposes of this work, we define the skewness in general as a measure of flare profile asymmetry. We assumed that a stellar flare should have a shorter rise time than its decay time. For this reason, the skewness of the profile should be a positive value. The bisector method gives information about the symmetry of the flare profile. Bisectors are defined by lines at the level of 10%, 20%, etc. of the maximum signal value. If the bisector method indicates an overly symmetrical profile (a value of about 0 or less), then the detection must be verified visually by a software user. Moreover, detections with a duration of less than 12 minutes (six points) are rejected.
After the analysis, separate data sets (three folders) with detections are created. The first one contains mostly real flares. The second set is for uncertain events that require visual inspection. The third set contains the detections rejected by the software.

Flare Injection-recovery Tests
Flare injection-recovery tests were performed to validate the results found via all three methods described in this paper. As in Günther et al. (2020), we decided to use two separate tests. For each of them, we selected randomly 1000 light curves of spectral type F, G, or K stars, as well as early M dwarfs and late M dwarfs. For the first test, we artificially added 10 timeseparated flares to each light curve. In the second test, we also inserted 10 flares, but they were only separated by 10 minutes up to a maximum of 1 hour. Each of the inserted flares was described by a single profile described above with randomly selected parameters. The range of this parameters was determined on the basis of the stellar flares we have observed. After proper preparation of the light curves, we tested our three-step software. In the case of both tests, we considered a flares as detected if its maximum was determined within the previously assumed duration. Figure 1 shows the test results we obtained. From the top three panels it can be seen that well detected flares usually have amplitudes greater than 0.01. Moreover, for M-type dwarfs, the mean recovery rate is lower. This is due to the lower brightness and SNR of the light curves of these stars. These dependencies are well illustrated by the three lower panels. It can be noticed that the recovery rate significantly decreases with the SNRs, TESS magnitudes, and effective temperatures of the stars.

Flare Energy
The maximum luminosity and the total energy of flares are estimated in two different methods. The first one requires a flare amplitude, a flare duration, a stellar luminosity, and a radius. Kowalski et al. (2016) showed using hydrodynamic simulations that a temperature of about 10,000 K is needed in order to reproduce correctly the flare emission on M stars in the range of white light. Thus, we assumed blackbody radiation  and an effective temperature of the flare (T flare ) of about 10,000 K (Shibayama et al. 2013;Mochnacki & Zirin 1980;Hawley & Fisher 1992).
The flare amplitude (C flare ¢ ) is defined as follows: where F star ¢ and F flare ¢ are the observed flux of the star and the flare, respectively, and A flare is the flare area, which are respectively given as: The bolometric flare luminosity (L flare ) and the energy (E flare ) could be calculated from, respectively: For a star for which information about its radius is available, it is possible to estimate the luminosity and energy emitted during the flare. According to Shibayama et al. (2013), the total energy uncertainty in this method is about 60%.
To estimate the energy of detected flares in the TESS detector bandpass we also used the modified method proposed by Kovári et al. (2007), and described by Vida & Roettenbacher (2018) and Vida et al. (2019). This method is based on the following integrating normalized flare intensity with subtracted background during the flare event: where t 1 and t 2 are the times of the start and end of the flare, respectively, and the integral ε TESS is the relative flare energy. Next, we had to estimate flux of the star in a selected interval of wavelengths (λ 1 , λ 2 ), which are the wavelengths of the TESS detector bandpasses. To estimate flux of a star we multiply the spectrum of the star ( ) l  , taken from ATLAS9 4 for stars g log( ), T eff (in spectra we assumed metallicity of the Sun and v turb = 2 km s −1 ), with the TESS response function S TESS (λ) multiplied by the star's area with radius R. Due to the use of the theoretical spectrum of stars from ATLAS9, there is no factor 4 in the following equation: To estimate the flare energy in the TESS detector bandpass E flare , we multiplied the star's flux F star in a selected interval of wavelengths by the relative flare energy ε TESS : Flare energies estimated in this way are not bolometric. For this reason, they have lower values than for the first method, which is based on Shibayama et al. (2013).

Problematic Cases
Our program for the automatic detection of stellar flares WARPFINDER may not work properly in a few cases. We observed a particularly large number of false detections on stars with spectral types earlier than F0. These are stars that are often characterized by rapid variation in their light curves. An example is RR Lyrae stars, whose periods of variation can be several hours. Moreover, the shape of the brightness change is similar to the flare profile assumed by us. It is characterized by a rapid growth phase and slow decay. Similar variability, classified as flares, is observed for δ Scuti, γ Doradus, and SX Phoenicis stars.
Flare detection is also difficult for very cool M spectral type stars. Due to their low brightness, light curves have a low SNR. Therefore, faint flares are not detectable by our software. This   conclusion is also confirmed by injection-recovery tests performed by us. Moreover, our flare catalog may include some false detections that are only observed once per star in all available sectors. This is due to a random spike in the signal accidentally classified as a flare. However, we decided not to exclude such stars from the analysis, so as not to remove information about real flaring stars.

Results
In this paper, we presented the results of the analysis of 329,176 stars from the first 39 sectors of TESS observations. We found 25,229 flaring stars and 147,683 individual flares.
The catalog with all detected flares is available in machine-readable form. Approximately 7.7% of all observed stars show flaring activity. This value varies depending on the spectral type of the studied stars and is greater than 50% for stars of spectral type M. On average, there are four flares per each flaring star observed during about 25 days (one TESS sector). We observed up to 69 flares on one star in one sector (TIC266744225; Figure 2). Moreover, for TIC150359500 we found 863 events during 23 sectors of observation. Figure 3 shows a histogram of the flare occurrence frequency per day. In most cases, one flare is observed on one star every ten days. There are, however, a few very active stars, where the average flare activity equals more than two events a day (TIC266744225, TIC142206123, and TIC441420236). A similar result was obtained by Davenport (2016), but for a   much larger sample of flares. In Table 1 we compare the number of stars observed in each of the 39 sectors with the number of flaring stars and individual flares. The hypothesis about a uniform distribution of stellar flares in sectors was rejected based on the performed Pearson's χ 2 test. Pearson's χ 2 test against the hypothesis of uniform distribution returns a value of 840 for 35 degrees of freedom, which exceeds the acceptable value for the assumed probability of the hypothesis. Moreover, there is an increase in the percentage of flaring stars during the TESS Extended Mission (Sectors above 27). This may be related to the change in the way of selecting observation targets. Now, about 80% of targets per sector come from Guest Investigator (GI) programs. In total, we found 9480 flaring stars during the first year of TESS observations (S01-S13, ecliptic decl. from −90°to −54°). For comparison, during the third year of the mission, when the Southern Hemisphere was re-observed (S27-S39), we found 10,661 flaring stars, including 8913 new objects.

Parameters of the Flaring Stars
The analysis carried out in this work allowed us to determine the distribution of the basic parameters of the flaring stars, such as effective temperature, spectral type, mass, and magnitude (Figures 4 and 5). All analyzed stars are marked in gray, and flaring stars are marked in blue. The analyzed data sample of flaring stars includes only stars with spectral types later than F0. For this reason, most of these stars' effective temperatures do not exceed 8000 K and 1.7 solar masses. As can be seen in Figure 4, the observed distributions of the presented parameters for the flaring stars are generally different than those for the entire sample of stars.
There is some similarity in the temperature distributions of the flaring stars, where the main maximum is about 3000 K and the second one about 6000 K. They are characteristic for the entire studied population. The distributions of the flaring stars' spectral types and masses are dominated by stars with spectral type M and masses of one solar mass or less. The distributions of the TESS, B, V, and U magnitudes for all stars and the flaring stars seem to be similar ( Figure 5). However, based on the two-sample Kolmogorov-Smirnov test we reject the hypothesis that these samples are of the same distribution. We used samples of the TESS, B, V, and U magnitudes for all stars and only for the flaring stars. The distribution of U magnitude of the flaring stars is the least similar to the entire sample distribution due to the fact that only about 10% of all stars in the sample observed by TESS have U magnitudes.
The parameters of the stars and the results obtained by us were examined using Benford's law (Benford 1938). This was done to check whether the data and the results are uniformly and randomly distributed. The tests show that the catalog data are not similar to a Benford distribution. The results obtained by us, such as the duration of the flares' or estimated flares energy, are consistent with Benford's law.
All stars observed by TESS with given luminosity and effective temperature were marked as light gray dots on the Hertzsprung-Russell (H-R) diagram ( Figure 6). The flaring stars with known luminosities, effective temperatures, and log(g) observed in Sectors 01-39 are marked with colors. The flaring stars with estimated luminosity and no log(g) are marked in dark gray. The luminosity of each star was estimated from its stellar radius and effective temperature. On the right is a colorbar of stellar log(g). The log(g) values of the flaring stars range from 3.3 to 5.1. The solid black lines represent the evolutionary tracks of stars with an initial mass of 0.5, 1, and 2 solar masses. Most of the flaring stars that we found are in the main sequence, but we detected flares on young stars and giants as well. To estimate the age of each star we use catalog data   from SIMBAD. We consider an object as young if the star is classified as pre-main-sequence, a young stellar object, T Tauri, or Wolf-Rayet star. The energies of the flares on young stars range from 8 × 10 31 to 4 × 10 35 erg. We compared the estimations of the energy of several super-flares on young stars with those from . For example, we determined the flare energy of TIC44678751 to be 2 × 10 35 erg and  found it to be 2 × 10 37 erg. There are similar big differences for other events, but we are not able to state clearly from what they may result. The two orders of magnitude difference cannot be explained only by the assumed flare temperatures (10,000 K in our case and 9000 K in ). The temperatures of the flaring stars mostly do not exceed 8000 K, but we found a few stars with higher effective temperatures. The results from our work (marked in blue in Figures 7, 8, 9, 10) were compared with works such as Doyle et al.  Table 2, that have estimated effective temperatures. The software prepared by us found more flaring stars in Sectors 1 and 2 than by Günther et al. (2020), mainly due to the different methods used for searching for flares. In both samples there are 963 flaring stars with the same TIC number and 7837 flares were found by both pipelines. However, it should be noted that our software analyses only events that have at least six observation points. The effective temperature distributions of the flaring stars agree nicely. Using the χ 2 goodness-of-fit test between the two distributions we obtain a value of about 0.90. There is no reason to reject the hypothesis at the 0.05 significance level.  Compared to  Figure 8), the temperature distributions of flaring stars are similar up to 8000 K (the χ 2 goodness-of-fit test value is about 0.93). We also chose the same sample that was presented by the authors in their work. The effective temperatures of most flaring stars that we have found are limited to this value due to the assumptions described earlier. There is a large difference in both samples' sizes, therefore on the vertical axis is the relative number of stars. Figure 9 shows the spread in spectral types within the solartype star samples from our work and Doyle et al. (2020; observations from Sectors 1-13). Again, the relative number of stars can be found on the vertical axis. The difference in numbers of stars is due to the sample selected by Doyle et al. (2020), in which spot modulation was searched. The different distributions of the spectral types of flaring stars may result from the different methods of searching for flares. In Doyle et al. (2020) the sample of the studied stars was small enough for each light curve to be checked visually, while our detections are mainly from using the automated software WARPFINDER.
In Figure 10 we compare the distributions of TESS magnitudes within the flaring stars observed in Sectors 1-20 from our work and . As in Figure 8, the slight differences in the TESS magnitude distributions are probably caused by the different methods of detecting flares used and the limitation of the stellar age up to 800 Myr.

Parameters of the Stellar Flares
Light curves of the detected flares were fitted with one of three types of profiles presented in the description of the flare profile method. We used a single profile, a double profile with the same value of B, and a double profile with a different value of B. More than half of the events fitted well with only the single profile. In other cases two profiles were needed for a good fit of the single flare. In this analysis, we did not handle the problem of fitting two or more overlapping events. Overlapping events are not fit well by the profiles we are using (Equations 1-3), so the χ 2 values are very high. For this reason, they are automatically rejected from the analysis. Based on the well-fitted flares (a value of the reduced χ 2 less than 3.0), we reconstructed three types of average profiles as a function of t 1/2 (Figure 11 and Figure 12), where t 1/2 is the full-time width at half the maximum flux. All flares were previously scaled to a relative time and amplitude. We then fitted the parameters A, B, C, and D of the obtained average profiles ( Table 2). A similar method was used by Kowalski et al. (2013) and Davenport et al. (2014). The distributions of the values of the reduced χ 2 for the analyzed stellar flares are shown in Figure 13.
The profiles differ in the average growth time (parameters B and C) and the decay time (parameter D). Longer and stronger events are usually fitted with two profiles. One of these profiles is supposed to be related to the direct heating of the photosphere by nonthermal electrons. It is characterized by very short growth and decay times (parameters B, C, and D). The second one, longer (with a significantly smaller D parameter) may be the result of back-warming processes (Isobe et al. 2007). It can be noticed that the contribution of the second component is particularly important for the profile with one B parameter. In both double profiles, the second component dominates in the decay phase.
The flare profiles that we obtained were compared with the flare model from Davenport et al. (2014), which is the result of the analysis of flares observed with Kepler for the M4 star GJ 1243. Figure 11 (right) shows our three profiles marked with black, blue, and green, and the reconstructed model of a flare proposed by Davenport et al. (2014;red). This analytic model was forced to go through 0 at t 1/2 = − 1 and 1 at t 1/2 = 0. The large differences between the profiles at the decay phase (t 1/2 > 2) may be related to the different methos of determining the end times of the stellar flares.
We examined the distributions of the basic parameters of the stellar flares. Figure 14 shows the distributions of the flare amplitudes, durations, and times of growth and decay. The amplitudes of most of the detected events do not exceed a value of 0.1 of the normalized flux, however, there are flares with amplitudes equal to several times the brightness of the star. The maximum of the distributions of all flare durations is approximately 50 minutes. The flares fitted with two profiles are characterized by a longer average duration. In this work, we limited the minimum of a flare duration to six data points (12 minutes). The longest observed events last up to several hours. Most of the detections with a duration of more than 500 minutes were false, so we did not include them in the further analysis. The average values of the growth times are much shorter than the decay times and are usually below 20 minutes. Examples of flares of long durations, long times of growth, and long decays are listed in Table 3.
Another result of the performed calculations is the estimation of the energy of most stellar flares. Figure 15 shows the cumulative energy distribution estimated using both of the previously discussed methods based on Shibayama et al. (2013;named v1; dark blue in Figures 15, 17, 20) and Kovári et al. (2007;named v2;light blue in Figures 15, 17, 20). Stars that have energies estimated by both methods must have such parameters as radius, effective temperature, and log(g). The sample of flares that have an energy estimated only by the first method is larger by about 3000 flares, because log(g) is not necessary for this calculation. The index of the power-law approximation is about 1.68 for the method based on Shibayama et al. (2013) and about 1.65 for the method based on Kovári et al. (2007) for flare energies from 5 × 10 33 to 10 36 erg. These results are in agreement with those from previous papers (Lacy et al. 1976;Hawley et al. 2014;Maehara et al. 2020). The flare frequency distributions (FFDs) of stars binned by the flare amplitude are presented in Figure 16. We examined how the slopes of the distributions vary for three samples of stars. One of them contained stars with masses less than 0.3 M e (violet), second one with masses in between 0.3 M e and 0.5 M e (blue), and the other with masses greater than 0.5 M e (beige). The slope values of the distributions are very similar to those in other works, and is about −1.4 ). It can be seen that the sample of stars with lower masses also has a shallower slope. The differences in slopes can be explained by theoretical works ) Figure 17 shows a histogram of flare energies. It could be noticed that the results obtained with the first method (Shibayama et al. 2013) are larger than those obtained with the second method (Kovári et al. 2007). The estimated energies of the stellar flares range from 10 31 to 5 × 10 36 erg, which means that most of the detected events are super-flares. Figure 18 additionally shows the dependence between the flare energies estimated using the method based on Shibayama et al. (2013;v1) and based on Kovári et al. (2007;v2). Our results were compared with the works of Doyle et al. (2020), Günther et al. (2020), and Howard et al. (2019). These studies also investigated stellar flares observed by the TESS satellite. The comparison shows that our results agree with the already known average energy of the detected stellar flares. Figure 19 shows the differences between the energy distributions of flares fitted with one or two profiles. Based on the Kolmogorov-Smirnov test we reject the hypothesis that they came from the same distribution. The flares with the highest energies (greater than 10 35 erg) are usually fitted with two profiles. Figure 20 shows the flare energy distributions calculated in our work and that from Doyle et al. (2020). The energies range from 10 32 to 10 36 erg and the maximum of the distribution from Doyle et al. (2020) is for energies about four times higher than in our work. This is due to the fact that we have modified the method of estimating the energy of the flares described by Kovári et al. (2007). The flare energies obtained in our sample are very similar to those of Günther et al. (2020; Figure 21). Small differences could be attributed to the different methods of searching for stellar flares. Moreover, in Günther et al. (2020) flares of shorter duration (from six minutes) than in our work were included. Figure 22 shows another comparison. Howard et al. (2019) estimated flare energies based on TESS observations in Sectors 1-6 and from Evryscope. The average energy of the flares from both data sets is significantly different. The authors explain that TESS observed much fainter, lower amplitude events. Our estimations are similar to the energies obtained by other authors from observations made by TESS. Based on all the presented results, it can be seen that the energies of the stellar flares vary from 10 31 erg to 10 36 erg, depending on the analyzed star sample. Table 4 additionally summarizes information about the stars with the most energetic flares found by our software. Among the super-flares listed here, we found only one that was already mentioned in other works. We estimated the energy of the flare from TIC79358659 at 10 35.88 erg, while in Doyle et al. (2020) it was 10 35.72 erg. We would like to note that the systematic overestimation of the flare energy is caused by too large  Notes. a Probably no overlapping events b Appear to be overlapping events estimates of the stellar radius from the MAST catalog. In the analyzed sample of stars, there were objects classified as mainsequence stars, but with radii of several or a dozen solar radii. For this reason, we limited our analysis only to stars with a radius smaller than 2 R e . Examples of super-flares on TIC175305230 and TIC232036408 with fitted profiles are shown in Figures 23 and 24. The flares were fitted with two profiles with one B. Assuming that a bolometric energy of 10 33 erg corresponds to an SXR scale X100 flare (Benz 2016), the detected events are more energetic than X10,0000 flares. Using the energy estimation method proposed by Shibayama et al. (2013), it is also possible to determine the areas of the stellar flares (Equation (7)). Figure 25 shows a distribution of the areas of the flares obtained in our research in parts per million [ppm] of the stellar disk. All flares are marked with gray, flares fitted with one profile with black, two profiles with one B with blue, and two profiles with B1 and B2 with red. The maximum of each of the distributions is about 2200 ppm. Based on the Kolmogorov-Smirnov test we rejected the hypothesis that they arise from the same distribution. For the fit with one profile, the average value of the flares areas is 0.31% of the entire stellar disk, and for the fits with one and two B parameters it is 0.35% and 0.27% respectively. Flares fitted with two profiles with one B parameter are characterized by, on average, larger areas than the other flares. Our results are in line with the general knowledge of the areas of solar and stellar flares observed in white light (Heinzel & Shibata 2018).
An analysis similar to those of Namekata et al. (2017) and Tu et al. (2021) was also carried out in our work. These authors compare the durations and energies of the flares observed on Figure 15. Cumulative distribution of the flare energies from TESS Sectors S01-S39. The energies, the fitted power function, and the power-law index estimated using the method based on Shibayama et al. (2013;v1) are marked in black. The energies, the fitted linear function, and the power-law index estimated based on Kovári et al. (2007;v2) are marked in blue.   Shibayama et al. (2013;v1) are marked in dark blue and those based on Kovári et al. (2007;v2) are marked in light blue. Figure 18. The dependence between the flare energies estimated using the method based on Shibayama et al. (2013;v1) and based on Kovári et al. (2007;v2). The red line represents the same energies in both methods.
the Sun and other solar-like stars. Based on the scaling laws, they estimated the magnetic field strengths and the lengths of the flare loops. The authors assume a constant pre-flare coronal density. Other versions of the scaling laws, in which the coronal density changes with the length of the flare loop, have also been proposed. Moreover, the coronal density of very active stars may be much greater than the solar coronal density. This can be concluded from their short rotation periods and the emission in the X-ray range. Therefore, the obtained results are only estimates of the strength of magnetic fields and the lengths of the flare loops. Scaling laws used in this work are given by the following equations (Namekata et al. 2017): where τ is the flare duration time, E is the estimated flare energy, B is the magnetic field strength, and L is the length of the flare loop. Figure 26 shows the observed relations between flare energy and duration from our sample of stars. On this graph, we have plotted theoretical lines for the different strength of magnetic fields (30, 60, 200, and 400 G) and the lengths of the flare loops (10 10 cm and 10 11 cm). The upper panel shows flares fitted with one profile, the middle one shows flares fitted with two profiles with one B, and the lower panel shows flares fitted with two profiles with B1 and B2. The location of the center of mass with its coordinates is marked with red. The centers of mass shifted toward greater flare energies and longer durations. Moreover, the crosses marked in blue mean the average energy of the flares in a given range of their duration (10-30 minutes, 30-50 minutes, etc.). The linear fit to these points with their parameters is also marked in blue. The slope of the line changes from 0.49 to 0.46 for the flares fitted with different profiles.  For flares fitted with two profiles, the dependence between the flare energy and the duration is stronger. Flares for which it is necessary to distinguish two components of the profile are usually stronger, longer, and more complicated. We assume that there are two mechanisms in such flares: direct heating of the photosphere with nonthermal electrons and back-warming processes. We assume that there are two mechanisms in such flares: direct heating of the photosphere with nonthermal electrons and back-warming processes of thermal origin. For example, observations of flares for the AU Mic star with the FUSE (Far Ultraviolet Spectroscopic Explorer) satellite shows an increasing continuum flux to shorter wavelengths, which can be interpreted as free-free emission from hot plasma (Redfield et al. 2002). Our results are similar to those of Namekata et al. (2017). The estimated energies and the durations of the flares agree very well with the analogous values for the super-flares observed by Kepler with a 30 minute cadence. The high values of the magnetic field strengths up to 400 G determined by the scaling laws are characteristic for shorter super-flares observed by Kepler with a cadence of 1 minute.
We assume a constant pre-flare coronal density. Using information about the flares' durations and its energies we calculated the magnetic field strengths and flare loop lengths. The average values of the magnetic field strengths are from 10 G to 200 G. The estimated lengths of the flare loops are from 10 10 cm to 2 × 10 11 cm. These values correspond to several radii of a given star (Figure 27). It could be concluded that the single-loop assumption is not sufficient for most of the observed flares. Stellar flares are probably associated with several or more flare loops. This conclusion should be correct, taking into account the observations of strong events occurring on the Sun (Aschwanden et al. 2007).

Discussion
We analyzed the first 39 sectors of TESS satellite data (about 330,000 stars). More than 140,000 flares were detected on over 25,000 stars. The flaring stars that we found make up about 7.7% of the entire data sample. This value varies greatly depending on the effective temperature or the spectral type of a star. Most of the found flaring stars are in the main sequence,   but we detected flares on young stars and giants as well. The effective temperatures for most of the flaring stars do not exceed 8000 K. The percentage of flaring stars could be even greater than 50% for the stars of spectral type M. In total, during the first year of TESS observations (S01-S13, ecliptic decl. from −90°to −54°), we found 9480 flaring stars. For comparison, during the third year of the mission, when the Southern Hemisphere was re-observed (S27-S39), we found 10,661 flaring stars, including 8913 new objects.
In our work, we used three types of flare profiles: single, double with one B, and double with B1 and B2. A single profile was sufficient for more than half of the events. The maximum of the duration distribution of all flares is approximately 50 minutes. The flares fitted with two profiles are characterized by a longer average duration. The average value of the durations of flares fitted with one profile is 66 minutes, and with two profiles is 98 minutes. The longest observed events last up to several hours. The average values of the growth times are much shorter than the decay times and are usually below 10 minutes.
Usually longer and stronger events are fitted with two profiles. One of these profiles is supposed to be related to the direct heating of the photosphere by nonthermal electrons. It is characterized by very short growth and decay times (parameters B, C and D). The second one, which is longer (with a significantly smaller D parameter) may be the result of backwarming processes (Isobe et al. 2007). In both double profiles, the second component dominates, especially in the decay phase.
We also estimated the energies of stellar flares based on two different methods (Shibayama et al. 2013 andmodified Kovári et al. 2007). In the first method, we assumed a flare temperature of 10,000 K, while in the second method no such assumption was needed. In most cases, the energy calculated using the stellar spectrum gives lower values. This is due to the fact that the energy estimated in this way is not a bolometric energy. The estimated energies of the stellar flares range from 10 31 to 10 36 erg, which means that most of the detected events are super-flares. Our results were compared with the works of  Figure 21). Small differences could be caused by the different methods of searching for stellar flares. Moreover, in Günther et al. (2020) flares of shorter duration (from six minutes) than in our work were included. Our estimations are similar to the energies obtained from TESS observations by other authors. The index of the power-law approximation of the cumulative distribution of the flare energies ( Figure 15) is about 1.7 for the method based on Shibayama et al. (2013) and about 1.5 for the method based on Kovári et al. (2007) for flare energies from 5 × 10 33 to 10 36 erg. These results are in agreement with previous papers (Lacy et al. 1976;Hawley et al. 2014;Maehara et al. 2020).
Using the energy estimation method of Shibayama et al. (2013), it is also possible to determine the areas of the stellar flares. The maximum of this distribution for all flares is about 2200 ppm (0.22%; Figure 25). Flares fitted with two profiles with one B parameter are characterized by, on average, larger areas than the other flares. For the fit with one profile, the average value of the flares' areas is 0.31% of the entire stellar disk and for fits with one and two B parameters it is 0.35% and 0.27% respectively. Our results are in line with the general knowledge of the areas of solar and stellar flares observed in white light (Heinzel & Shibata 2018).
The index of the power-law approximation of the relation between the flare energy and duration changes from 0.49 to 0.46 for the flares fitted with different profiles (Figure 26). For the flares fitted with two profiles, the dependence between the flare energy and the duration is stronger. Flares for which it is necessary to distinguish two components of the profile are usually stronger, longer, and more complicated. Our results are similar to those of Namekata et al. (2017). The estimated energies and the durations of the flares agree very well with the analogous values for the super-flares observed by Kepler with a 30 minute cadence. The high values of the magnetic field strengths up to 400 G determined by the scaling laws are characteristic for shorter super-flares observed by Kepler with a cadence of 1 minute. On the basis of the relationship between the length of the flare loops and the radius of the star (Figure 27), it could be concluded that the single-loop assumption is not sufficient for most of the observed flares. It looks as if stellar flares are probably associated with several or more flare loops. This conclusion should be correct, taking into account the observations of events occurring on the Sun.
Flares detected on cool stars have, on average, higher amplitudes due to the greater contrast between the flares' and the star's surfaces. Such an effect has already been noted in other studies that analyzed stellar flares (Jackman et al. 2021). Figure 28 shows the relation between the logarithm of the amplitude of the flares and the logarithm of the effective temperature of the stars. The linear fit to this data and the linear model parameters are shown in blue. The linear model parameters are shown in the upper-right corner. The large amplitude of a flare is not always related to the high energy emitted during its duration. For this reason, it is important to use methods to estimate the energy of stellar flares. Our results are consistent with previous studies on the topic based on the observations of TESS and Kepler. Figure 26. Observed relations between flare energy and duration with theoretical (dotted) lines for different strengths of the magnetic field and lengths of the flare loops. A constant pre-flare coronal density is assumed. In panel (a) the marked flares are fitted with one profile, in panel (b) the flares are fitted with two profiles with one B, and in (c) they are fitted with two profiles with B1 and B2. The discrete flare duration bins result from the time resolution of the observations, equal to 2 minutes. In each case, the red color shows the location of the center of mass with its coordinates. The blue color shows a linear fit to the data. More information is in the text.