Statistical Properties of Superflares on Solar-type Stars: Results Using All of the Kepler Primary Mission Data

We searched for superflares using the Kepler 30 minute (long) time cadence data (Koch et al. 2010) that were taken from 2009 May to 2013 May (quarters 0–17). We retrieved the data of the Kepler Date Release 25 (DR25; Thompson et al. 2016) from the Multimission Archive at the Space Telescope (MAST). We selected solar-type (G-type MS stars) stars based on the evolutionary state classification and effective temperature (T_eff) values listed in Berger et al. (2018). Berger et al. (2018) is the catalog of the 177,911 Kepler stars combining new stellar radius estimates (R_Gaia) from Gaia-DR2 parallaxes with the stellar parameters from the DR25 Kepler Stellar Properties Catalog (DR25-KSPC; Mathur et al. 2017), which incorporates the revised method of T_eff estimation (Pinsonneault et al. 2012). They reported the evolutionary state classifications (MS/subgiants/red giants/cool MS binaries) on the basis of these R_Gaia and T_eff values (see Table 5 of Berger et al. 2018). Using the values in Berger et al. (2018), we identified the MS stars with T_eff = 5100–6000 K as solar-type stars. In total, 49,305 solar-type stars are selected (see item (2) in Table 1). Among these, 11,601 stars with brightness variation amplitude (Amp) and rotation period (P_rot) values reported in McQuillan et al. (2014; see item (3) in Table 1) are finally used for the flare search process.

There are various methods of detecting flares from light curves (e.g., Shibayama et al. 2013; Davenport et al. 2014; Günther et al. 2020). In this study, we follow the flare detection method that is based on our previous studies (Maehara et al. 2012 & 2015; Shibayama et al. 2013), but we include some updates by referring to the methods in several studies (e.g., Gao et al. 2016; Yang et al. 2017). The advantage of the following method is that we do not need to assume any models of flare light-curve shapes to detect flares. Because of this, it is expected that we can detect most of the flares even if they have complicated shapes (e.g., a flare with multiple peaks) in the light curve.

Stellar rotation, which causes periodic variation in light curves, can prevent us from adequately detecting smaller flares (see Yang et al. 2017; Davenport et al. 2020; Feinstein et al. 2020). In particular, this effect becomes larger for detecting flares on rapidly rotating stars (P_rot ∼ 1 day), since the timescales of flares approximately equal the rotation periods. In order to reduce this rotation effect, we use a high-pass Butterworth filter. The passband edge frequency (W_p) is 2 × 1/P_rot, and the stop-band edge frequency (W_s) is 3 × 1/P_rot. If P_rot is smaller than 0.2 days, we set W_p = 2 × 1/0.2 and W_s = 3 × 1/0.2 day⁻¹, because the high-pass filter does not work. After the light curves are corrected with the high-pass filter method, we calculate the distributions of brightness variations between all pairs of two consecutive data points, as done in our previous studies (e.g., Shibayama et al. 2013). We use this distribution for the flare detection. The threshold of the flare detection was determined as in the to be three times the value at the value of the top 1% of this distribution (1% × 3). In more detail, first we calculated the value A with this criteria: if flux difference is larger than A, the data point is included in the top 1% distribution. Then, we define the flare detection threshold as A × 3. This threshold was chosen as a result of test runs so as not to falsely detect other brightness variations. We also use this to define the end time of the flare.

Figure 1 shows an example of our flare detection method using the light-curve data of KIC 003836772. The rotation period of this star is ∼0.6 days. Figure 1(a) is an original light curve without the high-pass filter, and panel (b) is the distribution of brightness variations between consecutive data points. Figures 1(c) and (d) are the same as panels (a) and (b) but with the high-pass filter. Red arrows (labeled “threshold”) in panels (b) and (d) show the threshold (1% × 3) in the distribution, and they are also shown as red bars in panels (a) and (c). Black arrows (labeled “flare”) in panels (b) and (d) indicate the maximum flux difference between two consecutive points during the flare shown in panels (a) and (c). Although the flare amplitudes in panels (a) and (c) are almost equal, the “threshold” value in panel (c) is smaller than that in panel (a), thanks to the high-pass filter. Without the high-pass filter, this flare is not detected as a flare, since the “threshold” value (red arrow in panel (b)) is smaller than the maximum flux difference during a flare (black arrow in panel (a)). As shown in this example, this high-pass filter method can help us to not overlook or falsely detect flares. However, since the rotational modulations are not completely removed, the detection threshold of flares can still be affected by the rotational modulations. In particular, rapidly rotating stars tend to have larger effects, and we discuss this point again in Section 3.6.

Zoom In Zoom Out Reset image size

The flare start time is defined as the time 30 minutes before the flux first exceeds the detection threshold, considering that the time cadence of the Kepler data used here is 30 minutes. We removed long-term brightness variations around the flare by fitting with the spline curves of three points in order to determine the flare end time. Three points are determined by the data points just before the beginning of the flare and 5 and 8 hr after the peak of the flare. The flare end time was determined by the time the flux residuals become smaller than the standard deviation of the distribution (σ_diff). Two consecutive data points must exceed the threshold so that the event is defined as a flare. We also excluded flare candidates with a shorter decline phase than the increase phase. This is because most light curves of general stellar flares in the optical show the impulsive brightness increase phase and subsequent, more gradual decrease phase (e.g., the “classical flares” shown in Figure 4 of Davenport et al. 2014). It has been reported that short-duration near-ultraviolet (NUV) flares show different shapes of light curves (e.g., flares whose brightness increase time is comparable to or shorter than the decay time), but the total duration is less than few minutes (Brasseur et al. 2019). We can assume it is almost impossible to detect these short-duration flares with the Kepler 30 minute time cadence data used in this study, though we need more studies on the detailed physical properties of these short-duration flares. Because of these things, we exclude flare candidates with a shorter decline phase than increase phase, considering that these events are caused by instrumental effects or unknown astrophysical events different from flares on solar-type stars.

The pixel scale of the Kepler CCDs is around 4″, and the typical photometric aperture for a 12 mag star is about 30 pixels (Van Cleve & Caldwell 2016). This brightness variation of neighboring stars can affect the flux of the target star, and we removed this effect as done in Shibayama et al. (2013) and Maehara et al. (2015). If a flare star is near the target and large-amplitude flares occur, Kepler will detect fake flares on the target star. To eliminate such events, we chose stars without neighbor stars within 12″ for the analysis using the Kepler Input Catalog. We also excluded pairs of flares that occurred at the same time, and the distance between the flares is less than 24″. By this selection, ∼64% of the detected superflares were excluded. We note here as a caution that it is still possible that there remain some faint stars (>19–20 mag, fainter than the typical limiting magnitude of the Kepler Input Catalog), even after this selection process. After eliminating the candidates of flares on neighboring stars, we checked the light curves of the flare candidates by eye. Finally, we checked the flux-weighted centroid of the flares in order to investigate whether the flux-weighted centroid in the quiescent phase and the flare phase are the same. If the difference of the flux-weighted centroid at the time of the flare is 3 times larger than the standard deviation of the difference of 10 moving averages, we judge that the flare does not occur on the target star. We compared the periodic flux-weighted centroid moving with periodic light-curve amplitudes of the targets in order to eliminate binaries. If the correlation coefficient of consecutive data points of brightness variation and the flux-weighted centroid is more than 0.5, we eliminated the targets from our flare detection.

We use the methods described in Shibayama et al. (2013) to estimate the total energy of each flare. We assumed that the spectrum of white-light flares can be described by blackbody radiation with an effective temperature (T_flare) of 10,000 K (Mochnacki & Zirin 1980; Hawley & Fisher 1992). Assuming that the star is a blackbody radiator, the relationship between the bolometric flare luminosity (L_flare), the flare effective temperature (T_flare), and the area of the flare (A_flare) is written as the following equation with the Stefan–Boltzmann constant (σ_SB):

$$\begin{eqnarray}&&{L}_{\mathrm{flare}}={\sigma }_{\mathrm{SB}}{T}_{\mathrm{flare}}^{4}{A}_{\mathrm{flare}}.\end{eqnarray} \tag{ 1 }$$

We estimate A_flare with the observed luminosity of the star (${L}_{\mathrm{star}}^{{\prime} }$) and flare (${L}_{\mathrm{flare}}^{{\prime} }$) using Equations (2)–(5) of Shibayama et al. (2013). The total bolometric energy of the superflare (E_flare) is written as an integral of L_flare(t) during the flare duration,

$$\begin{eqnarray}&&{E}_{\mathrm{flare}}={\int }_{\mathrm{flare}}{L}_{\mathrm{flare}}(t){dt}.\end{eqnarray} \tag{ 2 }$$

This estimate of E_flare is affected by several types of uncertainties. Errors of stellar effective temperature (T_eff) and stellar radius (R) affect E_flare values, and these are typically ∼3% and ∼7%, respectively (Berger et al. 2018). The determination errors of flare start/end points and quiescent levels also have effects on flare amplitude values and, consequently, E_flare values (typically ∼30%). If the blackbody temperature of a flare emission changes to 6000–7000 K, as in solar white-light flares, the flare energy can change by a factor of 0.5 (Namekata et al. 2017). Moreover, since both stellar quiescent radiation and flare emission may not be complete blackbody radiation, E_flare values may have an error of a few tens of percent. In total, E_flare values can have error values of a few tens of percent. We must note these potential error values when we discuss the exact E_flare values reported in the following sections. However, the overall statistical discussions of the relations between superflares and stellar properties (e.g., rotation period, starspot coverage) do not change, since the errors are smaller than 1 order of magnitude.

We detected 2344 superflares on 266 solar-type stars from Kepler 30 minute (long) time cadence data of ∼4 yr (quarters 0–17). Among them, we detected 29 superflares on 16 Sun-like stars (T_eff = 5600–6000 K and P_rot > 20 days). The number of superflares (N_flare), superflare stars (N_flarestar), and all stars analyzed (N_star) are listed in Table 1. The sample size (the number of stars × observation period) of solar-type stars in this study is about four times larger than that of Notsu et al. (2019), and the sample size on Sun-like stars is about 12 times larger. The number of superflares on Sun-like stars detected in this study is >12 times larger than Notsu et al. (2019), since we improved the methods to detect superflares (Section 2.1). We note again that our previous study, Notsu et al. (2019), did not include stars newly identified as solar-type stars in the updated Gaia-DR2 catalog (Berger et al. 2018). Notsu et al. (2019) applied Gaia-DR2 stellar radius values only to the stars that had been identified as solar-type stars in Shibayama et al. (2013), which used the original Kepler input catalog (Brown et al. 2011) for selecting solar-type stars. The effective temperature in the original Kepler input catalog (Brown et al. 2011) is about 250 K lower than that in the Gaia-DR2 catalog (Berger et al. 2018) for most solar-type stars. This is the main reason why the sample increase of Sun-like stars (T_eff = 5600–6000 K and P_rot > 20 day) is 12 times, while that of solar-type stars (T_eff = 5100–6000 K) is only four times.

The stellar and flare parameters of all Sun-like superflare stars detected in this study are listed in Table 2. The data of all 2344 superflares on 266 solar-type superflare stars detected in this study are available in the online-only table. Some example light curves of superflares on the Sun-like stars are shown in Figure 2, and all superflares on Sun-like stars are shown in Appendix A. In addition to superflares on solar-type (G-type MS) stars, we also detected superflares on subgiants, and these are briefly shown and discussed in Appendix C.

Zoom In Zoom Out Reset image size

Table 2.  Parameters of Sun-like Superflare Stars

KIC	Teffa	Ra	Protc	Ampc	mKeplerd	Nflaree	Emaxe	Aspotf	Flagg
	(K)	(R⊙)	(days)	(ppm)	(mag)		(erg)	(1/2 × A⊙)
005648294	5653 ± 198	${1.140}_{-0.082}^{+0.091}$	22.336 ± 0.085	2543.81	14.757	1	2.4 × 1034	4.3 × 10−3	1
005695372	5791 ± 203	${0.886}_{-0.060}^{+0.066}$	21.195 ± 0.130	8881.39	14.054	3	7.7 × 1033	8.9 × 10−3	⋯
0063476560b	5623 ± 197	${0.823}_{-0.057}^{+0.062}$	28.441 ± 0.598	4110.04	14.860	4	2.2 × 1034	3.6 × 10−3	1
006932164	5731 ± 201	${1.038}_{-0.093}^{+0.104}$	22.631 ± 0.077	8538.25	15.883	1	4.0 × 1034	1.3 × 10−2	⋯
007435701	5812 ± 203	${0.809}_{-0.057}^{+0.063}$	20.648 ± 0.053	8039.18	15.090	1	9.1 × 1033	6.7 × 10−3	⋯
007772296	5616 ± 197	${0.897}_{-0.063}^{+0.070}$	23.257 ± 0.083	4136.53	14.880	3	1.1 × 1034	4.3 × 10−3	2
007886115	5729 ± 201	${0.873}_{-0.060}^{+0.066}$	20.351 ± 0.110	5280.17	14.960	1	9.0 × 1033	5.2 × 10−3	⋯
008090349	5760 ± 202	${1.169}_{-0.083}^{+0.092}$	24.898 ± 3.822	2482.89	14.675	1	1.8 × 1034	4.4 × 10−3	⋯
008416788	5889 ± 206	${0.908}_{-0.061}^{+0.068}$	22.519 ± 0.042	3300.19	13.582	5	1.7 × 1034	3.4 × 10−3	⋯
009520338	5947 ± 208	${0.805}_{-0.058}^{+0.064}$	23.857 ± 0.093	4138.47	15.270	1	9.9 × 1033	3.4 × 10−3	1
010011070b	5669 ± 198	${0.770}_{-0.053}^{+0.059}$	23.970 ± 0.623	4423.92	14.949	2	1.6 × 1034	3.4 × 10−3	⋯
010275962	5782 ± 202	${0.789}_{-0.058}^{+0.064}$	26.117 ± 1.426	2323.69	15.296	2	1.9 × 1034	1.9 × 10−3	1
011075480	5953 ± 208	${0.795}_{-0.057}^{+0.063}$	21.651 ± 0.278	2001.46	15.108	1	9.3 × 1033	1.6 × 10−3	⋯
011141091	5756 ± 201	${0.879}_{-0.062}^{+0.068}$	23.446 ± 0.373	2345.32	14.982	1	1.5 × 1034	2.3 × 10−3	⋯
011413690	5886 ± 206	${0.818}_{-0.062}^{+0.069}$	24.238 ± 0.246	3232.87	15.725	1	3.0 × 1034	2.7 × 10−3	1
012053270	5971 ± 209	${0.744}_{-0.053}^{+0.058}$	21.796 ± 1.587	4338.80	15.017	1	3.1 × 1034	3.0 × 10−3	⋯

Notes. The data of the 2344 superflares on the 266 solar-type superflare stars detected in this study are available in the online-only table.

^aEffective temperature and stellar radius values are taken from Gaia-DR2 in Berger et al. (2018). Stellar radii are shown in units of solar radius (R_⊙). ^bThese two stars are also reported in Notsu et al. (2019). See further details in Appendix B. ^cThe rotation periods and average rotational variability are reported in McQuillan et al. (2014). ^dKepler-band magnitude (Thompson et al. 2016). ^eHere N_flare and E_max are the number of superflares and maximum energy of superflares of each Sun-like star, respectively. ^fArea of the starspots on each Sun-like star in units of the area of the solar hemisphere (1/2 × A_⊙ ∼ 3 × 10²² cm²). ^gFlag 1: It is possible that the real rotation period value of the star is different from those reported in McQuillan et al. (2014) and used in this table on the basis of our extra analyses (see Appendix A and Figures 16 and 19 therein). Then the rotation period value of these stars should be treated with caution, and it is possible that some of these stars are not really Sun-like stars, though more investigations are needed (see Appendix A for details). Flag 2: This star is suspected to have a binary (see Figures 17 and 20 in Appendix A). This star is not used in the statistical discussions after Section 3.3.

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

Download table as:  DataTypeset image

It is generally difficult to accurately detect rotation period values of stars as slow as the Sun (P_rot > 20 days), so we double-checked the rotation periods of all of the Sun-like superflare stars (detailed in Appendix A). As a result, some stars may have rotation period values smaller than the values reported in McQuillan et al. (2014). These five stars are flagged as “1” in Table 2, while they are used in the statistical discussions in the following part of this paper because we cannot finally judge only from the periodogram analyses in Appendix A. Flagged as “2” in Table 2, KIC 007772296 is suspected as a binary candidate. The light-curve and periodogram analyses suggest that this star also shows brightness variations with a period of ∼1 day, and this short period might be related to the orbital motion of a binary system (for details, see Appendix A). Because of this, KIC 007772296 is not used in the statistical discussions in the following part of this paper. We also note again that it is still possible that there remain some faint stars (>19–20 mag) in the Kepler photometric aperture (12″). As a result of this paper, we detected 2341 superflares on 265 solar-type stars, and among them, we detected 26 superflares on 15 Sun-like stars (T_eff = 5600–6000 K and P_rot > 20 days). These superflare stars are used in the statistical discussions in the following part of this paper.

Previous survey observations of stellar activity levels (e.g., X-ray-quiescent luminosity, Ca ii H and K index) have shown that the stellar magnetic level decreases as rotation period increases (Noyes et al. 1984; Güdel 2007; Wright et al. 2011). Since stellar age has a strong correlation with rotation, young rapidly rotating stars show higher activity levels, and slowly rotating stars, like the Sun, show lower activity levels. Following this, we expected superflare activities to depend on the rotation period. Notsu et al. (2019) suggested that the maximum superflare energy continuously decreases as the rotation period increases, and there is a difference of roughly 1 order of magnitude between the maximum flare energy on rapidly rotating stars (P_rot ∼ a few days) and slowly rotating stars (P_rot > 20 days). However, a large problem remains: there were some uncertainties on the upper limit of flare energy on Sun-like (T_eff = 5600–6000 K and P_rot > 20 days) stars because the number of Sun-like superflare stars in Notsu et al. (2019) was very small.

We then investigated this relation using a much larger number of superflares newly detected in this study (e.g., an ∼12 times larger sample of Sun-like stars, as described in Section 3.1). Figure 3 shows the relationship between the flare energy (E_flare) and the rotation period (P_rot). The data of the stars with a temperature range (T_eff = 5600–6000 K) close to the solar temperature (T_eff ∼ 5800 K) are plotted with red open squares, while those of late G-type MS stars (T_eff = 5100–5600 K) are blue plus signs. Similar to the results of Notsu et al. (2019), the upper limit of E_flare in each period bin has a continuous decreasing trend with the rotation period. There is at least an order-of-magnitude difference between the maximum flare energy on rapidly rotating stars (P_rot ∼ a few days) and slowly rotating stars (P_rot > 20 days). This decreasing trend, which is confirmed in Notsu et al. (2019) and this study, was not reported in our initial study, Notsu et al. (2013b). This is because some fraction (∼40%) of superflare stars in our initial studies (Maehara et al. 2012; Shibayama et al. 2013; Notsu et al. 2013b) are now found to be subgiants (Notsu et al. 2019), and this contamination of subgiants affected the statistics (see Appendix C for details).

Zoom In Zoom Out Reset image size

In Figure 3, we added the stellar age (t) using the gyrochronological relation (P_rot ∝ t^0.6; Ayres 1997)¹¹ for T_eff = 5600–6000 K in order to compare the age of the superflare stars with that of the Sun (t ∼ 4.6 Gyr). With this scale, as a result of Figure 3, we confirm again that old, slowly rotating Sun-like stars (T_eff = 5600–6000 K, P_rot > 20 days, and age ∼ 4.6 Gyr) can have superflares with ≲4 × 10³⁴ erg, while young, rapidly rotating stars have superflares with up to ∼10³⁶ erg. The scale of t is only plotted in the limited age range t = 0.5–5 Gyr for the following two reasons. (1) As for young solar-type stars with t ≲ 0.5–0.6 Gyr, a large scatter in the age–rotation relation has been reported from young cluster observations (e.g., Soderblom et al. 1993; Ayres 1997; Tu et al. 2015). (2) As for old solar-type stars beyond the solar age (t ∼ 4.6 Gyr), a breakdown of the gyrochronology relations was recently suggested (van Saders et al. 2016; Metcalfe & Egeland 2019).

We note that in Figure 3, there is an apparent negative correlation between the rotation period and the lower limit of the flare energy. However, this trend only reflects the detection limit affected by the rotation period, as discussed in Section 3.6. We also note that in Figure 3, the sample size of the stars in each P_rot bin is different. The number of samples of slowly rotating stars is larger than that of rapidly rotating stars. It is then possible that rare energetic flares approaching the upper limit value for each P_rot bin can be a bit more easily detected on slowly rotating stars than rapidly rotating stars, assuming the power-law flare frequency distributions shown in Section 3.7. This might be related to the tendency that the upper limit of flare energy is roughly constant in the range of P_rot > 10 days. For the same reason, the upper limit E_flare value of rapidly rotating stars (P_rot < 20 days) can be a bit larger (>10³⁶ erg). We will discuss this point in more detail, together with the flare frequency distributions, in Section 3.7.

Most superflare stars show large-amplitude brightness variations, and they suggest that the surfaces of superflare stars are covered by large starspots. Figure 4 shows a scatter plot of flare energy (E_flare) as a function of the spot group area (A_spot) of solar flares and superflare stars. The A_spot values are estimated the same way as in our previous studies (Shibata et al. 2013; Notsu et al. 2019). In these studies, we used the following equation (see Equation (3) of Notsu et al. 2019):

$$\begin{eqnarray}&&{A}_{\mathrm{spot}}=\displaystyle \frac{{\rm{\Delta }}F}{F}{A}_{\mathrm{star}}{\left[1-{\left(\displaystyle \frac{{T}_{\mathrm{spot}}}{{T}_{\mathrm{star}}}\right)}^{4}\right]}^{-1},\end{eqnarray} \tag{ 3 }$$

where A_star, T_spot, and T_star are the apparent area of the star (πR_star²) and the temperature of the starspot and the photosphere, respectively. In this study, T_eff values from Berger et al. (2018) are used for T_star values. The T_spot values are estimated from T_star values by using the empirical relation (Equation (4) of Notsu et al. 2019) deduced from Berdyugina (2005). The total energy released by the flare (E_flare) must be smaller than (or equal to) the magnetic energy stored around starspots (E_mag). Our previous papers (e.g., Shibata et al. 2013; Notsu et al. 2019) suggested that the upper limit of flare energy (E_flare) can be determined by the simple scaling law (see Equation (5) of Notsu et al. 2019):

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{flare}} & \approx & {{fE}}_{\mathrm{mag}}\approx \displaystyle \frac{{B}^{2}{L}^{3}}{8\pi }\approx 7\times {10}^{32}(\mathrm{erg})\\ & & \times \ \left(\displaystyle \frac{f}{0.1}\right){\left(\displaystyle \frac{B}{{10}^{3}{\rm{G}}}\right)}^{2}{\left(\displaystyle \frac{{A}_{\mathrm{spot}}/(2\pi {R}_{\odot }^{2})}{0.001}\right)}^{3/2},\end{array}\end{eqnarray} \tag{ 4 }$$

where f is the fraction of magnetic energy that can be released as flare energy, B and L are the magnetic field strength and size of the spot, and R_⊙ is the solar radius. We note that analyses of large solar flares in Emslie et al. (2012) suggest bolometric flare energies are ∼0.1 of the available magnetic energies (see Figure 4 of Emslie et al. 2012).

Zoom In Zoom Out Reset image size

In Figure 4, almost all the data points of the superflares are below the line of Equation (4) for B = 3000 G, and the upper limit of the flare energies tends to increase as the spot area increases. This means that the superflare is the magnetic energy release stored around the starspots, and the process is the same as that of solar flares. In Figure 5, we also compare the data of Sun-like stars with that of the Sun for reference.

Zoom In Zoom Out Reset image size

We note that the starspot group area estimated here can be smaller than the actual values if the stars have a low inclination angle or have starspots around the pole region (see also Notsu et al. 2015b & 2019 for details). Besides, the real size of the spot group when a superflare occurs can be larger than the size estimated from the brightness amplitude values, since we only use the time-averaged brightness amplitude (Amp) values reported in McQuillan et al. (2014). Some evolution/decay of starspot coverage, which can occur in a timescale of ∼100 days (e.g., Namekata et al. 2019, 2020a), is time-averaged. These can be the reason why some data points of superflares locate above the line of Equation (4) (for B = 3000 G) in Figure 4. As a result, we can conclude that the upper limit of the flare energy is consistent with the magnetic energy stored around the starspots, as also confirmed in Notsu et al. (2019).

In Section 3.2, we confirmed that the upper limit of the superflare energy has a decreasing trend with rotation period. In Section 3.3, we also confirmed that the upper limit of the flare energy is consistent with the starspot coverage, considering that magnetic energy is stored around starspots. Then, how is starspot coverage related to rotation period, and what implications exist for superflare energy?

Our previous paper, Maehara et al. (2017), investigated the statistical properties of starspots on solar-type stars by using the starspot size (A_spot) and rotation period (P_rot) estimated from the brightness variations of Kepler data, and Notsu et al. (2019) updated the results by using Gaia-DR2 stellar radius. The resultant values of A_spot and P_rot are shown in Figure 6, with the increased data of superflare stars newly identified in this paper. As also shown in Notsu et al. (2019), the stars showing more energetic superflares tend to have shorter rotation periods (younger ages) and larger starspot areas.

Zoom In Zoom Out Reset image size

Figure 6 shows that the largest area of starspots on solar-type stars in a given P_rot bin has a roughly constant or very gentle decreasing trend in the period range of P_rot ≲ 12 days (age: t ≲ 1.4 Gyr) for the solar-type stars with T_eff = 5600–6000 K and P_rot ≲ 14 days for those with T_eff = 5100–5600 K. However, in the period range of P_rot ≳ 12 days (age: t ≳ 1.4 Gyr) for the stars with T_eff = 5600–6000 K and P_rot ≳ 14 days for those with T_eff = 5100–5600 K, the largest starspot area of the stars steeply decreases as the rotation period increases. This decreasing trend of the maximum area of the starspots can be related to the decreasing trend of maximum flare energy confirmed in Section 3.2 (see Figure 3). This is because the maximum area of the starspots determines well the upper limit of flare energy, as indicated in Section 3.3. In the case of slowly rotating Sun-like stars with T_eff = 5600–6000 K and P_rot ∼ 25 days, Figure 6(b) shows that the maximum size of the starspots is a few percent of the solar hemisphere. These values correspond to 10³⁴–10³⁵ erg by Equation (4), and the upper limit of the superflare energy of these Sun-like stars in Figure 3 is roughly in the same range.

It is difficult to conclude whether A_spot values are constant or decrease in the short period range because of several factors. As indicated in Section 3.3, the A_spot values estimated from the brightness variations can be smaller than the actual values. For example, the A_spot values estimated from the brightness variations can be smaller than the actual values if the stars have a low inclination angle or starspots around the pole region (see also Notsu et al. 2015b & 2019). In addition, more active stars, such as rapidly rotating stars, can have multiple starspots on the surface, and this also causes the A_spot estimated from the brightness variations to be smaller than the total spot coverage of the star (Namekata et al. 2020a).

However, there is a difference in a bit stricter sense between the decreasing trends of the maximum superflare energy in Figure 3 and the maximum area of the starspots in Figure 6. The maximum superflare energy continuously decreases as the rotation period increases (the star becomes older) in Figure 3, but the maximum area of the starspots does not show such a continuous decreasing trend in Figure 6. On the other hand, the maximum area of the starspots is roughly constant or very gently decreases in the short period range (P_rot ≲ 14 days), but it steeply decreases as the period increases in the longer range (P_rot ≳ 14 days). This difference was also suggested in our previous paper, Notsu et al. (2019), but it was not as clear, since the number of superflare events was small. If this difference between the decreasing trends in Figures 3 and 6 is true, there can be a possibility that the flare energy is determined not only by the starspot area but also by other important factors, though the starspot area is still a necessary condition to determine the flare energy (see Section 3.3). Considering the correlation between the flare activity and the magnetic structure of sunspot groups (see Sammis et al. 2000; Toriumi & Wang 2019), one of the possible factors might be the effect of the magnetic structure of starspots. More complex spots can generate more frequent and energetic flares according to solar observations. If the magnetic structure (complexity) of spots also has a correlation with the rotation period, the upper limit of the flare energy can depend on rotation period, even if the starspot size is roughly constant. We need to conduct more detailed studies on starspot properties to clarify such possibilities (e.g., observations using exoplanet transits in Namekata et al. 2020a).

In addition, as suggested in Notsu et al. (2019), the constant and decreasing trends of maximum starspot coverage can be compared with the relation between X-ray flux and rotation period (e.g., Wright et al. 2011). The X-ray fluxes of solar-type stars are also known to show the constant regime (or so-called “saturation” regime) in the period range of P_rot ≲ 2–3 days, but they decrease constantly as the P_rot values increase in the range of P_rot ≳ 2–3 days. The changing point of this X-ray trend (P_rot ≲ 2–3 days) is different from that of maximum spot size values (P_rot ≲ 12–14 days). These similarities and differences can be interesting and helpful when considering the relation between stellar activity (including starspots, flares, and X-ray steady emissions) and rotation period in more detail, though a detailed study of this point is beyond the scope of this paper.

In Appendix B of Notsu et al. (2019), we investigated the potential differences between the number of stars with P_rot values in McQuillan et al. (2014) and the number of the “real” sample of Kepler field stars by using the gyrochronological relation. As also shown in Table 3, ∼69% (=(20,893–6527)/20,893) and ∼82% (=(28,412–5064)/28,412) of the solar-type stars with T_eff = 5100–5600 and 5600–6000 K, respectively, have no P_rot values in our sample, and they are not plotted in Figure 6. This is because the brightness variation amplitudes of these “inactive” solar-type stars are smaller than the detection limit. We need to pay attention to biases caused by these “inactive” solar-type stars when discussing the relation between the rotation period and superflare properties, such as the occurrence frequency of superflares. The number of stars having P_rot values in each P_rot bin of Table 3 does not show the actual P_rot distribution of the Kepler field. This is because the “inactive” stars with no P_rot values are expected to be dominated by old, slowly rotating stars, and this can be a big problem, especially for discussing the flare frequencies of slowly rotating stars.

In order to roughly evaluate the potential differences caused by the above points, our previous paper, Notsu et al. (2019), estimated the number fraction of solar-type stars in specific P_rot bins by using the empirical gyrochronology relation (Ayres 1997; Mamajek & Hillenbrand 2008). The number of stars with P_rot values larger than P₀ (N_star(P_rot ≥ P₀)) can be estimated from the duration of the MS phase (τ_MS), the gyrochronological age of the star (t_gyro(P₀)), and the total number of stars (N_all):

$$\begin{eqnarray}&&{N}_{\mathrm{star}}({P}_{\mathrm{rot}}\geqslant {P}_{0})=\left(1-\displaystyle \frac{{t}_{\mathrm{gyro}}({P}_{0})}{{\tau }_{{MS}}}\right){N}_{\mathrm{all}}.\end{eqnarray} \tag{ 5 }$$

We assumed that the star formation rate around the Kepler field has been roughly constant over τ_MS.

Using Equations (12)–(14) of Mamajek & Hillenbrand (2008), we roughly estimated the ages (t_gyro(P₀)) of solar-type stars with T_eff ∼ 5800 and 5350 K (B − V ∼ 0.65 and 0.80 from Equation (2) of Valenti & Fischer (2005)) and P_rot ∼ 5, 10, and 20 days. Using these t_gyro(P₀) values and Equation (5), we also estimated the fraction of stars as listed in Table 3. We assume that N_star(P_rot ≥ 20 days) is the same as N_star(P_rot = 20–40 days) when we compare it with N_P(P_rot = 20–40 days) in Table 3. This is because, according to the gyrochronology relation, the rotation period of solar-type stars cannot be larger than 40 days in the duration of the MS phase, and the relation has a breakdown in older, slowly rotating solar-type stars (e.g., t > 5.0 Gyr; van Saders et al. 2016; Metcalfe & Egeland 2019).

As seen from Table 3, there are differences between the number fractions of slowly/rapidly rotating stars in the Kepler sample from McQuillan et al. (2014; N_P) and those estimated from the gyrochronology relation (N_star). In the case of Sun-like stars with T_eff = 5600–6000 K, the fraction of stars with P_rot = 20–40 days among all sample stars has roughly a factor of 2 difference: N_P(P_rot = 20–40 days)/N_P,all ∼ 32% and N_star,all(P_rot = 20–40 days)/N_star,all ∼ 70%. This means that the flare frequency value of Sun-like stars (T_eff = 5600–6000 K and P_rot > 20 days) can become a factor of 2 smaller than the real value if we only consider the Sun-like stars in McQuillan et al. (2014). We have already shown these differences in Appendix B of Notsu et al. (2019) only as a caution. In this study, we consider the differences when discussing flare frequencies in the following sections. Through this, we aim to discuss more accurate values of flare frequencies, especially for slowly rotating Sun-like stars.

Before applying these “gyrochronology corrections” to flare frequency discussions in the following section, we note several possible errors that the gyrochronology can have. One example is that the age–rotation relation (gyrochronology relation) of young (t ≲ 0.5–0.6 Gyr) and old (t ≳ 5.0 Gyr) solar-type stars can have a large scatter (e.g., Soderblom et al. 1993; Ayres 1997; Tu et al. 2015) and breakdown (van Saders et al. 2016; Metcalfe & Egeland 2019), respectively, as mentioned in Section 3.2. We also assumed that the star formation rate around the Kepler field has been roughly constant over τ_MS in the above estimation, but this assumption is not necessarily correct. We should keep these potential errors in mind.

As shown in our previous study, Shibayama et al. (2013), the threshold to detect superflares used in our method (see Section 2.1) tended to be larger in rapidly rotating stars than in slowly rotating stars. We used the high-pass filter described in Section 2.1 to eliminate this tendency. In this section, we evaluate how this tendency is changed by the filter and how much this tendency remains in the filtered data to discuss the relations between rotation period and flare frequency in the following sections. As described in Section 2.1, a flare is detected if the assumed flare amplitude (FAmp_assumed) is larger than the flare detection threshold (Threshold). The Threshold value is calculated as three times the value at the value of the top 1% of the brightness distribution of each star (See Figure 1). Here we check how large-amplitude flares can be detected (FAmp_assumed ≥ Threshold) for each star. We count the number of stars that satisfy FAmp_assumed ≥ Threshold for several period ranges (P_rot = 0–5, 5–10, 10–20, and 20–40 days) as a function of FAmp_assumed: ${N}_{P}^{{{F}{\rm{Amp}}}_{\mathrm{assumed}}\geqslant {Threshold}}$. Then we can define the detection completeness (DC_filter) as

$$\begin{eqnarray}&&{{DC}}_{\mathrm{filter}}=\displaystyle \frac{{N}_{P}^{{{F}{\rm{Amp}}}_{\mathrm{assumed}}\geqslant {Threshold}}}{{N}_{P}},\end{eqnarray} \tag{ 6 }$$

where N_P is the total number of solar-type stars having P_rot values in the considered P_rot and temperature range shown in Table 3. Here DC_filter means the fraction of stars from which flares can be surely detected, since the flare amplitude is larger than the detection threshold (FAmp_assumed ≥ Threshold).

Figure 7 shows the relationship between DC_filter and the flare amplitude for several period ranges (P_rot = 0–5, 5–10, 10–20, and 20–40 days) for the filtered light-curve data (see Figures 1(c) and (d)). There is a tendency for the detection completeness DC_filter to become larger as the flare amplitude becomes larger. However, the tendencies are different between rapidly rotating (P_rot < 5 days) and slowly rotating (P_rot ≥ 5 days) stars. As for the stars with P_rot ≥ 5 days, the detection completeness DC_filter becomes almost ∼1.0 at the flare amplitude of ∼10⁻². This means that we can detect almost all of the superflares on the stars with P_rot ≥ 5 days, if the flare amplitudes are ≥ 10⁻². We note that this flare amplitude value of 10⁻² corresponds to the flare energy with 10^34.5 erg in the case of typical duration superflares in our sample. In contrast, for the rapidly rotating stars with P_rot < 5 days, the detection completeness DC_filter is still ∼0.7 at a flare amplitude of ∼10⁻². This means that we miss 30% of superflares with an amplitude value of 10⁻² occurring on stars with P_rot < 5 days, while we can detect 70% of superflares with an amplitude value of 10⁻² occurring on stars with P_rot < 5 days.

Zoom In Zoom Out Reset image size

This difference means that it is difficult to detect smaller superflares on rapidly rotating stars and corresponds to the apparent negative correlation between the rotation period and the lower limit of the flare energy seen in Figure 3. This difference also can be explained by how rapidly rotating stars tend to have a larger amplitude of the brightness variations (e.g., Figure 6) and the threshold values for the flare detection can be larger. In addition to Figure 7, using the filtered light-curve data, we also calculated DC_filter values for the original light-curve data without the high-pass filter process (see Figures 1(a) and (b)), DC_original. Then ΔDC is defined as the difference between the detection completeness from the high-pass filtered data (DC_filter) and without a high-pass filter (DC_original). As shown in Figure 8, the high-pass filter has caused the detection completeness to be improved by a few percent, especially for the rapidly rotating stars (P < 5 days), though large effects of the brightness variations remain, as already seen in Figure 7.

Zoom In Zoom Out Reset image size

We take into account these dependences of the detection completeness (or “missing rate of flares”) on the rotation period and flare amplitude when we discuss the relationship between the rotation period and superflare frequency in Section 3.7.

We then discuss the flare frequency distributions and their relation to rotation period by using the updated sample of superflares from the Kepler 4 yr data (quarters 0–17) in this study. Figures 9 and 10 represent the occurrence frequency distributions of superflares on solar-type stars with T_eff = 5600–6000 and 5600–6000 K, respectively, for the P_rot ranges of <5, 5–10, 10–20, and 20–40 days. These occurrence rate values of superflares are estimated for each P_rot bin from the number of observed superflares and stars and the total length of the observational period. The potential effects of the gyrochronology correction (Section 3.5) and the flare detection completeness correction (Section 3.6) are not taken into consideration in Figures 9(a) and 10(a). These are the same definitions of flare frequencies as in our previous studies (e.g., Shibayama et al. 2013; Notsu et al. 2019). In Figures 9(b) and 10(b), the frequency values are calculated by incorporating the effect of the gyrochronology correction (Section 3.5) and the flare detection completeness correction (Section 3.6). In order to calculate the frequency of superflares incorporating the effect of gyrochronology correction, we used N_star instead of N_P in Table 3 as a value of the number of observed stars in each P_rot bin. Through this, we aim to discuss more accurate values of flare frequencies, especially for slowly rotating Sun-like stars. In order to calculate the frequency of superflares incorporating the effect of the flare detection completeness correction, we correct the number of superflares in each energy and P_rot bin by using the relationship among DC_filter (detection completeness), flare amplitude, and P_rot shown in Figure 7. From the value of DC_filter in Figure 7, we can estimate how many flares are missed for each P_rot and E_flare (flare energy) range. When we estimate the flare frequency, we correct the number of superflares in each P_rot and E_flare by replacing the number of flares with “(the number of flares)/DC_filter.” In this process, we assume the average relationship among flare amplitude, flare duration, and E_flare for solar-type stars (see energy versus duration, discussed in Maehara et al. 2015). It is assumed that the first data point is the flare with the same value as the threshold, the second point is half of the threshold, and the flare does not continue longer than 1 hr. This does not include various flares, such as long or small flares or superflares with sharp declines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also estimate the potential errors of superflare frequency. The number of superflares found in the study is followed by the Poisson statistics, and the estimated σ is $\sqrt{{N}_{\mathrm{flare}}+1}$, where N_flare is the number of flares in each P_rot bin. We include these errors in Figures 9(a) and 10(a). Additionally, we include the following errors in Figures 9(b) and 10(b). The error from the sensitivity correction is from the number of stars that satisfies FAmp_assumed ≥ Threshold (${N}_{P}^{{{F}{\rm{Amp}}}_{\mathrm{assumed}}\geqslant {Threshold}}$). We estimate that σ is $\sqrt{{N}_{P}^{{{F}{\rm{Amp}}}_{\mathrm{assumed}}\geqslant {Threshold}}}$. The gyrochronology relation error is from the error values of Equations (12)–(14) of Mamajek & Hillenbrand (2008).

Comparing the results without and with corrections in Figures 9 and 10 shows two clear changes. First, because of the sensitivity correction, the smaller-amplitude events—which were originally missed, as shown in Figure 7—are now included in the calculation of the flare frequency, and the frequency value in the smaller-energy superflares (≲10³⁴ erg) becomes larger. The difference is larger, especially for rapidly rotating stars (P_rot < 5 days), since the detection completeness of rapidly rotating stars is smaller than that of slowly rotating stars, as shown in Figure 7. For example, the frequency value of flares with 10^33.75–10^34.0 erg on stars with T_eff = 5600–6000 K and P_rot < 5 days in Figure 10(a) (without correction) is <10⁻³⁵ erg⁻¹ yr⁻¹ star⁻¹, but that in Figure 10(b) (with correction) is >10⁻³⁵ erg⁻¹ yr⁻¹ star⁻¹. Second, because of the gyrochronology correction (using N_star instead of N_P in Table 3 as a value of the number of observed stars in each P_rot bin), the exact values of the flare frequencies change, while the shape of the overall distributions (e.g., the power-law distributions in the following) is not affected by this correction. For example, the frequency of slowly rotating Sun-like stars (T_eff = 5600–6000 K and P_rot = 20–40 days) becomes about half because of the gyrochronology correction: N_P(P_rot = 20–40 days)/N_P,all ∼ 32% and N_star,all(P_rot = 20–40 days)/N_star,all ∼ 70% in Table 3.

As our previous studies (e.g., Shibayama et al. 2013; Notsu et al. 2019) presented, these distributions can be described by a power law (dN/dE ∝ E^α), and rapidly rotating stars tend to have larger frequency values. The power-law indexes (α) of the distributions in Figures 9(b) and 10(b) are listed in Table 4. We estimated the power-law relation (dN/dE ∝ E^α) using all points from the P_rot bins shown in the figures, and the error is shown as the vertical error bars.

The indexes are roughly around α ∼ −2, which can be consistent with previous studies of superflares on solar-type stars (e.g., Shibata et al. 2013; Shibayama et al. 2013; Maehara et al. 2015; Notsu et al. 2019; Tu et al. 2020). We should note that the sensitivity correction in this study (see Figure 7) is a very simple one, only discussing the flare amplitude, and can be an incomplete correction, as we described above. Because of this, the power-law distributions in the smaller energy region (≲10³⁴ erg) can have some systematic errors that are not included in the error bars in Figures 9(b) and 10(b). We should keep this in mind when we see the power-law indexes listed in Table 4.

The upper-limit values of flare energy roughly depend on rotation period, as already seen in Figure 3. As also seen in the index values listed in Table 4, we can see that the power-law distribution becomes a bit steeper as the flare energy becomes larger. For example, the power-law index of superflares on Sun-like stars (T_eff = 5600–6000 K and P_rot = 20–40 days) in Figure 10(b) is α ∼ −2.7 if we calculate using the data of 10^33.75 erg < E_flare < 10^34.75 erg. However, if we use the data of 10^34.0 erg < E_flare < 10^34.75, the power-law index α is ∼−3.0. Similar tendencies can also be seen for the stars with P_rot < 20 days in Table 4, while we should note that the number of larger-energy flares is small, and the exact values of α should be treated with great caution. These tendencies might be related to the maximum flare energy cutoff of solar-type stars, since there is a correlation between the maximum spot coverage and the rotation period (Figure 6), and the spot size can explain the flare energy (Figure 4). In the case of old (age ∼ 4.6 Gyr), slowly rotating Sun-like stars with T_eff = 5600–6000 K and P_rot ∼ 25 days, as also described in Section 3.4, Figure 6(b) shows that the maximum size of the starspots is a few percent of the solar hemisphere. This corresponds to 10³⁴–10³⁵ erg on the basis of Equation (4), and the upper limit of the superflare energy of these Sun-like stars in Figures 3 and 10(b) is roughly in the same range.

Next, for reference, we see the flare frequency distributions of rapidly rotating stars with P_rot < 5 days in a bit more detail. Like Figures 9 and 10, Figures 11(a) and 12(a) are the flare frequency distributions of P_rot < 1, 1–3, and 3–5 days without the gyrochronology and sensitivity corrections, and Figures 11(b) and 12(b) are those with only the sensitivity correction. Since the gyrochronology relation has a larger scatter at a younger age (≲0.5–0.6 Gyr), as described in Section 3.2, the gyrochronology correction is not included in Figures 11(b) and 12(b). Although the original sample size of rapidly rotating stars is small (see N_P values in Table 3), we can still see the power-law distributions for very rapidly rotating stars in the range of E_flare ≳ 10^34.5 erg (e.g., P_rot < 1 day).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We then see the relation between flare frequency and rotation period by using the superflare data that include the above gyrochronology and sensitivity corrections. Figure 13 shows that the average flare frequency in a given period bin tends to decrease as the period increases in the range of P_rot longer than a few days, as we have already seen the same tendency in Figures 9(b) and 10(b). This result is basically the same as that presented in our previous study (Notsu et al. 2019), but the sample size is much larger, especially for slowly rotating stars (see Section 3.1). The frequency of superflares on young, rapidly rotating stars (P_rot = 1–3 days) is ∼100 times higher compared with old, slowly rotating stars (P_rot > 20 days), and this indicates that as a star evolves (and its rotational period increases), the frequency of superflares decreases. In contrast, the flare frequency is roughly constant in the range of P_rot ≲ 3 days. We note that in this range of P_rot ≲ 3 days, flare frequency distributions still show similar power-law distributions (Figures 11(b) and 12(b)). We can now interpret that this correlation between the rotation period (roughly corresponding to age) and flare frequency is consistent with the correlation between the rotation period (age) and the quiescent stellar activity level, such as the average X-ray luminosity (e.g., Noyes et al. 1984; Güdel 2007; Wright et al. 2011). The correlation between the rotation period and the X-ray luminosity of solar-type stars also shows a “saturation” regime in the range of P_rot ≲ 3 days and a “decreasing trend” in the range of P_rot ≳ 3 days.

Zoom In Zoom Out Reset image size

Our previous study (Shibayama et al. 2013; Maehara et al. 2015: Notsu et al. 2019) pointed out that the frequency distribution of superflares on Sun-like stars is roughly on the same power-law line. However, the definition of Sun-like stars in Shibayama et al. (2013) and Maehara et al. (2015; T_eff = 5600–6000 K and P_rot > 10 days) is different from that of this study (T_eff = 5600–6000 K and P_rot = 20–40 days). This indicates that many younger stars are included (e.g., stars with P_rot ∼ 10 days have an age of t ∼ 1 Gyr). Notsu et al. (2019) used the same definition of Sun-like stars as this study in order to discuss only the data of stars rotating as slowly as the Sun (P_rot ∼ 25 days and t ∼ 4.6 Gyr). However, as also described in this paper, the sample size of Notsu et al. (2019) was so small that we were not able to discuss the detailed properties accurately. Moreover, different from this study, our previous studies (Shibayama et al. 2013; Maehara et al. 2015; Notsu et al. 2019) did not include the effect of gyrochronology and flare detection completeness (sensitivity) correction when discussing the flare frequency distribution. Then, in Figure 14, we newly plot the frequency value of superflares on Sun-like stars with P_rot = 20–40 days (t > 3.2 Gyr) taken from Figure 10(b), in addition to the data of solar flares and superflares shown in Figure 17 of Notsu et al. (2019).

Zoom In Zoom Out Reset image size

The newly added data points of the superflares on Sun-like stars are roughly on the same power-law line. The exact value of the superflare frequency with P_rot = 20–40 days in this study is a bit smaller than that of Notsu et al. (2019), as clearly seen in Figure 14(b), mainly because the gyrochronology correction is newly included in this study. As also mentioned in Section 3.7, it might be possible that the frequency distribution of superflares becomes steeper around the upper limit of the superflare energy (E_flare ∼ 5 × 10³⁴ erg), but we cannot judge whether this possibility is true in the case of superflares on Sun-like stars because of the limited number of superflares detected on Sun-like stars in this study (see Table 2). We also note that some of the Sun-like superflare stars plotted in Figure 14 have flag 1 in Table 2. This means it is difficult to finally judge whether these stars really have rotation periods as slow as the Sun (P_rot = 20–40 days), and we need to be cautious (see Appendix A) when discussing the frequency of superflares on Sun-like stars in Figure 14. More investigations including spectroscopic measurements of rotation periods for all of these slowly rotating Sun-like superflare stars (e.g., Notsu et al. 2015b, 2019) are needed in the future.

As a result of Figure 14, we can roughly remark that superflares with energy E_flare ∼ 10^33.75, ∼10^34.0, ∼10^34.25, and ∼10^34.5 erg would be approximately once every ∼3000, ∼6000, ∼16,000, and ∼85,000 yr on old, slowly rotating Sun-like stars with P_rot ∼ 25 days and t ∼ 4.6 Gyr. These results of the superflare frequency distribution of Sun-like superflare stars (e.g., Figure 14) have shown the possibility that superflares whose energy is ∼7 × 10³³ erg (∼X700-class flares; see Figure 4) can occur on our Sun once every ∼3000 yr, those with ∼1 × 10³⁴ erg (∼X1000-class flares; see Figure 4) once every ∼6000 yr, and those larger than 1 × 10³⁴ erg can occur with longer timescales (e.g., once every ≳ 10,000 yr). This possibility of superflares on our Sun is now supported more strongly in this study because of a much larger number of superflares on slowly rotating Sun-like stars compared with our previous studies.

In contrast, the possibility of superflares with energy >10³⁴ erg has been questioned by several studies. Aulanier et al. (2013) suggested that it is unlikely that the current solar convective dynamo can produce the giant sunspot groups necessary for such large superflares (see Figure 4), since such giant sunspot groups have not been recorded in solar observations for a few hundred years (see Schrijver et al. 2012). Schmieder (2018) also supported the idea that with our Sun as it is today, it seems impossible to get larger sunspots and superflares with energy >10³⁴ erg. However, in Section 3.4, we have already seen that slowly rotating Sun-like stars (T_eff = 5600–6000 K, P_rot ∼ 25 days, and t ∼ 4.6 Gyr) have starspots with a size of ∼1% of the solar hemisphere (Figure 6). This spot size value is enough for ≳10³⁴ erg superflares on the basis of Equation (4), and the upper limit of the superflare energy of these Sun-like stars in this study (e.g., Figures 3 and 14) is roughly in the same range. The possible existence of such large starspots on Sun-like stars and the Sun is also suggested in the Kepler data analyses (Maehara et al. 2017; Notsu et al. 2019). Shibata et al. (2013) also suggested that even the current Sun can generate the large magnetic flux necessary for 10³⁴ erg superflares in a typical dynamo model.

Moreover, Maehara et al. (2017) and Notsu et al. (2019) also investigated the size–frequency distribution of these large starspot groups on Sun-like stars and that of sunspots and revealed that both sunspots and larger starspots could be related to the same physical processes. In Figure 18 of Notsu et al. (2019), they found that the occurrence frequency decreases as the area of sunspots or starspots increases, and the distribution of sunspot groups for large sunspots is roughly on this power-law line, although the appearance frequency of sunspots with spot areas of ∼1% of the solar hemisphere is about 10 times lower than that of starspots on Sun-like stars. The difference between the Sun and Sun-like stars might be caused by the lack of a “superactive” phase on our Sun during the past 140 yr (Schrijver et al. 2012). The upper limit of the starspot size values of Sun-like stars in Figure 18 of Notsu et al. (2019) is approximately a few percent of the solar hemisphere, and the appearance frequency of these spots is approximately once every 2000–3000 yr. As a result, these results of starspot frequency in Maehara et al. (2017) and Notsu et al. (2019) can be consistent with those of superflare frequency in this section (e.g., Figure 14).

Finally, it is interesting to note that several potential candidates of extreme solar flare events, which can be bigger than the largest solar flare in the past 200 yr (E_flare ∼ 10³² erg), have been reported from the data over the last 1000–2000 yr (e.g., see Usoskin 2017 and Miyake et al. 2019 for reviews). For example, significant radioisotope ¹⁴C and ¹⁰Be enhancements have been reported in tree rings and ice cores for the years AD 775, AD 994, and BC 660. They suggest extremely strong and rapid cosmic-ray increase events possibly caused by extreme solar flares (e.g., Miyake et al. 2012, 2013; Mekhaldi et al. 2015; O’Hare et al. 2019). These potential high-energy particle events can be caused by the superflares whose energies are at least larger than 10³³ erg (Cliver et al. 2020), although there is a large scatter between the flare energy and flux of energetic particles on the surface of the Earth (e.g., Takahashi et al. 2016), and more investigations are necessary. Besides, various potential low-latitude auroral and large sunspot drawings have also been reported from historical documents around the world, and they suggest the possibility that extreme solar flare events have occurred several times in the past ∼1000 yr (e.g., Hayakawa et al. 2017b & 2017b). In future studies, the superflare frequency information discussed in this study (e.g., Figure 14) should be compared with these radioisotope and historical data in detail, so that we can discuss the possibility of extreme solar flare events from various points of view.

Rubenstein & Schaefer (2000) argued that stars having a hot Jupiter-like companion are good candidates for superflare stars, since the hot Jupiter may play the role of a companion star in binary stars, such as RS CVn stars. The RS CVn stars are magnetically very active and produce many superflares. Since the Sun does not have hot Jupiters, Schaefer et al. (2000) concluded that our Sun has never generated superflares. In order to check the argument, we checked whether the superflare solar-type stars found in this study have any exoplanets, using the data retrieved from the NASA Exoplanet Archive¹² on 2020 May 28. As a result, the solar-type superflare stars found in this study have no confirmed exoplanets, while three candidate exoplanets¹³ and three false-positive exoplanets¹⁴ are found. This suggests that a hot Jupiter is not a necessary condition for superflares on solar-type stars. In addition, Shibata et al. (2013) also found from theoretical considerations that hot Jupiters do not play an essential role in the generation of magnetic flux in the star itself, if we consider only the magnetic interaction between the star and the hot Jupiter, whereas the tidal interaction seems to be a possible cause of enhanced dynamo activity, though more detailed studies will be necessary. This seems to be consistent with the above observational finding that a hot Jupiter is not a necessary condition for superflares on solar-type stars.