Superflares, chromospheric activities and photometric variabilities of solar-type stars from the second-year observation of TESS and spectra of LAMOST

In this work, we still use effective temperature ($T_{\mathrm{eff}}$) and surface gravity ($\log g$) as criteria to select solar-type stars. Following the criteria used by other works (e.g. Shibayama et al., 2013; Maehara et al., 2015; Tu et al., 2020), we also set the effective temperature rang from 5100 to 6000 K and surface gravity $\log g>4.0$. Note that Sun-like stars are also solar-type stars but with stricter criteria as $5600\mathrm{K}\leqslant T_{\mathrm{eff}}<6000\mathrm{K}$ and stellar periods $P>10$ days. Stellar effective temperature and surface gravity are obtained by using the final version of the TESS input catalogue (TIC v8; Stassun et al. (2019)). Then, 94 stars are excluded from our original data set as these stars are flagged as binary systems, after cross-matching with the Hipparocos-2 catalog (van Leeuwen, 2007). Next, we cross-match the data set with the Gaia Early Data Release 3 (Gaia-EDR3; Gaia Collaboration et al. (2020)). Note that stellar effective temperatures of sources in Gaia-EDR3 are all collected from Gaia-DR2 (Gaia Collaboration et al., 2018). Through cross-matching with Gaia-EDR3, we exclude those targets that may be contaminated by other unrelated brighter stars within 21^′′ distance. Those unrelated stars are nontargeting objects in TESS pixel resolution. Because TESS has a pixel resolution of 21^′′, which is much larger than that of Kepler, we should exclude those brighter stars, which may not be distinguished by the TESS cameras and contaminate our targets inside of 1 pixel. For more details, one can refer to Section 2.1 and Appendix A of Tu et al. (2020). We will also briefly discuss the contamination of M dwarfs in section 3.1. After these steps, an additional 209 stars are removed from our database. In total, we select 22,539 solar-type stars for further searching of superflare candidates.

Every single presearch data-conditioned (PDC) light curve of these solar-type stars from TESS has been checked by using the automatic searching algorithm. We have successfully used this method to process TESS first-year observations with 2 minute cadence. This method is obtained by detecting the brightness variations of a star through changes between pairs of consecutive points, which is expressed as

where $i$ stands for the $i$th data point. If $\left(F_{i}-F_{i-n-1}\right)>0$ and $\left(F_{i+1}-F_{i-n}\right)>0$, then $s=1$; otherwise, $s=-1$. Here $n$ describes the number of data points between two pairs of consecutive points. We set $n=2$ and $4$ to search each superflare candidate light curve for two times (Tu et al., 2020). Those points, of which $\Delta F^{2}$ is three times larger than the upper $1\%$ of the $\Delta F^{2}$ distribution, are selected as superflare candidates. Then, we use the quadratic function $F_{q}(t)$ to fit the light curves of these candidates for further calculation of superflare properties. One may refer to Section 2.2 of Tu et al. (2020) for more specific illustrations.

Comparing with the visual inspection procedure we used for the TESS first-year observation, in this work, we update some details for visually inspecting TESS pixel-level information. Specifically, in Figure 1, compared with our previous work, we add three additional panels to improve the accuracy of the visual inspection. The first one is shown in the upper right panel of Figure 1, which marks the specific moment of the detected superflare candidates in the whole light curve. In the lower panels, we not only show the pixel-level image at the flare peak moment but also attach the image at the preflare time. Therefore, we may easily be sure whether photometric apertures from the TESS pipeline are stable or not. For the pixel-level flux in Figure 1, we also attach the normalized flux for each corresponding pixel. This improvement reduces noises from other unrelated pixels and bolsters superflare signals. From this panel, it is obvious that the results can be described by the point-spread function (PSF). And it is more efficient to confirm those superflare candidates, which show the standard shape of rapidly rising and slowly decaying, as shown in the the upper left panel of Figure 1. Through visual inspections, we exclude many false-positive candidates and retain 1272 superflares for further research.

We estimate stellar luminosity through the Stefan–Boltzmann law

where the surface effective temperature and stellar radius are denoted as $T_{*}$ and $R_{*}$, and $\sigma_{\mathrm{sb}}$ is the Stefan–Boltzmann constant. The superflare energy is calculated by

where $F_{\text{flare}}$ is calculated by

Here $F_{\rm q}{(t)}$ is the flux fitted by quadratic function. The start and end times of each superflare are the first and last points, of which three times the photometric errors is less than $F_{\rm flare}(t)$,

where $F_{\rm error}(t)$ is the photometric error, which is provided by the TESS pipeline. Superflare duration is defined as the time interval between the start time and end time, which is also the integrating range of Equation (3). As in our previous work, we also use the Lomb-Scargle method to estimate stellar periods (Lomb, 1976; Scargle, 1982). The false-alarm probability is set to be $10^{-4}$ to search the most probable period for a star (e.g. Cui et al., 2019).

Many works have also used different methods to measure the amplitude of the up-down oscillation in the light curve. For example, He et al. (2015, 2018) used the rms of the stellar fluxes to estimate the effective fluctuation range of a light curve. The amplitude range of light curve can be understood as surface magnetic features rotating with the star (e.g. I\textcommabelowsık et al., 2020; Shapiro et al., 2020). In this work, we denote the photometric variability range as $R_{\rm var}$. We use the 95% and 5% ranked normalized flux as up and down limits at first (Basri et al., 2010). Then, the $R_{\rm var}$ is deduced by calculating the difference of the two limits as

where $F$ stands for the ranked normalized flux of the whole light curves of the star. Note that, in some literature, $R_{\rm var}$ values are multiplied by 100 and in units of percentile, which we do not practice in this work. For all 22,539 stars, $R_{\rm var}$ is calculated through the corresponding light curves.

LAMOST has been surveying the sky since 2012, mainly in the northern hemisphere of the sky (Cui et al., 2012; Deng et al., 2012; Zhao et al., 2012). It released over 10 millions low-resolution spectra (LRS) with wavelengths from 3700 to 9000 Å. In 2018, LAMOST started a new 5 yr medium-resolution spectroscopic (MRS survey. And in 2019, the LAMOST DR6 was released and provided MRS spectra that cover the wavelength range from 4950 to 5350 Å in the blue arm and 6300 to 6800 Å in the red arm (Liu et al., 2020). We use LRS data from 3890-4020 Å of LAMOST DR7 v1.1 and DR8 v0 in this work. The DR7 v1.1 covers all LRS data from 2012 to 2018, and DR8 v0 only covers data observed since 2019.¹¹1http://dr.lamost.org/data-release.html

We cross-match all 22,539 solar-type stars in our database with these two LAMOST data releases. We only keep those spectra according to three criteria. (1) The signal-to-noise ratio (S/N) of the spectral blue part should be larger than 10 to ensure accuracy of further analysis. (2) The distance between the targets of TESS and corresponding spectrum of LAMOST should be less than 1^′′. (3) As there might be more than one spectrum matching one TESS target, we only adopt the spectrum that is the nearest one around the target and the latest one observed. In total, 7454 solar-type stars have matched LAMOST spectra, of which 79 stars have superflares. Spectral data are now available for further estimation of chromospheric activity.

Eberhard & Schwarzschild (1913) found a sharp reversal at the K line of Ca from those stellar spectra, which were similar to the Sun. Then, Wilson (1968) used photometric results, which showed reliability for using Ca II H and K lines for indicating stellar chromospheric activity. To measure these two lines, Vaughan et al. (1978) defined a method to calculate the mean flux of the H and K lines. This method laid foundation for the best-known $S$-index, a reliable indicator of stellar chromospheric activity. In this work, we also use the LRS from LAMOST observations. Following Karoff et al. (2016), we use the reformed version of the $S$-index, which is expressed as

where 8 is the ratio of exposure time between channels of H and K lines and other referenced channels of the Mount Wilson spectral photometer (Karoff et al., 2016). Specifically, $\tilde{H}$ and $\tilde{K}$ are the mean fluxes of the triangle bandpass with $1.09\dot{\mathrm{A}}$ FWHM, so $\Delta\lambda_{\mathrm{HK}}=1.09{\rm Å}$. These band-passes center at 3968 and 3934 Å. $\tilde{R}$ and $\tilde{V}$ are mean fluxes of bandpasses centered at 3901 and 4001 Å, respectively. These two bandpasses are with 20 Å width, so $\Delta\lambda_{{RV}}=20{\rm Å}$. Karoff et al. (2016) proved that by taking $\alpha=1.8$ for the spectrum from LAMOST, their $S$-indexes have an identical distribution to that of the Isaacson and Fischer sample (Isaacson & Fischer, 2010). So, we also take the normalization constant $\alpha=1.8$ from Hall et al. (2007). We use $\log\sigma(S)=-\log{\rm S/N}-0.5$ to roughly estimate the uncertainty of the $S$-index, where ${\rm S/N}$ is the signal-to-noise ratio of the spectral blue part. This relation is obtained through $\chi^{2}$ minimization in the fitting relation between standard deviation of the $S$-index and ${\rm S/N}$ (Karoff et al., 2016).

We gather all superflare energy distributions in Figure 2, including M dwarf flares observed by TESS (Günther et al., 2020) and Kepler (Yang et al., 2017). For superflares on solar-type stars, we import superflares from Kepler (Shibayama et al., 2013) and TESS first-year observations (Tu et al., 2020). The normal distribution has been fitted for $E_{\rm flare}$ in log scale. The probability density function of the normal distribution is expressed as

where $\mu$ is the expectation of this distribution, and $\sigma$ is the standard deviation.

We have all solar-type stars in this work cross-matched with Gaia-EDR3 to select those TESS targets that may contain M dwarf candidates. However, they are not excluded from our database for the following reasons. (1) The M dwarfs are stars with 3000 K $<T_{eff}<$ 4000 K, and $\log g>$4.0 (e.g. Yang et al., 2017). However, the stellar parameter of the surface gravity ($\log g$) is not included in the Gaia-EDR3 catalog. So, those selected TESS targets only contain M dwarf candidates. We keep them in our database and mark them with flags, as shown in Tables 1 and 2. (2) As Kepler has better pixel resolution than TESS, in Figure 2(b), we compare the energy distributions of superflares on solar-type stars (Shibayama et al., 2013) and M dwarfs (Yang et al., 2017) from Kepler observations. Apparently, superflares on solar-type stars (gray histograms) basically have energies over 2 orders of magnitude larger than those of M dwarfs (yellow histograms). As shown in panel (a) of Figure 2, the superflares of this work still apparently show higher energies than the superflares on M dwarfs observed by TESS (Günther et al., 2020). This result indicates the low possibility for flares on M dwarfs contaminating our superflare database, even if they generate flares more actively than solar-type stars. (3) Here is the most important reason. We use a database of superflares on M dwarfs from the first 2 months of the TESS mission (Günther et al., 2020). The flare peak luminosity of 5398 superflares on M dwarfs is calculated using

where $F_{\rm{flare}}(t_{\rm peak})$ is the amplitude at the peak moment of the superflare on M dwarfs. The largest $L_{\rm{peak}}$ is equal to $7.22\times 10^{32}\,{\rm erg\,s^{-1}}$ from TIC 260506296, which was also specifically shown in the Figure 10 of Günther et al. (2020). In our database, the lowest quiescent stellar luminosity ($L_{*}$) of flare stars is equals to $9.67\times 10^{32}\,{\rm erg\,s^{-1}}$, which is higher than the largest flare peak luminosity ($L_{\rm{peak}}$) of M dwarfs. And only four nonflaring stars in our database have smaller stellar luminosities than the largest $L_{\rm{peak}}$ of the M dwarfs. In short, the superflares on M dwarfs may not affect the observation of solar-type stars, even if they cannot be distinguished in just 1 pixel. Therefore, the possibility that our database is contaminated by superflares on M dwarfs is very low. (4) As TESS has 21^′′ pixel resolution, we may not be able to distinguish M dwarf candidates apart from the main target. Statistically, only 7.23% (92 of 1272) superflares are from 9.97% (31 of 311) flare stars the may contain M dwarfs within 21^′′ radius. Specifically, all M dwarf candidates of these 31 flare stars are 1.5 mag fainter (0.25 flux of the Gaia band) than the corresponding flare star. Besides, through cross-matching with LAMOST data, none of these 31 flare stars have M dwarf candidates within 21^′′ radius. As these superflares do not occupy the main part of the whole database, they could not significantly influence the statistical results in the following sections.

As the bandpass filter of TESS (600-1000 nm) is less sensitive to shorter wavelengths than Kepler filter (420-900 nm), TESS is less capable of detecting relatively small flares than Kepler do (Doyle et al., 2019; Tu et al., 2020). We also get a similar conclusion from Figure 2(b). In this work, as shown in Figures 2(c) and (d), the data of the TESS second-year observation (red histograms) get $\mu\sim 34.67$ which is 0.08 higher than the $\mu\sim 34.59$ of the Kepler data (gray histograms), but 0.31 less than the $\mu\sim 34.98$ of the TESS first-year observation (blue histograms). However, TESS observes superflares with slightly higher energies than those of Kepler data. The difference of $\mu$ between superflares from TESS and Kepler is reduced from 0.39 (TESS first-year observation) to 0.08 (TESS second-year observation), which may be a result of the improved accuracy of visual inspection in this work, as described in Section 2.1 and Figure 1. The improvement bolsters flare signals by reducing noises from other unrelated pixels and makes those relatively less energetic superflares more easily to be seen. Also, TESS has upgraded their data processing pipeline.²²2https://heasarc.gsfc.nasa.gov/docs/tess/reprocessed-tess-data-coming-soon.html In this work, we use the reprocessed data for sectors 14-19 of TESS. Besides, in Section 3.6, the distribution of flare peak amplitudes shown in Figure 7(a) also indicates that the superflares in our database have slightly lower peak amplitudes than those of previous works, which pulls down the energy distribution of superflares in this work. Furthermore, in Section 3.8, the distribution of $R_{\rm var}$ of this work (red histograms) shown in Figure 12(c) has an overall lower $R_{\rm var}$ than that of our previous works (blue histograms). From the scatter plot of Figure 12(d), we may find a positive relation between maximum flare energy and $R_{\rm var}$. The overall lower $R_{\rm var}$ may also be the reason why the superflare energy in this work is slightly lower than that of Tu et al. (2020).

Table 4 lists counts of solar-type stars, superflares, and flare stars in each period bin. More intuitively, period distributions are shown in Figure 3(a). It shows that the period distributions are not changed significantly from the TESS first-year observation, as shown in Figure 3(b). For each period bin of Figures 3(a) and (b), the corresponding normalized flare star fraction (${N_{\rm fstar}}/{N_{\rm star}}$) is shown in panels (c) and (d). It is evident that the coverage of star spots is negatively related with stellar periodicity (e.g. Figure 15 of Notsu et al. (2019), and Figure 6 of Okamoto et al. (2021)). In other words, stellar magnetic activity is negatively related to stellar period. From panels (c) and (d) of Figure 3, we also see the decreasing trend of the fraction of flare stars with period $P>10\;{\rm days}$. So it might be the evidence that magnetic features like star spots are key to generate superflares (Shibata et al., 2013). For those stars with period $P>10\;{\rm days}$, the possibility of producing large magnetic features (e.g. star spots) is low, which may cause lower detectability of their stellar periods. So, the number of stars with period $P>10\;{\rm days}$ decreases dramatically, as shown in panels (a) and (b) of Figure 3.

As introduced in Section 2.2, we use the Lomb-Scargle method to estimate stellar periodicity, which is mainly dependent on the brightness variation of stellar light curves. Note that this method may be biased toward selecting stars with large star spots that cause a brightness variation of the stellar light curves. For stars without detectable variability, the periodicity may not be accurately estimated compared with the spectral method. So, it is possible that the real fraction of rapidly and slowly rotating stars in the observed fields may be different from that deduced from the Lomb-Scargle period estimation (Notsu et al., 2019). We refer to the result from Table 8 of Tu et al. (2020), which estimated the number of fractions through the empirical gyrochronology relation. Honestly, this relation is empirical and not totally reliable (e.g. Tu et al., 2015; van Saders et al., 2016). Appendix B of Tu et al. (2020) discussed this relation in detail; one could refer to it for more information. In order to compare, we set four data sets according to stellar effective temperatures and periods. The related results from Notsu et al. (2019), Tu et al. (2020), and this work are all collected in Table 5. From this table, we find that the fraction of $N_{\rm star}(P>10\;{\rm days})/N_{\rm all}$ in this work is basically the same as that from the TESS first-year observation (Tu et al., 2020). After comparing our results with those deduced by the gyrochronology relation and Notsu et al. (2019), we still have a lower fraction of $N_{\rm star}(P>10\;{\rm days})/N_{\rm all}$ than theirs. However, this will not significantly affect the following results related to stellar periodicity because (1) the empirical estimation may be not reliable (e.g. Tu et al., 2015; van Saders et al., 2016) and (2) the differences between the results of gyrochronology relation and ours are less than 1 order of magnitude. However, the observation mode of TESS is unique, which may limit the stellar period estimation due to some stars are observed for only 27 days. A similar discussion can also be found in the Appendix B of our previous work (Tu et al., 2020). Therefore, we expect that accurate periodic estimation that can be done for TESS targets, like that for Kepler targets (McQuillan et al., 2014).

In this work, we also calculate the flare frequency for a single flare star as

where $N_{\text{*flares}}$ is the count of superflares on a star, and $\tau_{{*}}$ is the effective observing time length for the star. We sort 311 flare stars according to their flare frequency ($f_{*}$) and number of superflares ($N_{\text{*flares}}$) in Table 2. For two stars, TIC 431107890 has the highest flare frequency, and TIC 394030788 exhibits the most superflares. The light curves of their most active parts are shown in Figure 4. Here we use two active stars (TIC 43472154 and TIC 364588501) from Tu et al. (2020) for comparison. Star TIC 431107890 has the highest flare frequency in this work, but it is still less active than TIC 43472154 which generates 16 superflares within just one observing sector. Star TIC 394030788 produces 76 superflares in just eight sectors, and this star is more active than the star TIC 364588501, which generates 63 superflares in 13 sectors. Statistically, we do not find any stars in this work that generate superflares as much as many superflares as TIC 364588501 produced (20 superflares) in the sector 5. However, TIC 394030788 from this work is still worth analyzing in detail, like Doyle et al. (2020) did for TIC 364588501.

In this work, the most energetic superflare is the one that releases about $1.24\times 10^{36}$ erg energy from TIC 420137030. However, this is 1 order of magnitude less than the largest superflare with energy $1.77\times 10^{37}$ erg from TIC 93277807 (Tu et al., 2020). So, the star TIC 93277807 is still interesting for exploring its possible superflare mechanisms. Also, we will discuss this topic in Section 3.8. We still expect that ground-based follow-up observations and TESS extended missions could provide much more valuable data about this star.

Stellar activities are significantly relevant to the habitability of planets, and star-planets interactions will improve our understanding of space weather (e.g. Airapetian et al., 2016; Riley et al., 2018; Linsky, 2019), and planet atmosphere changes by impact of stellar flares (Chen et al., 2020).

It has been argued that close-in Jovian plants may cause an enhancement of stellar magnetic activity. Cuntz et al. (2000) suggested the enhancement of stellar activity caused by the tidal and/or magnetic interactions between a host star and a hot Jupiter. Shkolnik et al. (2003) found that the enhancement of Ca II H and K emission is synchronous with Jupiter’s orbit in HD 179949. A similar synchronized Ca II enhancement were reported in other stars, e.g. $\upsilon$ And, HD 189733, HD 73256, and $\tau$ Boo (e.g. Shkolnik et al., 2005, 2008). Theoretically, Rubenstein & Schaefer (2000) proposed that superflares are caused by magnetic reconnection between fields of a host star and a hot Jupiter. They argued that superflares only happen on those stars hosting close-in giant planets, and our Sun could not exhibit superflares. Through MHD numerical simulations, Ip et al. (2004) suggested that flare-like activity with the energy comparable to that of the typical solar flares could be triggered by interaction of the magnetosphere of close-in giant planets with the coronal magnetic field of a host star.

In this work, we use the Exoplanet Follow-up Observing Program for TESS (ExoFOP-TESS³³3https://exofop.ipac.caltech.edu/tess/) to cross-match flare stars in our database. Only three targets (0.96%) of 311 flare stars may host planets, and 204 stars (0.91%) of 22,228 nonflare stars may host planets. We do not find a significant fraction excess of flare stars hosting planet candidates compared with nonflare stars. Besides, after adding 400 flare stars and 25,334 nonflare stars from TESS first-year observation (Tu et al., 2020), these two fractions are 0.89% for all 711 flare stars and 0.84% for all 47,563 nonflare stars, respectively. The above comparisons support that the idea that star-planet interactions may not be the main mechanism for generating superflares of solar-type stars.

We get three planet candidates listed in Table 6. From their properties derived by the automatic pipeline of TESS, these two candidates TOI 1254.01 and TOI 1425.01 are up to requirements of hot Jupiter (with orbiting period $P_{\rm orbit}<10\;{\rm days}$, and planetary radius $8.0\;{R_{\oplus}}\leqslant R_{\rm planet}\leqslant 32.0\;{R_{\oplus}}$ (Howard et al., 2012)). Hot Jupiter was believed to be an important factor of generating superflares (e.g. Rubenstein & Schaefer, 2000; Ip et al., 2004), but recent works show this mechanism is rare for superflares on solar-type stars (e.g. Maehara et al., 2012; Shibata et al., 2013; Tu et al., 2020). Besides, TOI 1429.01 is an ultrashort period (UPS) planet candidate, as its orbit period less than 1 day (Winn et al., 2018). A UPS planet is close to its hosting star, so star-planet interaction may evaporate its atmosphere. We only refer to those planetary properties from ExoFOP-TESS, which may be not accurate for further specific analysis, but these three planet candidates are worthy of being observed by further confirmations.

Notsu et al. (2019) used the plot of superflare energy and stellar period and found a decreasing trend between the upper limit of energy and stellar periodicity. This discovery was based on the catalog from McQuillan et al. (2014), which applied accurate period estimations for Kepler targets. In our previous work (Tu et al., 2020), we found that only the tail part (period more than a few days) of the plot showed this trend, as the TESS unique observation mode will limit the estimation of stellar periods.

In this work, from Figure 5(a), the descending tendency is only apparent for stars with periods of more than a few days (tail part). We find no strong decreasing trend with the mean of $\log E_{\rm flare}$ in nine period bins. Besides, the flare frequency is decreasing with periods over a few days, as shown in panel (b). Because the upper energy limit in each period bin from the work of Notsu et al. (2019) is just the maximum flare energy for the star, we use the maximum energy on the star instead of attaching all superflare energies. We also separate stars according to their surface temperatures, as described in the legends of Figures 5 (c) and (d). These two panels both show an upper-limit decreasing trend at the tail parts, as we see in panel (a). Gray squares indicate for the mean value of $\log E_{\rm flare}$ in each period bin, which shows a more apparent descending tendency.

Because TESS changes its observing field every $\sim 27$ days, the targets observed by TESS have different durations of observing time. We still use the same functions as we specifically described in our previous work, which is based on the frequency calculation of Maehara et al. (2012). We also subdivide those targets in our database into different subsets called as Set-$n$, where $n$ equals 1 - 13 and stands for the umber of sectors in which the targets have been observed. The distribution of Set-$n$ is shown in Figure 6. Specifically, counts of superflares, flare stars, and solar-type stars are shown in Table 7 according to their surface temperatures and periods. In Figure 6, we also show the Set-$n$ distribution in our previous work (Tu et al., 2020) as panel (b) for comparison. It seems that the Set-$n$ distributions are similar.

We use the following function to calculate flare frequencies (Tu et al., 2020).

where the subscript $n$ is same as the $n$ of Set-$n$. $N_{\rm flares}$ and $N_{\rm stars}$ represent counts of superflares and solar-type stars in the same Set-$n$. $\Delta E_{\text{flare}}$ gives the energy range in each bin, and

where 22.64 represents the mean of the effective observing time of one sector in the TESS second-year observation, as TESS does not fully observe in 27 days for one sector. The flare frequency is derived as

which represents the mean frequency of 13 Set-$n$.

In Figure 7(a), we compare the flare peak amplitude distribution of these two-year observations. Here $F_{\rm peak}$ is $F_{\rm flare}$ of Equation (4) at peak times. The solid line is higher than the gray dashed line, with peak amplitudes less than $10^{-2}$. Contrarily, the solid line is slightly lower than the gray dashed line, with peak amplitudes greater than $10^{-2}$. The results indicate that the detected superflares in this work generally have a lower $F_{\rm peak}$ than those of Tu et al. (2020). As peak amplitude is positively correlated with flare energy, if more superflares have lower peak amplitudes, it will pull down the overall energy distribution, just like what is shown in Figure 2(d). Meanwhile, the superflares with $F_{\rm peak}\gtrsim 10^{-2}$ are lower than those of previous work, which may even reduce the overall energy of their distribution. From the view of peak amplitude, the superflares in this work are overall weak.

In Figures 7 (b) and (c), we get similar results as in previous literature (e.g. Maehara et al., 2012; Shibayama et al., 2013; Maehara et al., 2015); (1) fast-rotating stars have higher superflare occurrence frequencies than slowly rotating ones, as shown in panel (b), and (2) cooler stars generate superflares more frequently than hotter stars, as shown in panel (c). Besides, in Figure 7(b), dotted line is basically overlapped with the solid line. This indicates that the solar-type stars observed by TESS are dominated by rapidly rotating stars. In this work, 19,101 stars have periods less than 10 days, which is 85% fraction of all solar-type stars. This same number of this fraction is also deduced from TESS first-year observations (Tu et al., 2020). In Figure 7(d), solid black lines are above dashed-dotted gray lines. This may also be explained by the large fraction of rapidly rotating stars in our database. In this work, 85% of the stars have a period of less than 10 days, and this number is 32% in Maehara et al. (2015). That stellar periodicity is related to their age, so young stars rotate faster than old stars (Skumanich, 1972; Barnes, 2003). Furthermore, young stars are more active and have higher superflare frequencies than old stars.

The flare frequencies of Sun-like stars are only deduced from four superflares in panel (d), so we may not be able to statistically get solid conclusions from that.

Flare frequency distribution as a function of energy can be applied by a power-law relation as

Through linear regression, we get a power-law index $\gamma\sim-1.76\pm 0.11$ for all superflares, and $\gamma\sim-2.30\pm 0.66$ for superflares on slowly rotating stars. For all superflares, $\gamma\sim-1.76\pm 0.11$ is lower than that of the TESS first-year observation ($\gamma\sim-2.16\pm 0.10$) and Shibayama et al. (2013, $\gamma\sim 2.2$). We also combine the superflare database from our previous work with data from this work and get $\gamma\sim-1.91\pm 0.07$ for all superflares, as shown in Figure 8(a). For superflares on slowly rotating stars, $\gamma\sim-1.94\pm 0.31$. One thing needs to be addressed: that combining 2 yr superflares may not be a strict behavior, as the TESS pipeline and accuracy of our visual inspection are all improved (refer to Section 2.1). However, we do not apply the new method and TESS newly released data for reprocessing TESS first-year observation, as they may not significantly affect the robustness of conclusions in this work. Anyway, we can still roughly get overall flare frequencies as a function of energies from TESS 2 yr observations.

For comparison with solar flare frequencies, we use different kinds of solar flares from some databases and show them in Figure 8(b). It is obvious and fascinating that, as shown in Figure 8(b), solar nanoflares (Aschwanden et al., 2000), solar microflares (Crosby et al., 1993), solar flares (Shimizu, 1995), and superflares on solar-type stars from Kepler (Shibayama et al., 2013; Maehara et al., 2015), and from TESS (Tu et al., 2020, and this work) basically distribute around $dN/dE\propto E^{-1.8}$ (black dashed-dotted line). This may indicate same mechanism is shared by these different types of flares. As the flare energy increases, their frequency decreases. It is obvious that TESS gets overall higher superflare frequencies than Kepler after we use the combined superflare database from the TESS 2 yr observations. The reason why superflare frequencies from this work (red solid lines) are higher than those of Maehara et al. (2015, gray dash-dotted lines) is addressed in the above paragraph. In total, the superflare frequency for Sun-like stars is apparently higher than that of Shibayama et al. (2013). However, only 11 superflares are observed by TESS 2 yr missions. We expect that more superflares on Sun-like stars can be detected by TESS extended missions, and other practicable methods can more accurately estimate stellar periods for TESS targets, which will strengthen the reliability of selecting superflares on Sun-like stars. Interestingly, similar power-law distributions are also found in different astronomical phenomena, i.e., gamma-ray bursts (Wang & Dai, 2013; Yi et al., 2016; Lyu et al., 2020), fast radio bursts (Wang & Yu, 2017; Cheng et al., 2020), and black holes (Li et al., 2015; Wang et al., 2015; Yan et al., 2018; Yuan et al., 2018).

Statistically, we find that superflares and solar flares may have similar physical mechanisms, as discussed in Section 3.6. Here we also apply power-law correlation fitting between the flare energy and duration. In Figure 9(a), we get the power-law correlation

with $\beta=0.42\pm 0.01$. Here $\beta=1/3$ is derived from the magnetic reconnection theory (Maehara et al., 2015), which is successfully applied to solar flares. Here $\beta=0.42\pm 0.01$ is the same as that of our previous work (Tu et al., 2020). And Maehara et al. (2015) got $\beta=0.39\pm 0.03$. For solar flares, Tu & Wang (2018) got $\beta=0.33\pm 0.001$. So, superflares generally get a higher $\beta$ value than one=third. This reinforces our conclusion that solar flares and superflares may have similar mechanisms not only from their occurrence frequency distributions but also from the close correlation between superflare energies and durations. However, the obviously higher value of $\beta$ may also point out that there must be some differences between mechanisms of solar white flares and superflares on solar-type stars.

If we consider magnetic strength as a variable, the Equation 15 will be reformed as

where $B$ and $L$ are the coronal magnetic field strength and flare length scale, respectively. Namekata et al. (2017) derived these two relations in detail, and we also use the coefficients of $B$ and $L$ also from their work to draw scale lines in Figure 9(b). Solar white-light flares of Namekata et al. (2017), superflares from Maehara et al. (2015), and our previous work are all imported into this figure. Compared with the red dots in the figure, superflares from this work (blue circles) have a larger overlapping area with those (red squares) from Maehara et al. (2015). Even if the TESS filter is less sensitive to weak superflares, as we have discussed in Section 3.1, there are still some less energetic superflares with $E_{\rm flare}\lesssim 10^{34}$ erg, which we have selected through the improved visual inspection. Overall, superflares have a higher coronal magnetic field strength and larger length scale than those of solar white-light flares.

Recently, statistical research on superflares on giant stars has been done by Kővári et al. (2020). They used Kepler and a small part of the TESS data and searched 60 superflares on the star KIC 2852961. It is also interesting to compare superflares on solar-type stars with those on giant star KIC 2852961. These superflares are shown by black squares in Figure 9(b). They seem to be an extension part of the superflares on solar-type stars with greater flare energies and longer durations. From the scaling line of the coronal magnetic field strength ($B$) and flaring length scale ($L$), we find that the superflares on KIC 2852961 basically have a coronal magnetic strength between 30 G and 60 G, which is smaller than that of superflares on solar-type stars. From the view of flaring length scale, superflares on solar-type stars are basically under the limitation of solar diameter ($R_{\rm sun}$). However, superflares on giant stars basically have a much larger $L$ than solar-type stars. The radius of this giant star is reported as $R_{*}=13.1{R_{\odot}}$. According to Figure 9(b), the flaring lengths of superflares on KIC 2852961 are less than the diameter of this star ($R_{\rm star}=2R_{*}$). These results also strengthen the conclusion from Kővári et al. (2020) that the largest flares of KIC 2852961 have flare loop sizes on the order of $R_{*}$.

In Section 2.3 and 2.4, we have introduced the$S$-index and $R_{\rm var}$, including their calculation methods. In this work, we have 7454 solar-type stars matched with spectra obtained by LAMOST. According to their observing time (before or after 2016), we classify these stars into two subsets, as described in the legend of Figure 10(a). From the histograms in Figure 10(a), they do not show obvious differences. The Kolmogorov–Smirnov test also shows that their similarity is over the 99% confidence level. The result indicates that the observation of LAMOST is more and more stable over many years, and we could make full use of LAMOST to push our understanding of our topic addressed in this work. Besides, the $S$-index represents the chromospheric activity of a star, which may not be significantly changed within a couple of years. This idea was proved by Zhang et al. (2020a). They used 577 stars, of which each star matched more than one record of LAMOST data. After analyzing those spectra, they did not find significant variations of the $S$-index of a star within a couple of years. For the above reasons, in this work, we choose to still use LAMOST data from late 2011 to early 2020, even though only 228 (3%) of 7454 spectra are observed by LAMOST after 2019 July, when TESS began its second-year mission.

In Figure 10(a), two ranges of solar $S$-index are imported. The shaded area represents ranges imported from Zhang et al. (2020a), who calibrated the solar $S$-index in the scale of LAMOST. The dashed vertical lines are from Henry et al. (1996), who collected $S$-indexes at solar maximum and minimum solar cycles. From the distribution in Figure 10(a), when we use the upper limit of the $S$-index range, over 52.58% of the stars have larger $S$-indexes than that of the Sun.

In Figure 11, we classify solar-type stars into three subsets according to (1) whether the star generates superflares or not and (2) how much energy the superflare have released. The value $4.68\times 10^{34}$ erg is the mean value ($\mu$) of the normal distribution fitting in Figure 2(a). The correlation between $S$-index and $R_{\rm var}$ is suggested to be true by Notsu et al. (2015) and Karoff et al. (2016). They concluded that the larger $R_{\rm var}$ may indicate that larger magnetic features (e.g. star spots) are covering on the stellar surface; meanwhile, their magnetic activities (indicated by a larger $S$-index) are active. Zhang et al. (2020a) used more than 2000 spectra and concluded that the asymmetry of the magnetic features plays an important role in creating photometric variability. In Figure 11(a), the stars without any superflares (blue circles) vaguely show this correlation, but it is not as obvious as that shown in Notsu et al. (2015) and Karoff et al. (2016). Here the unbalanced distributions of the $S$-index and $R_{\rm var}$ might be the main reason for the seeming dispersion of this correlation. From Figure 11(b) and (c), scatters are more concentrated in the area with an $S$-index ranging from 0.2 to 0.3 and $R_{\rm var}$ ranging from $10^{-3}$ to $10^{-2}$. Besides, we have 7375 nonflare stars in our database, which is more than the 1400 scatters in Figure 6 of Karoff et al. (2016). For flare stars (red and yellow circles), the correlation is much more obvious. As we can still see the positive tendency of this correlation from nonflare stars, we think the magnetic features (e.g. star spots) might be the main reason for the overall positive trend of this correlation. And in small $S$-index intervals, dispersion of the $R_{\rm var}$ value may be caused by asymmetry of the magnetic features on the stellar surface. Moreover, it is also reasonable to consider inclination angle between the rotation axis and observing direction (e.g. Notsu et al., 2013; Ikuta et al., 2020) as another kind of asymmetry to affect $R_{\rm var}$. Besides, we also see from Figure 11(a) that those flare stars are not with $R_{\rm var}\gtrsim 10^{-1}$, while their $S$-indexes become larger. Especially for those flare stars with $S\gtrsim 0.4$, unlike what is expected according to the positive correlation between $R_{\rm var}$ and $S$-index, we do not see larger $R_{\rm var}\gtrsim 0.1$. We speculate that there seems to be a saturation-like feature at $R_{\rm var}\sim 0.1$. However, as there are only 79 flare stars that do not make this saturation-like feature statistically convincing, we will specifically discuss this saturation-like feature from Figure 12 in the following.

In Figures 11 (b) and (c), as only 79 superflare stars have been observed by LAMOST, fitting the normal distribution to the two subsets may not be rational. We use the mean and median values of their distributions to represent the overall level of distribution, which is shown as short straight solid or dashed lines. Specially, from the $S$-index distribution, flare stars basically have higher $S$-indexes than nonflare stars. Red histograms have higher mean and median values of the $S$-index than yellow histograms, which indicates that flare stars with higher flare energies ($E_{\rm flare}>4.68\times 10^{34}$ erg) tend to possess greater chromospheric activity than those flare stars with lower flare energies ($E_{\rm flare}<4.68\times 10^{34}$ erg). Apart from that, the majority of flare stars have more active chromospheric activity than the Sun.

We also see that flare stars have greater photometric variability ($R_{\rm var}$). In order to draw conclusions from more convincing results, we use all solar-type stars in this work to get a distribution of $R_{\rm var}$ as shown in Figure 12(a). The same classification for three subsets in Figure 11 is also practiced. Overall, flare stars generally have larger $R_{\rm var}$ than nonflare stars. In Figure 12(a), we do not use the normal distribution to fit the red histograms, because they do not morphologically show the shape of normal distributions. However, many more flare stars in this subset are concentrated in the higher $R_{\rm var}$ range than yellow histograms. In the panel (b) of Figure 12, we plot the stellar $R_{\rm var}$ versus the maximum flare energy of the star. Gray squares represent mean values of $\log E_{\rm flare,max}$ (maximum flare energy in logarithm scale) in 11 $R_{\rm var}$ bins. From the gray squares, the increasing trend is much more obvious. We also combine 400 superflare stars from our previous work (Tu et al., 2020), and this positive tendency is also shown by gray squares in panel (d) of Figure 12.

Nonflare stars generally have lower values of $S$-index and $R_{\rm var}$ than flare stars, which indicates that their magnetic activity may be less active and may not produce any energetic magnetic features (e.g. larger star spots) to generate superflares (Shibata et al., 2013). For flare stars, their maximum superflare energy is correlated with stellar $R_{\rm var}$, i.e., the greater $R_{\rm var}$ they have, the bigger magnetic features they may host, and the more energetic a superflare can be generated. And larger magnetic features may also increase the probability of the star generating superflares, which may explain the positive correlation between flare frequency for a single star and its $R_{\rm var}$ as shown in Figure 15(a). Besides, the tendency of many more scatters being located on the lower right side of Figure 13 can also be explained by this improved probability. Specific discussion of Figures 13 and 15(a) can be found in Sections 3.9 and 3.10, respectively. Besides, the above results indirectly reinforce that solar-type stars can generate superflares through magnetic activity (Sections 3.6 and 3.7), as long as their magnetic features are active enough. From the view of the possibility for the Sun to generate superflares, the above results enhance the conclusion that those flare stars are much more active than the Sun, as the flare stars in our combined database (711 superflare stars) have larger values of $R_{\rm var}$ than those of the Sun.

From Figure 12(a), we may also see that red histograms saturate around $R_{\rm var}\sim 10^{-1}$. We also draw the line in panel (b) and find that for flare stars, there is an upper limit at $R_{\rm var}\sim 10^{-1}$. Larger flares may be generated from stars with greater $R_{\rm var}$, according to the positive relation between $R_{\rm var}$ and maximum flare energy shown in Figures 12(b) and (d). But those superflares are not generated from stars with $R_{\rm var}\gtrsim 10^{-1}$. The 22,228 nonflare stars (blue histograms) shown in panel (a) also have extremely fewer counts over $R_{\rm var}\sim 10^{-1}$. So, this indicates that for general solar-type stars, the surface magnetic features (e.g. star spots) and their asymmetry distributions may not be able to cause $R_{\rm var}$ greater than $\sim 10^{-1}$. Furthermore, there are still a few nonflare stars with even larger $R_{\rm var}$ than $R_{\rm var}\sim 10^{-1}$, as shown in Figure 11. Potentially, they also generate superflares that are not recorded by TESS. But this possibility may not affect our statistical conclusions. It is also possible that those stars are ellipsoidal binaries, but we think that they will not take a large portion of our data set. From Kepler observations, only about 1.3% of all Kepler targets are ellipsoidal binary systems (Kirk et al., 2016). By the way, from the $R_{\rm var}$ distributions in Figure 12(c), the $R_{\rm var}$ of the flare stars in this work (red histograms) is lower than that of Tu et al. (2020). As we have discussed in Section 3.1, these results may be the reason why the superflare energies of this work are slightly lower than those of Tu et al. (2020), as shown in Figure 2(d).

Under the assumptions that we observe the star with an observing inclination angle of $90^{\circ}$, and the distribution of its magnetic features (e.g. star spots) is extremely asymmetrical, according to our results, it is hard for a solar-type star to be covered by brighter (e.g. faculae) or darker (e.g. star spots) magnetic features with an area over 10% of the stellar surface ($R_{\rm var}\sim 10^{-1}$ and $10^{-1}\times 100$). It must be noted that, if we consider the limb-darkening law and the temperature of stellar spots as factors for measuring photometric variation (Maehara et al., 2020), the actual coverage of magnetic features for some stars may be larger than 10%. Considering the limb-darkening coefficients of Claret (2000), the effective stellar disk is underestimated by around 15% - 17% for the solar-type stars of TESS. Temperature of stellar spots is estimated by correlation with stellar surface temperature (Berdyugina, 2005; Maehara et al., 2020). Then, the saturation-like feature is still apparent but with a little bit higher value around 15%. From this aspect, magnetic features like star spots cannot be subtly generated or congregated even larger. So, it may prohibit the maximum energy of the superflare from getting even bigger. It may be the reason for superflare energies saturating at around $\sim 10^{37}$ erg (Wu et al., 2015). We also get a similar result with a limitation of flare energy around $\sim 10^{36}$ erg from the TESS observation. Superflares from Wu et al. (2015) are systematically with higher energies than ours. Anyway, the saturation of superflare energy may be caused by the coverage limitation of stellar magnetic active regions.

This saturation-like feature of brightness variation amplitude may not be caused by the selection effect of different telescopes. Using Kepler data, Maehara et al. (2012) and Notsu et al. (2013) also found that the upper limit of the brightness variation amplitude of superflare stars is around $10^{-1}$. However, they mainly focused on the scale correlation between brightness variation and flare peak amplitudes and did not discuss this saturation-like feature. Furthermore, Rebull et al. (2016) used K2 data, and also showed this feature for stars from the Pleiades cluster, as shown in Figure 12 of their work. It is necessary to figure out what physical mechanism caused this saturation-like feature. This is also an interesting topic for the study of stellar activities. A bunch of simulation works have focused on reconstructing light curves for observed stars by using artificial functions of star spots and faculae (e.g. Basri & Shah, 2020; Ikuta et al., 2020; I\textcommabelowsık et al., 2020; Shapiro et al., 2020). These simulations are all based on morphologically artificial functions to describe the generation and evolution of star spots and faculae. These functions are generally without any internal physical mechanism. Fang et al. (2016) tried to figure out spot coverage of stars from Pleiades cluster using LAMOST data. Because there are not enough data with successful spot configurations, it is hard to see a clear spot coverage limitation from their results. In the future, we hope numerous observations could statistically configure out the upper limit of star-spot coverage distribution, and numerical simulations could also consider real physical patterns (e.g. Schüssler & Vögler, 2006; Rempel & Cheung, 2014) of star spots to give rational explanations of this saturation-like feature.

Besides, the star TIC 93277807 with $1.77\times 10^{37}$ erg superflares has $R_{\rm var}=0.0735$, which is also less than $10^{-1}$. We show TIC 93277807 as a blue star in Figures 12 (b) and (d). Star TIC 93277807 is far away from the majority of superflare stars. It increases the necessity for analyzing TIC 93277807 with details in the future. Star TIC 93277807 deviates from other flare stars, and its energy is almost more than 1 order of magnitude higher than the others. This superflare can be treated as an extreme event. According to the Dragon-King hypothesis (Sornette, 2009; Sornette & Ouillon, 2012), this kind of event may belong to another population, rather than the general distribution, where relatively small events (e.g. the majority of superflares) belong. The extreme event may be generated by different mechanisms, e.g. processes of synchronization and amplification. Aschwanden (2019) found that Dragon-King events are rare in solar and stellar flares. From their results, a single model is enough to describe small, large, and extreme flares. However, they also believed that, for solar and stellar flares, processes of synchronization and amplification may not be possible. In the future, it will also be attractive to solve what caused the most enormous superflare on TIC 93277807.

Apart from the above discussions from Figure 12 of maximum superflare energy versus $R_{\rm var}$, we also make a diagram of energy and $R_{\rm var}$ for all superflares. In Figure 13, 1216 superflares from the TESS first year observation (Tu et al., 2020) are represented by blue circles, and red circles stand for the 1272 superflares in this work. Comparing with Figures 12(b) and (d), the saturation-like feature of $R_{\rm var}\sim 0.1$ is still prominent in Figure 13. The superflares from one solar-type star with relatively smaller energy are vertically distributed under the maximum superflare of the star. For TIC 93277807, there were four superflares searched last year and represented by blue stars in this diagram. The second-largest superflare of TIC 93277807 is still away from the majority superflares of other solar-type stars. This result makes this star and its superflares more intriguing to deeply study in the future. We also make a diagram of maximum superflare energy versus $S$-index in Figure 14. As $R_{\rm var}$ may be influenced by an inclination effect, the $S$-index seems to be a better indicator of stellar activity. But, from Figure 14, we do not find any apparent correlation between $S$-index and maximum superflare energy. Only the upper limits of maximum flare energy from those scatters with $S\lesssim 0.4$ may be positively correlated with the $S$-index. A larger $S$-index indicates a more active stellar chromosphere, which may increase the possibility for the star to generate superflares. More apparently, this opinion can also be supported by the positive correlation between flare frequency for a single star and $S$-index, as shown in Figure 15(e). However, there are only 79 superflare stars observed by LAMOST that may not give statistically convincing results. In the future, we expect that more spectral information on superflare stars can be observed by sky survey projects to test these results.

Some works have studied the correlation between star spots and energies of superflares (e.g. Shibayama et al., 2013; Maehara et al., 2017; Notsu et al., 2019; Okamoto et al., 2021). Through assuming the relation between the size of the star spot and superflare energy, their works indicated that the same physical processes are found from stellar superflares and solar flares. Here we try to deduce a simple correlation between superflare energy and $R_{\rm var}$. The variation of normalized stellar light curves and size of star spots can be connected by

where $\Delta{F}/{F}$ describes the amplitude variation of the star (Notsu et al., 2019, etc.), which can also be replaced by the photometric variability ($R_{\rm var}$). Here $A_{\rm star}$ is the apparent stellar size ($\pi{R_{*}}^{2}$), and $T_{\rm spot}$ and $T_{\rm star}$ are the temperatures of the stellar surface and star spot, respectively. The energy of the superflare can be derived by

where $L$ and $B$ are the size and magnetic field strength of the star spots, respectively (Notsu et al., 2019, etc.), and $f$ is the fraction of magnetic energy ($E_{\rm mag}$) that is released through superflares.

From research of solar eruptive events, Emslie et al. (2012) found that the bolometric energy of solar flares takes, on average, $\sim 6\%$ of the magnetic energy in the active region. Typically, $f$ is usually assumed to be on the order of $10\%$ (e.g. Shibata et al., 2013; Maehara et al., 2015; Notsu et al., 2019). In this work, we also take $f=0.1$ in Equation (19). The magnetic field of spots in Equation (19) is absolutely different from the coronal magnetic field strength in Section 3.7 and Figure 9(b). Here we take $B=1000\,{\rm G}$ according to the idea that the typical magnetic field strength of sunspots is on the order of $1000\,{\rm G}$ (e.g. Shibata et al., 2013; Maehara et al., 2015; Notsu et al., 2019). We also add $B=3000\,{\rm G}$ as an upper-limit magnetic field strength of spots (Notsu et al., 2019). As our database only includes solar-type stars, we assume $A_{\rm star}=\pi{R_{\odot}}^{2}$, $T_{\rm star}=5800$ K (solar surface temperature), and $T_{\rm spot}=4000$ K in Equation (18). By considering the inclination effect (with inclination angle $i$; Notsu et al., 2013), the correlation of $E_{\rm flare}$ and $R_{\rm var}$ can be written as

These correlations ($E_{\rm flare}\propto R_{\rm var}^{3/2}$) are shown by red and blue dashed lines (with inclination angles labeled by $i=2^{\circ}$ and $90^{\circ}$) in Figures 12(b) and 12(d) and 13.

From Figures 12(d) and 13, there are many scatters obviously located on the lower right side of the plot, which is also suggested by Notsu et al. (2019). This tendency can be explained by the idea that a greater $R_{\rm var}$ indicates bigger magnetic features on the star, which may increase the probability for the star to generate superflares as discussed in Section 3.8. As shown in Figure 15(a), the correlation between flare frequency for a single star and $R_{\rm var}$ may also support this opinion. Besides, from Figures 12(b) and (d) and 13, almost all superflares are below the line of $i=2^{\circ}$ ($B=1000\,{\rm G}$). The case of $i=2^{\circ}$ is where we observe the star from a nearly stellar pole-on direction, which is the upper limit of Equation (20). So, from the above results, we confirm that superflare energy is basically released by the energy stored around the star spots. These results strengthen our conclusions in Section 3.8 that the connection between magnetic features (e.g. large star spots) and superflares is subsistent, and they also strengthen the consistency of stellar superflares and solar flares. Similar discussions can also be found in Section 3.6 and 3.7.

Here we study the flare frequency for a single star ($f_{*}$ from Equation (10)), and try to explore the relations between $f_{*}$ and other stellar properties. Due to the unique observing mode of TESS, stars are not observed with the same time interval. As we use the effective observing time to roughly estimate $f_{*}$ through Equation (10), and many stars are just observed by one sector (Figure 6), we just show statistically obtained relations and try to connect with empirical theorems of stellar evolution. We also use the mean flare frequency in each bin of different properties to depict their relations, as shown by gray squares in each panel of Figure 15.

In Figure 15(a), those flare stars have a relative higher $R_{\rm var}$ in the range between $10^{-3}$ and $10^{-1}$. The mean $f_{*}$ is significantly increased with $R_{\rm var}$ higher than $10^{-2}$. Apart from that, in Figure 15(e), for those superflare stars with $S$-indexes, the positive relation seems more obvious, as shown by gray squares. The distribution of these flare stars are also indicating that stars with higher $S$-indexes generate superflares more frequently. As discussed in Sections 3.8 and 3.9, these results reinforce that magnetic activities are positively correlated with the generation of superflares.

For Figures 15 (b)-(d), the maximum energy ($E_{\rm flare}$), duration, and flare peak amplitude ($F_{\rm peak}$) of superflares are all collected. Significant positive relations between these three parameters and $f_{*}$ can be found not only from the gray squares but also from their scatter distribution. Statistically, if the star can generate a much larger superflare, it will frequently produce more relative, smaller superflares, which will apparently increase the flare frequency of the star. So, the flare frequency ($f_{*}$) for single stasr can be positively correlated with these three parameters.

From Figure 15 (f)-(h), we try to explore whether there are relationships between $f_{*}$ and stellar evolution. As solar-type stars evolving, it will rotate slower (Noyes et al., 1984; Barnes, 2003), and become more luminous and hotter (Hurley et al., 2000). We use stellar periodicity, effective temperature and luminosity as indicators of star evolution. In Section 3.6, we have found that slowly rotating stars have lower flare frequencies than rapidly rotating stars. In Figure 15(f), a similar result is also shown by gray squares and scatter distribution gathering in the short-period range. So, it also indicates that young stars are rapidly rotating stars, more active and with higher flare frequency. From the gray squares in panels (g) and (h), we also get slightly negative correlations between $f_{*}$ and surface temperature and luminosity. Compared with less evolved stars, it seems that the stars with larger $T_{\rm eff}$ and $L_{*}$ generate superflares less frequently. From the scatter distribution, there are some hotter and more luminous stars with greater $f_{*}$ from scatters in panels (g) and (h) of Figure 15, which makes the negative correlation less obvious than that in panel (f).