Do Kepler Superflare Stars Really Include Slowly Rotating Sun-like Stars?—Results Using APO 3.5 m Telescope Spectroscopic Observations and Gaia-DR2 Data

In Section 3 and Appendix A.1, we summarize the results of our new APO 3.5 m spectroscopic observations of 18 superflare stars.

We reported that more than half (13 stars) among these 18 stars have no obvious evidence of being in a binary system. Combined with the results of the 50 stars that we observed with Subaru telescope in Notsu et al. (2015a), we have conducted spectroscopic observations of 64 superflare stars in total (Table 6). Four stars (KIC 4742436, KIC 4831454, KIC 9652680, and KIC 11610797) among these 64 stars were observed in both studies. As a result, 43 stars among the total 64 solar-type superflare stars are classified as "single" stars. However, we need to remember here that we cannot completely exclude the possibility that some of these 43 "single" superflare stars have companions since observations and analyses in this study are limited, as we have also described in detail in Section 4.1 of Notsu et al. (2015a).

For example, as for the target stars with "multiple" observations in this study, only those showing large RV changes would likely be detected in the randomly spaced observations. Targets with longer-period orbits would require more observations spaced accordingly. Thus, only the short-period systems likely would be captured, and even then some would be missed by accident owing to poorly spaced observations (e.g., the case that the time differences between the observations correspond to the orbital period of the binary system). We must note these points, but we consider that more detailed analyses of binarity are not really necessary for the overall discussion of stellar properties of superflare stars in this paper.

Two target stars, KIC 11551430 and KIC 7093428, are classified as visual binary stars on the basis of the slit viewer images, Figures 19 and 20 in Appendix A.1. Pixel count data of these two stars suggested that superflares occur on the primary star KIC 11551430A as for the visual binary system KIC 11551430, while flares occur on the fainter companion stars KIC 7093428B or KIC 7093428C as for the system KIC 7093428. Measurement results of rotation velocity v sin i and chromospheric line intensities in this study (listed in Tables 4 and 5) support this suggestion. As for KIC 11551430, the primary star KIC 11551430A rotates much more rapidly (v sin i ∼ 18.8 km s⁻¹) than the companion star KIC 11551430B (v sin i < 4 km s⁻¹). Moreover, the primary KIC 11551430A shows strong chromospheric emissions, while the companion KIC 11551430B does not show any strong emissions (Figures 26–28). The primary KIC 11551430A can have much higher probability of generating superflares since superflare stars are generally well characterized with rapid rotation velocity and high chromospheric activity levels. In contrast, the primary KIC 7093428A has no strong chromospheric emissions (Figures 26–28), and this primary star does not show any properties as superflare stars. It is therefore highly possible that flares occur on the fainter companion stars KIC 7093428B or KIC 7093428C mentioned above.

As seen in Table 1, stars identified as binary stars tend to show a larger number of flares. For example, KIC 11551430 and KIC 4543412, which show RV shifts, have N_flare = 89 and 13, respectively. All the stars identified as "single" stars have N_flare < 10 in Table 1. As a result, among the 187 superflares (on 23 stars) from Kepler 1-minute cadence data in Maehara et al. (2015), at least 127 events are found to occur on binary stars. This means that a large portion of the superflare data of the Kepler 1-minute cadence sample (Maehara et al. 2015) are contaminated by binary stars (e.g., close binary stars). The data from Maehara et al. (2015) are therefore not enough for investigating the possibility of superflares on Sun-like stars, and we need more investigations by increasing the number of single superflare stars (e.g., See also the first paragraph of Section 6.3).

We estimated the atmospheric parameters (T_eff, log g, and [Fe/H]) of the target superflare stars in Appendix A.2. The measured values of the 13 single target stars of this observation are T_eff = 5400–6300 K, log g = 3.9–4.8, and [Fe/H] = (−0.3) – (+0.3), respectively (Figure 1). This indicates that the stellar parameters of these 13 single target superflare stars are roughly in the range of solar-type (G-type main-sequence) stars, though the stars with log g ≲ 4.0 are possibly subgiant G-type stars (for discussions on evolutionary state of the target stars, see also Section 5.1). Compared with the 34 single target superflare stars that we observed in Notsu et al. (2015a), the target stars tend to have a bit hotter T_eff values, and most of them are "solar-analog" stars (early-type G-type main-sequence stars). No clear "metal-rich" or "metal-poor" stars are included in the target stars of this observation.

Zoom In Zoom Out Reset image size

As described above, we report the values of projected rotational velocity v sin i, stellar radius R (R_Gaia and R_spec), and rotation period P_rot from the brightness variation of the single target superflare stars (listed in Table 4). Using the R and P_rot values, we can estimate the rotational velocity (v_lc):

$$\begin{eqnarray}&&{v}_{\mathrm{lc}}=\displaystyle \frac{2\pi R}{{P}_{\mathrm{rot}}},\end{eqnarray} \tag{ 1 }$$

as also described in Section 4 of Notsu et al. (2015b). As for the R values, we used R = R_Gaia as a first priority and R = R_spec only for the stars without R_Gaia values. As for the errors of v_lc, we took the root sum squares of the two types of errors from R (see the descriptions in Appendix A.3) and P_rot. We here assume that the possible errors of P_rot are about 10%, as done in Notsu et al. (2015b). In this assumption, we very roughly consider typical possible differences between equatorial rotation period and the measured period caused by the solar-like differential rotations. For example, the solar latitudinal differential rotation has a magnitude of 11% of the average rate from equator to midlatitudes (≈45° latitude; see Benomar et al. 2018). We must note here that the differences can be much larger if we observe starspots near the pole region, and there can be solar-type stars with much larger differential rotation magnitudes (Benomar et al. 2018). Moreover, it is sometimes hard to distinguish the correct P_rot values with 1/n × P_rot (n = 2, 3, ...) values in the periodogram, especially when there are several starspots on the surface of the stars (see Figure 2 of Notsu et al. 2013b). As a result, we must note here with caution that P_rot can have much larger error values, though the detailed investigations of each error value are beyond the scope of the overall discussions of v sin i versus v_lc in Figure 2. The resultant values of v_lc are listed in Table 4.

Zoom In Zoom Out Reset image size

We then plot v sin i as a function of v_lc in Figure 2. Not only the target superflare stars of this study (13 "single" stars and two binary stars categorized as "RV" and "VB" in Appendix A.1) but also the 34 single superflare stars that we have observed in Notsu et al. (2015b) are plotted for reference. The v_lc values of the latter 34 superflare stars are updated from the original ones in Notsu et al. (2015b). We newly recalculated v_lc (see Equation (1)) of each target star by using the R_Gaia value if the star has an R_Gaia value reported in Berger et al. (2018).

In Figure 2, v sin i tends to be smaller than v_lc, and such differences should be explained by the inclination angle effect, as mentioned in the previous studies (e.g., Hirano et al. 2012; Notsu et al. 2015b). On the basis of v sin i and v_lc, the stellar inclination angle (i) can be estimated by using the following relation:

$$\begin{eqnarray}&&i=\arcsin \left(\displaystyle \frac{v\,\sin \,i}{{v}_{\mathrm{lc}}}\right).\end{eqnarray} \tag{ 2 }$$

In Figure 2, we also show four lines indicating i = 90°(v sin i = v_lc), i = 60°, i = 30°, and i = 10°. First of all, almost all the stars (except for KIC 6291837, with v sin i ∼ 7.9 km s⁻¹ and v_lc ∼ 5.0 km s⁻¹) in Figure 2, the relation "v sin i ≲ v_lc" is satisfied. This is consistent with the assumption that the brightness variation is caused by the rotation since the inclination angle effect can cause this relation "v sin i ≲ v_lc" if v_lc really corresponds to the rotation velocity (i.e., v = v_lc). This is also supported by the fact that the distribution of the data points in Figure 2 is not random, but is concentrated between the lines of i = 90° and i = 60°. The distribution is expected to be much more random if the brightness variations have no relations with the stellar rotation. In addition, stars that are distributed in the lower right side of Figure 2 are expected to have small inclination angles and to be nearly pole-on stars. Later in Section 5.2, we see these inclination angle effects from another point of view with the scatter plot of flare energy and starspot size.

Summarizing the results in this section, we can remark that rotation velocity values from spectroscopic results (v sin i) and those from Kepler brightness variation (v_lc) are consistent, and this supports that the brightness variation of superflare stars is caused by the rotation. This remark was already suggested in Notsu et al. (2015b), but the conclusions are more strongly confirmed by the new spectroscopic observations and Gaia-DR2 stellar radius (R_Gaia) values.

We measured the core intensity and flux values of chromospheric lines (Ca ii λ8542, Hα λ6563, and Ca ii H and K lines) in Appendices A.5–A.7. These lines have been widely used for investigating stellar chromospheric activities (e.g., Wilson 1978; Linsky et al. 1979a; Noyes et al. 1984; Duncan et al. 1991; Baliunas et al. 1995; Hall 2008) and are good indicators of stellar average magnetic fields (e.g., Schrijver et al. 1989; Notsu et al. 2015b). Only the Ca ii λ8542 line is mainly used in our previous studies using the Subaru telescope (Notsu et al. 2013a, 2015b; Nogami et al. 2014). In this study, Ca ii H and K lines are used with high-dispersion spectra for the first time in our series of studies of the Kepler solar-type superflare stars.¹⁸

In Figure 3, we compare the Ca ii λ8542 and Ca ii H and K index values (r₀(8542) and $\mathrm{log}{{ \mathcal R }}_{\mathrm{HK}}^{+}$) with the amplitude of the brightness variation of Kepler data ($\langle \mathrm{BVAmp}\rangle$). The r₀(8542) index is the normalized intensity at the center of the Ca ii λ8542 line (see Section A.5 for details). The $\mathrm{log}{R}_{\mathrm{HK}}^{+}$ index is a universal and "pure" Ca ii H and K activity indicator introduced by Mittag et al. (2013) and is defined as ${{ \mathcal R }}_{\mathrm{HK}}^{+}={{ \mathcal F }}_{\mathrm{HK}}^{+}/\sigma {T}_{\mathrm{eff}}^{4}$ (Equation (18)), where ${{ \mathcal F }}_{\mathrm{HK}}^{+}$ is the Ca ii H and K surface flux (unit: erg cm⁻² s⁻¹) with photospheric and "basal" flux contributions removed (see Appendix A.7 for details). As described in Figure 30 in Appendix A.7, this $\mathrm{log}{R}_{\mathrm{HK}}^{+}$ index has a good correlation with the r₀(8542) index, and $\mathrm{log}{R}_{\mathrm{HK}}^{+}$ can be more sensitive to the difference in the lower activity level region compared with r₀(8542). $\langle \mathrm{BVAmp}\rangle$ values are calculated by taking the average of the amplitude value of each Quarter (Q2–Q16) data, as we have done in Notsu et al. (2015b). The resultant values of $\langle \mathrm{BVAmp}\rangle$ of the target stars of our APO 3.5 m observations are listed in Table 5. In Figure 3, we also plot solar values as a reference, as done in Notsu et al. (2015b), and we can see that, as for most of the target superflare stars, Ca ii λ8542 and Ca ii H and K index values are higher than the solar values. This suggests that these superflare stars have higher chromospheric activities compared with the Sun.

Zoom In Zoom Out Reset image size

More importantly, there is a rough positive correlation between Ca ii index values and $\langle \mathrm{BVAmp}\rangle$ in Figure 3. Assuming that the brightness variation of superflare stars is caused by the rotation of a star with starspots, the brightness variation amplitude ($\langle \mathrm{BVAmp}\rangle$) corresponds to the starspot coverage of these stars. As mentioned above, Ca ii index values are good indicators of stellar average magnetic field (or total magnetic flux). Then, we can remark that there is a rough positive correlation between the starspot coverage from Kepler photometric data and the stellar average magnetic field from spectroscopic data. All the target stars expected to have large starspots on the basis of their large amplitude of the brightness variation show strong average magnetic field compared with the Sun. In other words, our assumption that the amplitude of the brightness variation corresponds to the spot coverage is supported, since the average magnetic field is considered to be caused by the existence of large starspots. These results have already been confirmed in Notsu et al. (2015b) with the Ca ii λ8542 line, but in this study we confirmed the same conclusion with the following two updates: (i) larger number of stars with Ca ii λ8542 line data (Figure 3(a)), and (ii) analyses using Ca ii H and K lines (Figure 3(b)).

We also estimated Li abundances A(Li) of the target superflare stars in Appendix A.8. Li abundance is known to be a clue for investigating the stellar age of solar-type stars (e.g., Skumanich 1972; Takeda et al. 2010, 2013; Honda et al. 2015; Notsu et al. 2017). The Li depletion is seen in the stars with T_eff ≲ 5500 K, while the stars with T_eff ≳ 6000 K show no Li depletion. This is because as the star becomes cooler, the convection zone in the stellar atmosphere evolves and the Li is transported to a deeper hotter zone, where Li is easily destroyed (p, α reactions: ⁷Li, T ≥ 2.5 × 10⁶ K; ⁶Li, T ≥ 2.0 × 10⁶ K). The depletion of Li in the stellar surface caused by convective mixing increases as time passes, and we can remark that young stars tend to have high A(Li) values.

In Figure 4, we plot A(Li) as a function of T_eff of the superflare stars in comparison with the ordinary stars in the Hyades (taken from Takeda et al. 2013). Some of the superflare stars show high Li compared with the stars in the Hyades (Figure 4), and such stars are suggested to be younger than the Hyades. The age of the Hyades is estimated to be 6.25 × 10⁸ yr (e.g., Perryman et al. 1998). It is reasonable that such young stars tend to have high activity levels and produce superflares. However, more than 10 target stars do not show higher A(Li) values compared with the Hyades, and a few of them are as low as the solar value. These results suggest that superflare stars include many young stars but also include old stars like our Sun, as also suggested in Honda et al. (2015).

Zoom In Zoom Out Reset image size

In this section, we investigate again the statistical discussions of Kepler solar-type superflare stars in our previous studies (Maehara et al. 2012, 2015, 2017; Notsu et al. 2013b; Shibayama et al. 2013), by incorporating Gaia-DR2 stellar radius estimates (reported in Berger et al. 2018) and the results of our spectroscopic observations (Subaru and APO 3.5 m observations). In these previous studies, we reported 1547 superflare events on 279 solar-type (G-type main-sequence) stars from Kepler 30-minute (long) time-cadence data (Shibayama et al. 2013) and 189 superflares on 23 solar-type stars from Kepler 1-minute (short) time-cadence data (Maehara et al. 2015; line 1 of Tables 7 and 8). As also used in Appendix A.3, Berger et al. (2018) reported the catalog of the 177,911 Kepler stars with the stellar radius estimates (R_Gaia) estimated from Gaia-DR2 parallax values. A total of 245 out of the 279 solar-type superflare stars in Shibayama et al. (2013) and 18 of the 23 stars in Maehara et al. (2015) have R_Gaia values in this Gaia-DR2 catalog, respectively (Line 2 of Tables 7 and 8).

Table 7.  Number of Solar-type Superflare Stars (${N}_{\mathrm{star}}$) and Superflares on These Stars (${N}_{\mathrm{flare}}$) Found from Kepler 30-minute Time-cadence Data (See Figures 5(a) and 6)

	${N}_{\mathrm{star}}$	${N}_{\mathrm{flare}}$
(1) Original data (Shibayama et al. 2013)	279	1547
(2) Stars having RGaia values (Berger et al. 2018) among (1)	245	1402
(3) Stars identified as main sequence among (2) (see Figure 5(a))	136 (55.5%)	496
(4) Stars identified as subgiants among (2) (see Figure 5(a))	108 (44.1%)	905
(5) Stars identified as red giants among (2) (see Figure 5(a))	1 (0.4%)	1
(6) Stars with Teff,DR25 = 5100–6000 K among (3)	106	419
(7) Stars that were originally early K dwarfs in Candelaresi et al. (2014) but are newly identified as solar-type stars with Teff,DR25 and RGaia	36	178
(8) (6)+(7)	142	597
(9) Stars having ΔF/F values in McQuillan et al. (2014) among (8)	113	527

Download table as:  ASCIITypeset image

We plot these superflare stars on the stellar radius (R) versus T_eff diagram in Figures 5(a) and (b). In these figures, all the Kepler stars reported in Berger et al. (2018) are also plotted for reference with the evolutionary state classifications (main sequence (MS)/subgiants/red giants/cool main-sequence binaries). In Figure 5(a), among the 245 stars found as solar-type superflare stars from Kepler 30-minute time-cadence data in Shibayama et al. (2013), 136 stars (55.5%) are classified as main-sequence stars, while 108 stars (44.1%) are classified as subgiants and only one star (KIC 4633721) as a red giant (lines (3)–(5) of Table 7). Originally in Shibayama et al. (2013), we used the T_eff and log g values in the initial Kepler Input Catalog (KIC; Brown et al. 2011) and selected solar-type (G-type main-sequence) stars with the definition of 5100 K ≤ T_eff ≤ 6000 K and log g ≥ 4.0. T_eff and log g values in the initial KIC have large error values of ±200 K and 0.4 dex, respectively (Brown et al. 2011), and the reliability of each value can be low. In Figure 5(a), in addition to the updated stellar radius R_Gaia from Gaia-DR2 data, T_eff values are also updated with the latest DR25-KSPC (Mathur et al. 2017), which incorporates the revised method of T_eff estimation (see Pinsonneault et al. 2012). As a result, a large fraction (more than 40%) of the stars that were originally identified as solar-type superflare stars in Shibayama et al. (2013) are now classified as subgiant or red giant stars.

Zoom In Zoom Out Reset image size

In Figure 5(b), among the 18 stars found as solar-type superflare stars from Kepler 1-minute time-cadence data in Maehara et al. (2015), 13 stars are classified as main-sequence stars, while 5 stars are classified as subgiant stars (lines (3)–(5) of Table 8). In the case of these 18 stars, the revised KIC catalog (Huber et al. 2014) was used when we selected solar-type stars in Maehara et al. (2015), and it is possible that the errors of T_eff were somewhat smaller, compared with the case of the stars from Shibayama et al. (2013) in Figure 5(a).

In Figure 5(c), we plot 50 superflare stars found from Kepler 30-minute cadence data that we have spectroscopically observed using the Subaru telescope (Notsu et al. 2015a, 2015b). As for the radius value in the vertical axis of these two figures, the R_Gaia value is used if it exists. If the star has no R_Gaia value in Berger et al. (2018), R_spec values (in Table 4 of Notsu et al. 2015a) are used for single stars and R_DR25 values (from DR25-KSPC in the above) are used for binary stars. As for the temperature value in the horizontal axis, T_eff values estimated from spectroscopic data (in Table 4 of Notsu et al. 2015a) are used for single stars, while T_eff,DR25 values are used for binary stars. As a result, 24 stars are classified as main sequence, while the remaining 10 stars¹⁹ are classified as subgiants, among all 34 "single" target stars. Among the 16 "binary" target stars, only 5 stars are classified as main-sequence stars, while the remaining 11 stars²⁰ are classified as subgiants.

In Figure 5(d), we plot the target stars of the APO 3.5 m spectroscopic observation of this study, which are equal to all 23 superflare stars found from Kepler 1-minute cadence data (Maehara et al. 2015). The radius and temperature values in the vertical and horizontal axes are plotted in the basically same way as in Figure 5(c), but R_spec values are used for "RV" binary stars if the stars have no R_Gaia values, and T_eff estimated from spectroscopic data are used for all "RV" binary stars. As a result, 9 stars are classified as main sequence, while 4 stars²¹ are classified as subgiants, among all 13 "single target stars. Among the five "binary" target stars, three stars are classified as main sequence, while two stars²² are classified as subgiants. Moreover, among the remaining five stars for which we did not get spectroscopic data with enough S/N in this study (see Section 2),²³ three stars are classified as main sequence, while two stars²⁴ are classified as subgiants. In the following sections, we update statistical studies by using the data of the stars that we can newly classify as solar-type (G-type main-sequence) stars.

Most superflare stars show large-amplitude brightness variations, and they suggest that the surfaces of superflare stars are covered by large starspots (see Section 4.4). Figure 6 shows the scatter plot of flare energy (E_flare) as a function of the spot group area (A_spot). The values of solar flares in Figure 6 are the same as those in our previous studies (Notsu et al. 2013b; Shibata et al. 2013; Maehara et al. 2015). We estimate bolometric energies (white-light flare [WLF] energies) of solar flares from GOES soft X-ray (SXR) flux values. We use the relation that WLF energy (E_WL) is proportional to SXR flux (F_SXR), E_WL ∝ F_SXR, on the basis of the results and detailed descriptions in Section 4.1 of Namekata et al. (2017). This relation is supported by the observational comparisons between WL and SXR data during flares (e.g., Figure 4 of Namekata et al. 2017), which are related to the well-known relation between hard X-ray flux and SXR flux during a flare ("Neupert effect"; Neupert 1968). As a result, we here assume that bolometric energies of B-, C-, M-, X-, and X10-class solar flares are 10²⁸, 10²⁹, 10³⁰, 10³¹, and 10³² erg, respectively.

Zoom In Zoom Out Reset image size

In Figure 6, A_spot of superflare stars were estimated from the normalized amplitude of light variations (ΔF/F) by using the following equation (Notsu et al. 2013b; Shibata et al. 2013):

$$\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 is the apparent area of the star and T_spot and T_star are the temperature values of the starspot and photosphere of the star. In Figure 6(a), original values of superflares on solar-type stars reported in Shibayama et al. (2013) (from Kepler 30-minute cadence data) and Maehara et al. (2015) (from Kepler 1-minute cadence data) are plotted²⁵ for comparisons with the results of this study. In these previous papers, we defined the amplitude (ΔF/F) as the brightness range normalized by the average stellar brightness, in which the lower 99% of the distribution of the brightness difference from average, except for the flares, is included. A_star and T_star values were based on the stellar radius and temperature values used in these previous papers,²⁶ and we assumed T_spot = 4000 K.

In this study, we newly updated both of the E_flare and A_spot values by using the latest T_eff,DR25 and R_Gaia values described in Section 5.1. First, from the superflare stars that we previously reported from Kepler 30-minute and 1-minute time-cadence data (Shibayama et al. 2013; Candelaresi et al. 2014; Maehara et al. 2015) we selected the stars classified again as solar-type stars (main-sequence stars with T_eff = 5100–6000 K) on the basis of T_eff,DR25 and R_Gaia values (see Figure 5; line (6) of Tables 7 and 8). Not only "subgiants" and "red giants" but also stars identified as "cool main-sequence binary" stars in Figure 5 are not included in the "solar-type" stars classified in this study.

As for the selection from Kepler 30-minute cadence data, we used the data not only from Shibayama et al. (2013) (described in Figure 5(a)) but also from Candelaresi et al. (2014). We investigated superflares on solar-type stars in Shibayama et al. (2013) and those on G-, K-, and M-type stars in Candelaresi et al. (2014). Exactly the same Kepler data set (the first 500 days: Quarters 0–6) is used for the flare surveys in both these studies, and the superflare data of Shibayama et al. (2013) (solar-type stars) exactly correspond to the subset of those of Candelaresi et al. (2014) (G-, K-, and M-type stars). Then in this paper, we newly selected solar-type superflare stars not only from the original solar-type superflare stars in Shibayama et al. (2013) (see Figure 5(a)) but also from the stars originally identified as K- and M-type stars in Candelaresi et al. (2014). A total of 36 superflare stars originally identified as early K-type stars (with T_eff = 4900–5100 K) in Candelaresi et al. (2014) are newly categorized as solar-type stars with the revised T_eff,DR25 and R_Gaia values (line (7) of Table 7). We note that these superflare stars were not included in Figure 5(a) but are included in the following discussions of solar-type superflare stars in this paper. As a result, 142 solar-type stars (with 597 superflares in total) from Kepler 30-minute cadence data (Shibayama et al. 2013; Candelaresi et al. 2014; line (8) of Table 7) and 10 solar-type stars (with 51 superflares in total) from 1-minute cadence data (Maehara et al. 2015) are selected (line (6) of Table 8).

We then recalculated E_flare and A_spot for these selected data. The recalculation of E_flare was done by applying these updated T_eff,DR25 and R_Gaia values to the equations presented in Section 2.3 of Shibayama et al. (2013). As for A_spot, we used the new methods presented in Maehara et al. (2017). ΔF/F values taken from McQuillan et al. (2014) were used, and we applied the following relation on the temperature difference between photosphere and spot (T_star − T_spot) deduced from Berdyugina (2005):

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}T({T}_{\mathrm{star}}) & = & {T}_{\mathrm{star}}-{T}_{\mathrm{spot}}\\ & = & \ 3.58\times {10}^{-5}{T}_{\mathrm{star}}^{2}+0.249{T}_{\mathrm{star}}-808.\end{array}\end{eqnarray} \tag{ 4 }$$

As a result, we got the resultant values of the updated E_flare and A_spot values for 113 stars (with 527 flares in total) from Kepler 30-minute cadence data (Shibayama et al. 2013; Candelaresi et al. 2014; line (9) of Table 7) and eight stars (with 48 flares in total) from 1-minute cadence data (Maehara et al. 2015; line (7) of Table 8). We note that only the stars having ΔF/F values in McQuillan et al. (2014) are selected here. These resultant superflare values are plotted in Figure 6(b).

In Figure 6, the majority of superflares occur on the stars with large starspots, though there is a large scatter. Flares are sudden releases of magnetic energy stored around starspots (see Shibata & Magara 2011). 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 paper (e.g., Shibata et al. 2013) suggested that the upper limit of E_flare can be determined by the simple scaling law:

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{flare}} & \approx & {{fE}}_{\mathrm{mag}}\approx \displaystyle \frac{{B}^{2}{L}^{3}}{8\pi }\approx \displaystyle \frac{{B}^{2}}{8\pi }{A}_{\mathrm{spot}}^{3/2}\\ & \approx & \ 7\times {10}^{32}(\mathrm{erg})\left(\displaystyle \frac{f}{0.1}\right){\left(\displaystyle \frac{B}{{10}^{3}{\rm{G}}}\right)}^{2}{\left(\displaystyle \frac{{A}_{\mathrm{spot}}}{3\times {10}^{19}{\mathrm{cm}}^{2}}\right)}^{3/2}\\ & \approx & \ 7\times {10}^{32}(\mathrm{erg})\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{ 5 }$$

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. The black solid and dashed lines in Figure 6 represent Equation (5), and almost all the solar flare data are located below these lines. As for the original superflare data from our previous papers (Figure 6(a)), many of the superflares are located below the solid line, but some of them are located above this line. Our previous papers (e.g., Notsu et al. 2013b, 2015b) considered that these stars above this line are expected to have low inclination angle or have starspots around the pole region. In contrast, as for the new data of superflares updated in this study (Figure 6(b)), the number of flares above this line (Equation (5) with B = 3000 G) looks to be decreased, and a much higher fraction of superflares are located below this line. In order to see this point more clearly, we plot only the maximum energy flare of each superflare star in Figures 7(a) and (b) and compare them using histograms in Figure 8. The number of stars above the black solid line (the upper-limit line from Equation (5) with B = 3000 G and the inclination angle i = 90°) is significantly decreased. With these updates, it is more strongly supported that the upper limit of the energy released by the flare is not inconsistent with the magnetic energy stored around the starspots.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we see this point (the difference of the number of superflare stars above Equation (5) between the previous studies and this study) a bit more in detail by incorporating the results of our spectroscopic studies in Section 4. In Figure 2 in Section 4.3, we compared the projected rotational velocity (v sin i) with the stellar rotational velocity (v_lc) on the basis of our spectroscopic observations (Notsu et al. 2015b; this study), and stellar inclination angle (i) can be calculated using Equation (2). We plot again the data of single superflare stars in Figure 9. In this figure, we newly classified solar-type stars by using spectroscopically measured temperature (T_eff,spec) and stellar radius values (R_gaia or R_spec) in Figures 5(c) and (d). White open symbols and black filled ones are solar-type stars (main-sequence stars with T_eff = 5100–6000 K). The former points are the stars having A_spot values deduced from ΔF/F in McQuillan et al. (2014) so that we can plot the data points in Figure 10(b), but the latter ones are the stars without ΔF/F values in McQuillan et al. (2014) so that we cannot plot the data in Figure 10(b). Orange and green symbols correspond to the main-sequence stars with a bit hotter (T_eff = 6000–6300 K) and a bit cooler (T_eff = 5000–5100 K) temperature values, respectively. Red symbols are the stars classified as subgiants (see Section 5.1). As a result, many of the low inclination angle stars are not classified as solar-type stars in this study.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Then, in Figures 10(a) and (b), we overplotted the maximum energy flare data of these spectroscopically observed superflare stars on the data of Figures 7(a) and (b), respectively. In Figure 10(a), the starspot size and flare energy data of the overplotted stars are taken from those used in Figure 7(a).²⁷ As for Figure 10(b), we update starspot size and flare energy values of the overplotted stars by recalculating them with T_eff,spec and stellar radius values (R_gaia or R_spec) used in Figures 5(c) and (d). The equations used for these recalculations are the same as those used for Figures 6(b) and 7 in the above. Notsu et al. (2015b) suggested that many of the data points located on the upper left side of the A_spot versus E_flare diagram tend to have low inclination angle values. This tendency can also be seen in Figure 10(a), but most of these low inclination angle stars are now classified as subgiant stars (i.e., red points in this figure). In the updated Figure 10(b), the data points are located more in the lower right side of this diagram. There is only one star above the black solid line (Equation (5)). However, the number of stars is decreased in Figure 10(b) compared with Figure 10(a), and we need to increase the number of target stars with more spectroscopic observations in the future.

Finally, in Figures 11(a) and (b), we plot again A_spot and E_flare values (taken from Figure 6(b)), separating the superflare data on the basis of the stellar temperature values. The data of the stars with a temperature range (T_eff = 5600–6000 K) close to the solar temperature are plotted in Figure 11(b), while those of late G-type main-sequence stars (T_eff = 5100–5600 K) are in Figure 11(a). There are no big differences between Figures 11(a) and (b). As a result, almost all the data points of superflares are located below the line of Equation (5), and we confirmed again that the upper limit of the energy released by flares is not inconsistent with the magnetic energy stored around the starspots.

Zoom In Zoom Out Reset image size

Previous observations of solar-type stars (e.g., X-ray, UV, Ca ii H and K) have shown that the stellar magnetic activity level decreases as rotation period increases (e.g., Noyes et al. 1984; Güdel 2007; Wright et al. 2011). Since stellar age has a good correlation with rotation, young, rapidly rotating stars show the higher activity levels, and it was expected that they show more energetic flares more frequently compared with slowly rotating stars like the Sun. In our previous paper (Notsu et al. 2013b), we investigated the relationship between the superflare energy (E_flare) and the rotation period (P_rot) and suggested that the maximum superflare energy in a given rotation period bin does not have a clear correlation with the rotation period (Figure 12(a) of this paper), while the average flare frequency in a given period bin has a correlation with the rotation (Figure 7(b) of Notsu et al. 2013b and Figure 2 of Maehara et al. 2015). This suggestion is important since, against the above-mentioned previous expectations, energetic superflares with ∼10³⁵ erg can occur on solar-type stars rotating as slow as the Sun (P_rot ∼ 25 days), even though the frequency is low (once in a few thousand years), compared with rapidly rotating stars.

Zoom In Zoom Out Reset image size

Next, we investigate again this relation by using superflare values updated with Gaia-DR2 stellar radius (R_Gaia) data. In Section 5.2, we newly classified solar-type stars on the basis of T_eff,DR25 and R_Gaia and recalculated flare energy E_flare. Figure 12(b) shows the relationship between this updated E_flare value and the rotation period P_rot of each solar-type superflare star. P_rot values plotted here are taken from McQuillan et al. (2014). Unlike the results of our previous studies (Figure 12(a)), Figure 12(b) suggests that the upper limit of E_flare in each period bin has a continuous decreasing trend with rotation period. For example, stars rotating as slow as the Sun (P_rot ∼ 25 days) show superflares up to 10³⁵ erg, while rapidly rotating stars with P_rot = 1–3 days have more energetic superflares, up to 10³⁶ erg.

Flare activity also depends on stellar temperature (Candelaresi et al. 2014; Davenport 2016; Van Doorsselaere et al. 2017), and even among solar-type stars, cooler stars can have higher flare activities (see Figure 5 of Shibayama et al. 2013). Then, the data of superflares plotted in Figure 13 are separated into panels (a) and (b) on the basis of the stellar temperature values. The data of the stars with a temperature range (T_eff = 5600–6000 K) close to the solar temperature are plotted in Figure 13(b), while those of late G-type main-sequence stars (T_eff = 5100–5600 K) are in Figure 13(a). In Figure 13(b), we also added the scale of stellar age (t) on the basis of the gyrochronology relation of solar-type stars (t ∝ ${P}_{\mathrm{rot}}^{0.6}$; Ayres 1997) in order to compare the age of superflare stars with that of the Sun (t = 4.6 Gyr). 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 solar age (t = 4.6 Gyr), a breakdown of gyrochronology relations has been recently reported (van Saders et al. 2016; Metcalfe & Egeland 2018). With this scale, for example, we can see that Sun-like stars with P_rot ∼ 25 days are more than four times older (t ∼ 4.6 Gyr) than the stars with P_rot ∼ 10 days (t ∼ 1 Gyr).

Zoom In Zoom Out Reset image size

As a result of Figures 13(a) and (b), we confirmed again the suggestions from Figure 12(b). The upper limit of superflare energy in each period bin depends on the rotation period in both Figures 13(a) and (b), and there is a one-order-of-magnitude difference between the maximum flare energy on young stars (t < 0.5 Gyr) and that on Sun-like old stars (t ∼ 4.6 Gyr). We cannot judge whether there are any clear differences between these two figures because of the low number statistics. It is necessary to increase the number of superflare data with future observations (see also Section 6.3).

Then, we see the relation between flare frequency and rotation period by using the updated superflare data. Figure 14 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. The frequency is averaged from all of the superflare stars in the same period bin. This result is basically the same as that presented in Notsu et al. (2013b) with the original superflare data. 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 ∼ 25 days), and this indicates that as a star evolves (and its rotational period increases), the frequency of superflares decreases. 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 previous measurement of the stellar activity level such as the average X-ray luminosity (e.g., Noyes et al. 1984; Güdel 2007; Wright et al. 2011). Summarizing the results in this section, superflares with their energy ≲5 × 10³⁴ erg can occur on old, slowly rotating Sun-like stars (T_eff = 5600–6000 K, P_rot ∼ 25 days, and t ∼ 4.6 Gyr), even though the frequency and maximum flare energy are lower compared with young, rapidly rotating stars and cooler stars.

Zoom In Zoom Out Reset image size

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. Here we update these values by using Gaia-DR2 stellar radius, as also done for superflare stars in the above. As done in Maehara et al. (2017), we used P_rot and brightness variation amplitude ΔF/F values reported in McQuillan et al. (2014). We newly classified 49,212 solar-type stars (N_data in Table 9 in Appendix B) among the stars in McQuillan et al. (2014) by using T_eff,DR25 and R_Gaia values, as done for superflare stars in Section 5.2. P_rot and ΔF/F values of 11,594 stars (see N_P(all) values in Table 9) are detected among these 49,212 stars. Then, we recalculated A_spot values from ΔF/F values, using these T_eff,DR25 and R_Gaia values. The method of A_spot calculation is exactly the same as those used for superflare stars in Section 5.2. The resultant values are plotted in Figure 15. The data in Figure 15 are separated into panels (a) and (b) on the basis of the temperature values. The data of the stars with a temperature range (T_eff = 5600–6000 K) close to the solar temperature are plotted in Figure 15(b), while those of late G-type main-sequence stars (T_eff = 5100–5600 K) are in Figure 15(a). The stars showing superflares tend to have shorter rotation periods (younger ages) and larger starspot areas.

Zoom In Zoom Out Reset image size

Figure 15(b) shows that the largest area of starspots on Sun-like stars with T_eff = 5600–6000 K in a given period bin has roughly a constant or very gentle decreasing trend around A_spot = 5 × 10⁻²–1 × 10⁻¹A_1/2⊙ in the period range of P_rot ≲ 12 days (age: t ≲ 1.4 Gyr). However, in the period range of P_rot ≳ 12 days (t ≳ 1.4 Gyr), the largest starspot area on them clearly decreases as the rotation period increases. As for the late G-type main-sequence stars with T_eff = 5100–5600 K in Figure 15(a), the similar steep trends can be seen for maximum size of starspots, but the exact values are a bit different. This steep decreasing trend starts around P_rot ∼ 14 days (Figure 15(a)). In the longer period range, the maximum area of spots on the stars with T_eff = 5600–6000 K at a given period bin is roughly half compared with those of the stars with T_eff = 5100–5600 K (see the black dashed line in Figure 15(b)). These differences support that cooler stars can generate larger magnetic flux more effectively thanks to the development of the convection zone. We discuss such temperature differences of the P_rot versus A_spot relations including more cooler (K- and M-type) dwarf stars in more detail in our next study (H. Maehara et al. 2019, in preparation).

The trends similar to those presented here were already reported in Figure 1(a) of Maehara et al. (2017), but we did not discuss the temperature differences. Moreover, the steep decreasing trend in the longer period range is now much clearer compared with this previous paper. It might be because the errors of A_spot values decrease thanks to the updates using Gaia-DR2 stellar radius (R_Gaia) values, and the potential contamination of subgiants in the previous data can also be eliminated. Then, the steep decreasing trend of maximum spot area in the longer period range is more strongly supported.

In Section 5.3, we reported that maximum energy of superflares in a given period bin decreases as P_rot increases (Figures 13(a) and (b)). This decreasing trend of maximum flare energy can be related to the decreasing trend of maximum area of starspots in the longer-period regime in Figures 15(a) and (b) described in this section. This is because the maximum area of starspots determines well the upper limit of flare energy (see Section 5.2). For example, Figure 15(b) shows that the maximum size of starspots on old (t ∼ 4.6 Gyr) Sun-like stars with T_eff = 5600–6000 K and P_rot ∼ 25 days is ∼1% of the solar hemisphere. This corresponds to flare energy of 10³⁴–10³⁵ erg on the basis of Equation (5), and the upper limit of superflare energy of such Sun-like stars in Figure 13(b) is roughly in the same range.

However, there is a difference in a bit strict sense between the decreasing trend of the maximum superflare energy in Figures 13(a) and (b) and that of the maximum area of starspots in Figures 15(a) and (b). The maximum superflare energy roughly continuously decreases as the rotation period increases (the star becomes older) in Figures 13(a) and (b), but the maximum area of starspots does not show such continuous decreasing trend in Figures 15(a) and (b). The maximum area of starspots is constant in the short-period range (P_rot ≲ 12–14 days), but it steeply decreases as the period increases in the longer range (P_rot ≳ 12–14 days). This difference can suggest a possibility that the flare energy is determined not only by the starspot area but also by the other important factors, though the starspot area is still a necessary condition determining flare energy (see Section 5.2). By analogy with the correlation between the flare activity and the magnetic structure of sunspot groups (see Sammis et al. 2000; Maehara et al. 2017), one of the possible factors might be the effect of the magnetic structure of starspots. More complex spots can generate more frequent and more 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 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 (see also the final paragraph of Maehara et al. 2017). We must note here that we also need more superflare observations to more quantitatively discuss the above difference between Figures 13 and 15, since the number of superflare events, especially in Figure 13(b), is small.

In addition, the existence of constant and decreasing trends of maximum starspot coverage can be compared with the relation between SXR flux and rotation period (e.g., Wright et al. 2011). The SXR 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 SXR trend (P_rot ∼ 2–3 days) is different from that of maximum spot size values (P_rot ∼ 12–14 days). Although a detailed study on this point is beyond the scope of this paper, these similarities and differences can be interesting and helpful when considering the relation between stellar activity (including starspots, flares, X-ray steady emissions) and rotation in more detail.

Figure 16 represents the occurrence frequency distributions of superflares on solar-type stars with T_eff = 5600–6000 K ("the temperature range close to the solar value" in the above sections) and various rotation period (P_rot) values, derived from the superflare data updated earlier in this study. As also presented in our previous studies (e.g., Shibayama et al. 2013) with the original data of superflares, we can see that there are the power-law distributions (dN/dE ∝ E^−α with the index α ≲ 2), and rapidly rotating stars tend to have larger frequency values. The upper-limit values of flare energy roughly depend on rotation period, as already seen in Figure 13. However, error values are relatively large, especially for slowly rotating stars, because of the small number of analyzed data. For example, one blue circle in the energy range E_flare = 10^35.75–10^36.0 erg in Figure 16, which is a bit far from the other blue circles, corresponds to only one superflare event. We must note here that these values are treated with caution, and we need to increase the number of slowly rotating superflare stars by including more new data in the future. We also note here that, as in Figure 14, the flare frequency value in the vertical axis of Figure 16 is calculated by using the number of solar-type stars in each P_rot bin (N_P(P_rot) in Table 9 in Appendix B) detected in McQuillan et al. (2014). The potential errors from this N_P(P_rot) value, which are discussed in Appendix B, should be kept in mind when we discuss the dependence of flare frequency on the rotation period on the basis of Figure 16. For example, flare frequency of Sun-like stars with P_rot = 20–40 days could become a factor of two smaller in the largest error case, and more investigations with new data in the future are strongly needed.

Zoom In Zoom Out Reset image size

Our previous papers (Shibayama et al. 2013; Maehara et al. 2015) pointed out that the frequency distribution of superflares on Sun-like stars and those of solar flares are roughly on the same power-law line. However, the definition of Sun-like stars in these previous studies (T_eff = 5600–6000 K and P_rot > 10 days) includes many stars much younger than the Sun. For example, stars with P_rot ∼ 10 days have an age of t ∼ 1 Gyr (see Figures 13(b) and 15(b)). It might be better to use only the data of stars rotating as slowly as the Sun (P_rot ∼ 25 days and t ∼ 4.6 Gyr). Then, in Figure 17 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 16, in addition to the data of solar flares and superflares shown in Figure 4 of Maehara et al. (2015). As a result, the newly added data point of superflares on Sun-like stars is roughly on the same power-law line, though the exact value of superflare frequency of stars with P_rot = 20–40 days (t > 3.2 Gyr) is a bit smaller than those of stars with P_rot > 10 days (t > 1 Gyr). From this figure, we can roughly remark that superflares with energy >10³⁴ erg would be approximately once every 2000–3000 yr on old Sun-like stars with P_rot ∼ 25 days and t ∼ 4.6 Gyr, though the error value is relatively large because of the small amount of data of slowly rotating stars.

Zoom In Zoom Out Reset image size

We note here 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 recent 1000–2000 yr (e.g., Usoskin 2017 for review). For example, significant radioisotope ¹⁴C enhancements have been reported in tree rings for the years AD 775 and AD 994, and they suggest extremely strong and rapid cosmic-ray increase events possibly caused by extreme solar flares (e.g., Miyake et al. 2012, 2013). Various potential low-latitude auroral drawings have also been reported from the historical documents around the world, and they suggest the possibility that extreme solar flare events have occurred several times in the recent ∼1000 yr (e.g., Hayakawa et al. 2017a, 2017b). In the future studies, the superflare frequency information as in Figure 17 should be compared with these radioisotope and historical data in detail.

In addition to the relation between the rotation period and starspot area in Section 5.4, in the following we also investigate the size–frequency distribution of large starspot groups on Sun-like stars and that of sunspots. Our previous paper (Maehara et al. 2017) already conducted this analysis, but here we aim to investigate again whether both sunspots and larger starspots can be related to the same physical processes, by including the updates using Gaia-DR2 stellar radius values.

Figure 18 shows the comparison between the cumulative appearance frequency of starspots on Sun-like stars (T_eff = 5600–6000 K and P_rot = 20–40 days) and that of sunspot groups. The estimation method of the cumulative appearance frequency values is basically the same as that used for Figure 5 of Maehara et al. (2017), and this was already described in detail in that paper. Only one difference is that we use the updated A_spot values described in the above Section 5.4. The occurrence frequency decreases as the area of sunspots or starspots increases. The appearance rate of both the sunspot groups and that of starspots with the area >10^−2.5A_1/2⊙ is approximately once in a few years. The cumulative appearance frequency of starspots on Sun-like stars can be fitted by a power-law function with the power-law index of −2.9 ± 0.5 (red dashed line) for spot areas between 10^−2.5A_1/2⊙ and 10^−1.0A_1/2⊙. This power-law index is a bit steeper compared with those originally presented in Maehara et al. (2017).

Zoom In Zoom Out Reset image size

According to Bogdan et al. (1988), the size distribution of individual sunspot umbral areas shows lognormal distribution. Although the overall size–frequency distribution of sunspot groups in Figure 18 also shows a similar lognormal distribution, the distribution of sunspot groups for large sunspots is roughly on this power-law line for sunspot areas between 10^−3.25A_1/2⊙ and 10^−2.5A_1/2⊙. The appearance frequency of sunspots with spot areas of ∼10⁻²A_1/2⊙ is about 10 times lower than that of starspots on Sun-like stars. This 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). These results shown here are basically the same as those originally presented in Maehara et al. (2017), and we confirmed the similarity between the size distribution of sunspots and that of starspots. This supports that both sunspots and larger starspots could be related to the same physical processes. The upper limit of starspot size values of Sun-like stars in Figure 18 is approximately a few × 10⁻²A_1/2⊙, and the appearance frequency of these spots is approximately once every 2000–3000 yr.

In Sections 2–4, we described the results of our spectroscopic observations of Kepler solar-type superflare stars. We newly conducted APO 3.5 m telescope spectroscopic observations of the 18 solar-type superflare stars found from Kepler 1-minute (short) time-cadence data (Maehara et al. 2015).

• (i)  
  More than half (43 stars) are confirmed to be "single" stars, among the 64 solar-type superflare stars in total that have been spectroscopically investigated so far in this APO 3.5 m observation and our previous Subaru/HDS (Notsu et al. 2015a) observations (Table 6).
• (ii)  
  The "v sin i" (projected rotational velocity) values are consistent with the rotational velocity values v_lc estimated from the brightness variation period of Kepler data (Figure 2).
• (iii)  
  There is a positive correlation between the amplitude of the brightness variation and Ca ii (Ca ii λ8542 and Ca ii H and K lines) index values (Figure 3), and this suggests that there is a rough positive correlation between the starspot coverage from Kepler photometric data and the stellar average magnetic field.
• (iv)  
  The results of (ii) and (iii) support the idea that the quasi-periodic brightness variation of Kepler solar-type superflare stars is caused by the rotation with large starspots.
• (v)  
  Lithium abundances of superflare stars suggest that Kepler solar-type superflare stars include many young stars but also include old stars like our Sun (Figure 4).

Then in Section 5, we investigated the statistical properties of Kepler solar-type superflare stars originally described in our previous studies (Maehara et al. 2012, 2015, 2017; Notsu et al. 2013b; Shibayama et al. 2013), by incorporating Gaia-DR2 stellar radius estimates (reported in Berger et al. 2018) and updating the parameters (e.g., flare energy, spot size).

• (i)  
  More than 40% of the original solar-type superflare stars in our previous studies are now classified as subgiant stars (Figure 5).
• (ii)  
  The bolometric energy released by flares (E_flare) is not inconsistent with the magnetic energy (E_mag) stored around the large starspots (Figures 6(b) and 11).
• (iii)  
  Our previous studies suggested that the maximum superflare energy in a given rotation period bin does not have a clear correlation with the rotation period, and superflares up to 10³⁵ erg could occur on slowly rotating solar-type stars (Figure 12(a)). This study suggests that the maximum superflare energy continuously decreases as the rotation period increases (as the star becomes older; Figures 12(b) and 13). Superflares with energies ≲5 × 10³⁴ erg occur on old, slowly rotating Sun-like stars (T_eff = 5600–6000 K, P_rot ∼ 25 days, and age t ∼ 4.6 Gyr) approximately once every 2000–3000 yr (Figures 13(b) and 17). In contrast, superflares up to ∼10³⁶ erg can occur on young, rapidly rotating stars (P_rot ∼ a few days and t ∼ a few hundred megayears; Figure 13), and the flare frequency of such young, rapidly rotating stars is ∼100 times higher compared with the above old, slowly rotating Sun-like stars (Figures 14 and 16).
• (iv)  
  The maximum area of starspots does not depend on the rotation period and has a roughly constant or very gentle decreasing trend around A_spot = 5 × 10⁻²–1 × 10⁻¹A_1/2⊙ (A_1/2⊙ ∼ 3 × 10²² cm²: solar hemisphere) when the star is young and rapidly rotating. However, as the star becomes older and its rotation slows down, it starts to have a steep decreasing trend at a certain point: P_rot ∼ 12 days (t ∼ 1.4 Gyr) for the stars with T_eff = 5600–6000 K (Figure 15(b)), and P_rot ∼ 14 days for the stars with T_eff = 5100–5600 K (Figure 15(a)). The maximum size of starspots on slowly rotating Sun-like stars is ∼1% of the solar hemisphere, and this is enough for generating superflares with the energy ≲5 × 10³⁴ erg described in (iii).
• (v)  
  These decreasing trends of the maximum flare energy (in (iii)) and the maximum starspot area (in (iv)) can be related to each other since the superflare energy can be explained by the starspot magnetic energy as in (ii). However, there is also a difference between the two: the maximum area of starspots starts to steeply decrease at a certain P_rot value (as in (iv)), while the maximum flare energy continuously decreases as the rotation slows down (as in (iii)). This can suggest a possibility that the flare energy is determined not only by the starspot area but also by other important factors (e.g., spot magnetic structure).
• (vi)  
  The size distribution of starspots on Sun-like stars (T_eff = 5600–6000 K and P_rot ∼ 25 days) between A_spot = 10^−2.5A_1/2⊙ and 10^−1.0A_1/2⊙ is roughly located on the extension line of the distribution of sunspot groups between A_spot = 10^−3.25A_1/2⊙ and 10^−2.5A_1/2⊙ (Figure 18). The upper limit of starspot size values on slowly rotating Sun-like stars would be ∼a few × 10⁻²A_1/2⊙, and the appearance frequency of these spots is approximately once every 2000–3000 yr.

In this paper, we have reported the current updates of our series of studies on superflares on solar-type stars. However, we need more studies to clarify the properties of superflare stars on Sun-like stars and to answer the important question, "Can our Sun have superflares?" For example, our spectroscopic observations so far have observed 64 solar-type superflare stars, but the number of Sun-like superflare stars (T_eff = 5600–6000 K, P_rot ∼ 25 days, and t ∼ 4.6 Gyr) that have been investigated spectroscopically and found to be "single" stars are now only a few (∼1). In the figures in Section 5, the number of old, slowly rotating Sun-like superflare stars are now very small, and the current statistical discussions are not enough. In this study, as for superflares found from Kepler 30-minute cadence data, we only used the data originally found as solar-type superflare stars from Kepler data of the first ∼500 days (Quarter 0–6) in our previous study (Shibayama et al. 2013; Candelaresi et al. 2014). In our next study, we plan to increase the amount of data by using the whole Kepler data set of ∼1500 days (Quarter 0–17).²⁸ Moreover, the data from the next missions such as TESS (Ricker et al. 2015) and PLATO (Rauer et al. 2014) are expected to bring us more superflare data on solar-type stars in the near future. These two missions are helpful also for spectroscopic follow-up observations of superflare stars since these two missions observe more nearby stars (with their wider spatial coverage), compared with Kepler.

In addition to the statistical properties of superflares that we discussed in this paper, there are so many topics that should be investigated in the near future to fully understand the physics of superflares and their effects in the related research fields, for example, mechanisms of white-light continuum emissions of superflares (e.g., Namekata et al. 2017; Heinzel & Shibata 2018; Kowalski & Allred 2018), chromospheric line profiles during superflares (e.g., Houdebine et al. 1993; Hawley et al. 2007; Honda et al. 2018), stellar mass ejections (e.g., CMEs) during superflares (e.g., Aarnio et al. 2012; Leitzinger et al. 2014; Hudson 2015; Osten & Wolk 2015; Harra et al. 2016; Takahashi et al. 2016; Vida et al. 2016; Moschou et al. 2017; Crosley & Osten 2018), the impacts of superflares on planets (e.g., Segura et al. 2010; Airapetian et al. 2016; Atri 2017; Lingam & Loeb 2017; Riley et al. 2018), detailed comparisons with a history of solar activity over ∼1000 yr (e.g., Miyake et al. 2012, 2013; Hayakawa et al. 2017a, 2017b; Usoskin 2017), quasi-periodic pulsations (QPPs) of the brightness during superflares (e.g., Pugh et al. 2016; McLaughlin et al. 2018), complexities of superflare-generating starspots (e.g., Maehara et al. 2017; Toriumi et al. 2017), how starspots are distributed on the surface of superflare stars (e.g., Doyle et al. 2018; Roettenbacher & Vida 2018), how starspots and plages are colocated in the case of active stars (e.g., Morris et al. 2018), and lifetimes and formation/decay processes of large starspots in solar-type stars (e.g., Shibata et al. 2013; Giles et al. 2017; Namekata et al. 2019). The superflare stars newly found from the above new missions (TESS and PLATO) would be helpful for the detailed observations to investigate these topics (e.g., spectroscopic or multiwavelength observations of superflares themselves, and long-term changes of stellar activities).

For the first step of our analyses, we checked the binarity of each superflare star, as we did in Notsu et al. (2015a). First, we checked slit viewer images of the APO 3.5 m telescope. Two stars (KIC 11551430 and KIC 7093428) have visual companion stars, as shown in Figures 19(a) and 20(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As for KIC 11551430, we took spectra of both components (the brighter primary star, KIC 11551430A, and the fainter companion star, KIC 11551430B) of the visual binary system separately, as shown in Figures 19(b) and (c). We checked again pixel count data of superflare events on KIC 11551430 (Figures 19(d) and (e)), which are originally from Maehara et al. (2015), and confirmed that the center positions of the brightness during superflare events are roughly the same as those during quiescent state. This is consistent with the possibility that superflares occur on the primary star KIC 11551430A (see Figures S1 and S2 of Maehara et al. 2012).

As for KIC 7093428, we can see one primary brighter star (KIC 7093428A) and two fainter companion stars (KIC 7093428B and KIC 7093428C) in Figure 20(a). We only took spectra of the main brighter star (KIC 7093428A) as shown in Figure 20(b), since companions stars (B and C) are too faint. We also checked again pixel count data of superflare events on KIC 7093428 (Figures 20(c) and (d)), which are originally from Maehara et al. (2015). We can see that the center positions of the brightness during superflare events are shifted compared with those during quiescent state. This suggests the high possibility that superflares occur on the fainter companion stars KIC 7093428B or KIC 7093428C (see Figures S1 and S2 of Maehara et al. 2012), not on the primary G-type star KIC 7093428A. Then KIC 7093428 has "VB" in the second column of Table 1.

As described in Section 2, we selected 23 target stars and took spectra of 22 stars (minus the too faint target star KIC 10745663; K_p = 14.3 mag). Among the 22 observed stars, the data quality of our spectroscopic data of the fainter four stars (KIC 6032920, KIC 10528093, KIC 10646889, and KIC 9655134) is not high enough for the detailed following discussions. We plot examples of photospheric lines in Figure 21 for reference. We only remark here that they do not show any clear signs of binarity, on the basis of Figure 21 and the slit viewer images of the APO 3.5 m telescope. In the following discussions of this paper, we only treat the 18 target stars that we got spectra of with enough S/N.

Zoom In Zoom Out Reset image size

Next, we investigated the line profiles and found that two stars (KIC 11128041 and KIC 10338279) show double-lined profiles. In this process, we checked by eye the profiles of the many spectral lines, and the double-lined spectra of these two stars are shown in Figure 22. Since these double-lined profiles are caused by the overlap of the radiation of multiple stars, we regard these two stars as double-lined spectroscopic binary stars. These two stars have "SB2" in the second column of Table 1.

Zoom In Zoom Out Reset image size

Third, we investigated time variations of the line profiles between the multiple observations that are expected to be caused by the orbital motion in the binary system. This investigation was for the target stars that we observed multiple times (16 stars). In Section 2, we already measured the RV of all the target stars that were not classified as visual binary stars or double-lined spectroscopic binary stars, and these values are listed in Table 2. We have also conducted spectroscopic observations of these four stars (KIC 4742436, KIC 4831454, KIC 9652680, and KIC 11610797) using the Subaru telescope (Honda et al. 2015; Notsu et al. 2015a, 2015b), and we also use these data in the RV investigation here (they are also listed in Table 2). As a result, KIC 11551430A, KIC 4543412, and KIC 11128041 show RV changes as shown in Figure 22, and these RV changes are large enough compared with RV errors of APO data (≲1 km s⁻¹). We confirmed this typical RV error value (≲1 km s⁻¹) by comparing our spectroscopic data of the single comparison stars that we observed multiple times in this study and are listed in Table 3. These three stars have "RV" in the second column of Table 1.

In total, we regard five superflare stars as binary stars among the 18 target stars of which we got spectra with enough S/N. The remaining 13 superflare stars do not show any evidence of binarity within the limits of our analyses, so we treat them as "single stars" in this paper. They have "no" in the second column of Table 1. Spectra of photospheric lines, including Fe i λλ6212, 6215, 6216, and 6219 lines, of the 13 "single" superflare stars are shown in Figure 23. We observed nine stars among the 13 "single" stars multiple times, and we made co-added spectra of these nine stars by conducting the following two steps. First, we shifted the wavelength value of each spectrum to the laboratory frame on the basis of the RV value of each observation listed in Table 2. Next, we added up these shifted spectra to one co-added spectrum. The co-added spectra are mentioned as "comb" in Table 2, and only co-added spectra of these nine stars are used in Figure 23. Only the co-added spectra are used for the detailed analyses in the following sections of this paper, when we analyze the spectral data of the nine stars that we observed multiple times.

Zoom In Zoom Out Reset image size

In addition to these nine "single" stars, we also made co-added spectra with the same methods for the two binary target stars (KIC 11551430A and KIC 4543412) that show RV shifts but do not show double-lined profiles. They are also shown in Figure 23. For reference, we also use these data when we estimate stellar parameters in the following.

We estimated the effective temperature T_eff, surface gravity log g, microturbulence v_t, and metallicity [Fe/H] of the target superflare stars, by using the method that is basically the same as the one we have used in Notsu et al. (2015a, 2017). We measured the equivalent widths of ∼200 Fe i and Fe ii lines and used the TGVIT program developed by Takeda et al. (2002, 2005). For reference, Rich et al. (2017) also applied this method to their ARCES spectroscopic data, which were taken with the same wavelength resolution value (R ∼ 32,000) as our data, and they confirmed that the resultant values are consistent with the other previous studies. The resultant atmospheric parameters (T_eff, log g, v_t, and [Fe/H]) are listed in Table 4.

We then compare these resultant atmospheric parameters with the values reported in DR25-KSPC (Mathur et al. 2017). We show the results in Figure 24, where the data points are classified with colors on the basis of the methods used to derive the parameters in DR25-KSPC. Our values seem to be in good agreement with the DR25-KSPC values, especially for T_eff (Figure 24(a)). log g has some dispersion (Figure 24(b)), but this small dispersion does not cause essential problems when we consider whether the target stars are main-sequence stars or not. As for [Fe/H], if we only consider the spectroscopic values in DR25-KSPC ("DR25-SPE" in Figure 24(c)), both values have better agreement, and the difference does not affect the overall discussions (e.g., whether the target stars are "metal-rich/poor" or not). Through this, we can remark that our spectroscopically derived values are good sources with which to discuss the actual properties of stars, as also mentioned for our previous Subaru spectroscopic data in Section 4.3 of Notsu et al. (2015a).

Zoom In Zoom Out Reset image size

As already listed in Table 1, we have stellar radius values R_Gaia of Kepler stars deduced from Gaia-DR2 parallax values (Berger et al. 2018). These values are listed again in Table 4, but five stars among our observed 18 target stars have no R_Gaia values reported in Berger et al. (2018).

In Figure 25 we compare these R_Gaia values with log g values estimated from our spectroscopic studies (Notsu et al. 2015a; this study). These two values are roughly correlated, and R_Gaia values look more sensitive to the boundary region between main-sequence stars and subgiants (around log g = 4.0–4.5 and R_Gaia = 1.5–2.0 R_⊙). Then in Section 5, we mainly use these stellar radius R_Gaia values to strictly define main-sequence stars.

Zoom In Zoom Out Reset image size

We also estimated stellar radius (R_spec) values from the stellar atmospheric parameters (T_eff, log g, and [Fe/H]), by applying the latest PARSEC isochrones (Bressan et al. 2012; Marigo et al. 2017).²⁹ In this process, we used basically the same method as the one used in our previous study (Notsu et al. 2015a). We selected all the data points that had possible sets of T_eff, log g, and [Fe/H] from the PARSEC isochrones, taking into account the error values of T_eff and log g (ΔT_eff and Δlog g). For the six stars that have no suitable isochrones within their original error range of T_eff and log g values, we then took into account larger error values as mentioned in footnote c of Table 4. We note that the resultant values of these six stars can have relatively low accuracy. We then selected the maximum and minimum R_spec value of each star and determined the resultant R_spec value as a median of the maximum and minimum values. The error values of R_spec listed in Table 4 are ≲20% for most of the stars, as also mentioned in Notsu et al. (2015a).

In the following, we use R_Gaia values as a first priority, and we use R_spec values only for the stars without R_Gaia values.

We measured v sin i (stellar projected rotational velocity) of the target stars by using the method that is basically the same as that in our previous studies (Notsu et al. 2015a, 2017). This is originally based on the one described in Takeda et al. (2008). We took into account the effects of macroturbulence and instrumental broadening by considering a simple relationship among the line-broadening parameters (see Takeda et al. 2008), which can be expressed as

$$\begin{eqnarray}&&{v}_{{\rm{M}}}^{2}={v}_{\mathrm{ip}}^{2}+{v}_{\mathrm{rt}}^{2}+{v}_{\mathrm{mt}}^{2}.\end{eqnarray} \tag{ 6 }$$

Here v_M is the e-folding width of the Gaussian macrobroadening function $f{(v)\propto \exp [(v/{v}_{{\rm{M}}})}^{2}]$, including instrumental broadening (v_ip), rotation (v_rt), and macroturbulence (v_mt). We derived v_M by applying an automatic spectrum-fitting technique (Takeda 1995), assuming the model atmosphere corresponding to the atmospheric parameters estimated in Appendix A.2. In this process, we used the MPFIT program contained in the SPTOOL software package developed by Y. Takeda. We applied this fitting technique to the 6212–6220 Å region (shown in Figures 22 and 23) to derive v sin i values. This region has also been used in our previous studies (Notsu et al. 2015a, 2017).

The instrumental broadening velocity v_ip was calculated using the following equation (Takeda et al. 2008):

$$\begin{eqnarray}&&{v}_{\mathrm{ip}}=\displaystyle \frac{3\times {10}^{5}}{2R\sqrt{\mathrm{ln}2}},\end{eqnarray} \tag{ 7 }$$

where R(=λ/Δλ) is the wavelength resolution of the observation. For estimating the R value adopted here, we conducted Gaussian-fitting to emission lines in the 6180–6240 Å region of the Th–Ar spectrum data. This region is around the 6212–6220 Å region, where we conduct the above fitting process. We finally got R = 32,500, and we applied this to Equation (7) in the following. The macroturbulence velocity v_rt was estimated by using the relation v_mt ∼ 0.42ξ_RT (Takeda et al. 2008). The term ξ_RT is the radial-tangential macroturbulence, and ξ_RT was estimated using the relation reported in Valenti & Fischer (2005):

$$\begin{eqnarray}&&{\xi }_{\mathrm{RT}}=\left(3.98-\displaystyle \frac{{T}_{\mathrm{eff}}-5770\,{\rm{K}}}{650\,{\rm{K}}}\right).\end{eqnarray} \tag{ 8 }$$

As we described in Notsu et al. (2015a), the choice of macroturbulence equation is often important, but we only use this Equation (8) to estimate the resultant values of this study, as we have also done in our previous studies (Notsu et al. 2015a, 2017). v_rt was then derived with the above equations, and we finally got v sin i with the relation v_rt ∼ 0.94 v sin i (Gray 2005). The resultant v sin i values of the target superflare stars are listed in Table 4. As error values of v sin i, we consider the systematic uncertainty with changing ξ_RT up to ±25%, as described in Hirano et al. (2012).

The above estimation method of v sin i has been developed to be suitable for the Subaru/HDS spectroscopic data with the high spectral resolutions of R ∼ 55,000–100,000 (Hirano et al. 2012; Notsu et al. 2013a, 2015a). It is not so appropriate that we apply this method without any modifications to our APO 3.5 m/ARCES spectroscopic data (the spectral resolution of only ∼32,500). It is difficult to estimate v sin i values as low as 2–3 km s⁻¹ only with APO spectra. On the other hand, in this study we do not need to estimate v sin i values with the high precision of 2–3 km s⁻¹ for the overall discussions in this study with Figures 2 and 9 (e.g., whether v sin i ≲ v_lc is roughly achieved, and whether the target stars can have low inclination angle values or not). There are twice as many (or more) differences of the spectral resolution values between this APO observation and the previous Subaru/HDS observations (Hirano et al. 2012; Notsu et al. 2013a, 2015a), and the instrumental velocity v_ip has twice (or more) larger values between them. Then, in this study (e.g., Table 4) we only report a rough upper-limit value "v sin i < 4 km s⁻¹" for the slowly rotating stars with v sin i < 4 km s⁻¹. We also need to use the v sin i values of mildly slowly rotating stars (v sin i ∼ 5 km s⁻¹) with caution.

The observed spectra of the target superflare stars around Ca ii λ8542 and Hα λ6563 are shown in Figures 26 and 27, respectively. We measured the r₀(8542) and r₀(Hα) indices, which are the residual core fluxes normalized by the continuum level at the line cores of Ca ii λ8542 and Hα λ6563, respectively. As we have already introduced in Notsu et al. (2013a), these indices are known to be indicators of stellar chromospheric activity (e.g., Linsky et al. 1979a; Takeda et al. 2010). As the chromospheric activity is enhanced, the intensity of these indicators becomes large since a greater amount of emission from the chromosphere fills in the core of the lines. The values of the r₀(8542) and r₀(Hα) indices of the target superflare stars are listed in Table 5. For reference, r₀(8542) values of the 28 comparison stars are also listed in Table 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The observed spectra of the target superflare stars around Ca ii H λ3968 and Ca ii K λ3934 lines are shown in Figure 28. As described in Section 2, these spectra are normalized by using spectra of early-type standard stars, as also done in Morris et al. (2017). The emission in the cores of the Ca ii H and K lines is a widely known indicator of the stellar chromospheric activity (e.g., Hall 2008 for review), and this is more sensitive to activity level changes compared with the Ca ii λ8542 line (see Takeda et al. 2010, 2012). Ca ii H and K emission is often measured as "S-index," which is the flux in the emission features normalized by two pseudocontinuum regions on either side of the absorption features (e.g., Vaughan et al. 1978; Duncan et al. 1991; Isaacson & Fischer 2010; Mittag et al. 2013; Karoff et al. 2016; Morris et al. 2017). We measured the S-index values of the target stars in the rest of this section.

Zoom In Zoom Out Reset image size

The S-index value can vary among other instruments for the same intrinsic flux, as described in Section 2.1 of Isaacson & Fischer (2010). Following the method used in Isaacson & Fischer (2010) and Morris et al. (2017), we calibrated S-index values of our APO observation data (S_APO) against the S-index values already calibrated to the Mount Wilson Observatory (MWO) sample (S_MWO).³⁰ In this calibration process, we use our observation data of the 28 bright solar-type comparison stars described in Section 2 and listed in Table 3. All of these 28 stars have S-index values calibrated to the MWO sample (S_MWO values in Table 3) that are reported in Isaacson & Fischer (2010) on the basis of CPS program spectroscopic observations.

We first measured the S-index from our APO data (S_APO) of these comparison stars using the following equation (see Morris et al. 2017):

$$\begin{eqnarray}&&{S}_{\mathrm{APO}}=\displaystyle \frac{{aH}+{bK}}{{cR}+{dV}},\end{eqnarray} \tag{ 9 }$$

where H and K are the recorded counts in 1.09 Å FWHM triangular bandpasses centered on the Ca ii H and K lines at 3968.469 and 3933.663 Å, respectively. V and R are two 20 Å wide reference bandpasses centered on 3901.07 and 4901.07 Å. The values of a, b, c, and d should be selected so that S_APO has roughly equal flux contribution from the H and K emission lines and roughly equal flux contribution from the R and V pseudocontinuum regions in the APO spectra. Here we finally had a = b = c = d = 1 after trial and error in this study.

Then in Figure 29, we compared these APO S-index values (S_APO) listed in Table 3 with the S_MWO values reported in Isaacson & Fischer (2010), and we investigated the two constants of the following equation (see Morris et al. 2017):

$$\begin{eqnarray}&&{S}_{\mathrm{MWO}}={C}_{1}{S}_{\mathrm{APO}}+{C}_{2}.\end{eqnarray} \tag{ 10 }$$

Since S-index varies over time for each star in the sample, the linear correlation between S_APO and S_MWO can have some intrinsic spread as in Figure 29. Applying the least-squares method to the data in Figure 29, finally we got C₁ = 22.005 and C₂ = 0.009. Such larger (e.g., C₁ ∼ 22) conversion factors of S-indices have commonly appeared in studies using a relatively high dispersion spectrograph, for example, Morris et al. (2017) using APO/ARCES data (the same instruments and same method as our study; C₁ ∼ 21.1), Karoff et al. (2013) using NOT/FIES data (R ∼ 25,000 and C₁ ∼ 16.6), and Isaacson & Fischer (2010) using Keck/HIRES data (R ∼ 52,000 and C₁ ∼ 31.5). As described in the final paragraph of Section 3.6 of Karoff et al. (2013), the relatively large C₁ value in our study can be mainly due to the higher spectral resolutions of APO data, compared with the original Mount Wilson observations (Duncan et al. 1991). Using Equations (9) and (10) and these C₁ and C₂ values, we finally estimated the calibrated S-index value of each spectrum (S_HK ≡ S_MWO) from each S_APO value. The S_HK values from the observations of comparison stars are listed in Table 3, and the resultant S_HK values of the target superflare stars from their spectral data (Figure 28) are listed in Table 5.

Zoom In Zoom Out Reset image size

S-index is known to be stellar temperature dependent (e.g., Noyes et al. 1984), and it is a purely empirical quantity. Hence, it can be advantageous to calculate chromospheric fluxes. Then, in the following, we calculated Ca ii H and K flux values from the S-index values with stellar colors, using the formulation described in Mittag et al. (2013). We summarize the method in the following.

The stellar surface fluxes emitted in the Ca ii H and K lines (${{ \mathcal F }}_{\mathrm{HK}}$) are expressed with S-index (S_HK) as

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{HK}}=\displaystyle \frac{{{ \mathcal F }}_{\mathrm{RV}}}{\alpha }{S}_{\mathrm{HK}},\end{eqnarray} \tag{ 11 }$$

where ${{ \mathcal F }}_{\mathrm{RV}}$ is the surface flux in both continua (V and R region mentioned in Appendix A.5) and the factor α is a historical dimensionless conversion factor. Following Mittag et al. (2013), we use the value α = 19.2. ${{ \mathcal F }}_{\mathrm{RV}}$ is estimated from stellar color B − V using the following empirical equation applicable to main-sequence and subgiant stars:

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{{ \mathcal F }}_{\mathrm{RV}}}{19.2}\right)=8.25-1.67(B-V).\end{eqnarray} \tag{ 12 }$$

As for the target Kepler superflare stars in Table 4 and the 15 comparison stars (among all the 28 comparison stars) with remarks (1) or (2) in Table 3, which were observed in Notsu et al. (2017), B − V values were calculated from T_eff and [Fe/H] using Equation (2) of Alonso et al. (1996). As for the remaining 13 stars, B − V values reported in Table 1 of Isaacson & Fischer (2010) are used, and T_eff values used in the following are derived from B − V values by using Equation (2) of Valenti & Fischer (2005).

We then estimated chromospheric excess flux of Ca ii H and K lines (${{ \mathcal F }}_{\mathrm{HK}}^{{\prime} }$) by subtracting the photospheric flux contribution (${{ \mathcal F }}_{\mathrm{HK},\mathrm{phot}}$) in the line center:

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{HK}}^{{\prime} }={{ \mathcal F }}_{\mathrm{HK}}-{{ \mathcal F }}_{\mathrm{HK},\mathrm{phot}}.\end{eqnarray} \tag{ 13 }$$

${{ \mathcal F }}_{\mathrm{HK},\mathrm{phot}}$ were derived from B − V by using the following equations in Mittag et al. (2013):

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{{ \mathcal R }}_{\mathrm{HK},\mathrm{phot}} & = & -4.898+1.918{\left(B-V\right)}^{2}\\ & & -\ 2.893{\left(B-V\right)}^{3}\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{HK},\mathrm{phot}}=\displaystyle \frac{{{ \mathcal F }}_{\mathrm{HK},\mathrm{phot}}}{\sigma {T}_{\mathrm{eff}}^{4}}.\end{eqnarray} \tag{ 15 }$$

${{ \mathcal F }}_{\mathrm{HK}}^{{\prime} }$ is frequently converted to the flux-related stellar activity index $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ (Linsky et al. 1979b):

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{HK}}^{{\prime} }=\displaystyle \frac{{{ \mathcal F }}_{\mathrm{HK}}-{{ \mathcal F }}_{\mathrm{HK},\mathrm{phot}}}{\sigma {T}_{\mathrm{eff}}^{4}}=\displaystyle \frac{{{ \mathcal F }}_{\mathrm{HK}}^{{\prime} }}{\sigma {T}_{\mathrm{eff}}^{4}}.\end{eqnarray} \tag{ 16 }$$

With this normalization, we can compare the activity level of stars with different effective temperatures and colors. The resultant values of ${{ \mathcal F }}_{\mathrm{HK}}^{{\prime} }$ and ${{ \mathcal R }}_{\mathrm{HK}}^{{\prime} }$ of the target superflare stars are listed in Table 5.

Schrijver (1987) introduced the concept of "basal flux" (${{ \mathcal F }}_{\mathrm{HK},\mathrm{basal}}$), which is residual flux remaining in the core of the Ca ii H and K lines of inactive stars after their photospheric line contribution (${{ \mathcal F }}_{\mathrm{HK},\mathrm{phot}}$) is removed. Mittag et al. (2013) investigated ${{ \mathcal F }}_{\mathrm{HK}}^{{\prime} }$ values of a very large number of main-sequence, subgiant, and giant stars and measured the "lower" boundary of the ${{ \mathcal F }}_{\mathrm{HK}}^{{\prime} }$ distribution as a function of B − V. Mittag et al. (2013) then reported the empirical scaling between ${{ \mathcal F }}_{\mathrm{HK},\mathrm{basal}}$ and B − V:

$$\begin{eqnarray}&&\mathrm{log}{{ \mathcal F }}_{\mathrm{HK},\mathrm{basal}}=a+b(B-V)+c{\left(B-V\right)}^{2},\end{eqnarray} \tag{ 17 }$$

where coefficients a, b, and c vary depending on B − V and luminosity classes (V, IV, or III), and they are listed in Table 3 of Mittag et al. (2013). By subtracting this basal flux component ${{ \mathcal F }}_{\mathrm{HK},\mathrm{basal}}$ from ${{ \mathcal F }}_{\mathrm{HK}}^{{\prime} }$, we can get new "pure" and universal activity indicators ${{ \mathcal F }}_{\mathrm{HK}}^{+}$ and $\mathrm{log}{{ \mathcal R }}_{\mathrm{HK}}^{+}$ defined by

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}_{\mathrm{HK}}^{+} & = & \displaystyle \frac{{{ \mathcal F }}_{\mathrm{HK}}^{{\prime} }-{{ \mathcal F }}_{\mathrm{HK},\mathrm{basal}}}{\sigma {T}_{\mathrm{eff}}^{4}}\\ & = & \ \displaystyle \frac{{{ \mathcal F }}_{\mathrm{HK}}-{{ \mathcal F }}_{\mathrm{HK},\mathrm{phot}}-{{ \mathcal F }}_{\mathrm{HK},\mathrm{basal}}}{\sigma {T}_{\mathrm{eff}}^{4}}\\ & = & \ \displaystyle \frac{{{ \mathcal F }}_{\mathrm{HK}}^{+}}{\sigma {T}_{\mathrm{eff}}^{4}}.\end{array}\end{eqnarray} \tag{ 18 }$$

This $\mathrm{log}{{ \mathcal R }}_{\mathrm{HK}}^{+}$ index allows comparisons of the activity levels of stars with different luminosity classes and different temperatures on the same scale (Mittag et al. 2013). The resultant values of ${{ \mathcal F }}_{\mathrm{HK}}^{+}$ and ${{ \mathcal R }}_{\mathrm{HK}}^{+}$ of the target superflare stars and the comparison stars are listed in Tables 5 and 3, respectively.

In Figures 30(a)–(c), we compared the resultant values of the Ca ii H and K indices (S_HK index, $\mathrm{log}{{ \mathcal F }}_{\mathrm{HK}}^{+}$ index, and $\mathrm{log}{R}_{\mathrm{HK}}^{+}$ index) with those of the r₀(8542) index (normalized intensity at the center of the Ca ii λ8542 line). As seen in these figures, Ca ii λ8542 index and Ca ii H and K indices have good correlations also for Kepler solar-type superflare stars, as expected from the previous studies (e.g.,Takeda et al. 2012; Karoff et al. 2016). In particular, the $\mathrm{log}{R}_{\mathrm{HK}}^{+}$ index (Figure 30(c)) has the best correlation with the Ca ii λ8542 index among these three Ca ii H and K indices. As for this $\mathrm{log}{R}_{\mathrm{HK}}^{+}$ index, stellar temperature (color) dependence and the contribution from photospheric/basal fluxes are removed (Equation (18)). In addition, Figure 30(c) might suggest that the $\mathrm{log}{R}_{\mathrm{HK}}^{+}$ index is much more sensitive to the difference in the lower activity level region (r₀(8542) ≲ 0.3), compared with the r₀(8542) index, as also suggested in Takeda et al. (2012). Thus, we mainly use this $\mathrm{log}{R}_{\mathrm{HK}}^{+}$ index when we discuss the measurement results of Ca ii H and K lines in Section 4.4.

Zoom In Zoom Out Reset image size

The observed spectra of the target superflare stars around the Li i λ6708 line are shown in Figure 31. We measured the Li abundances, A(Li), of the stars using these spectra and the atmospheric parameters estimated in Appendix A.2. We used an automatic profile-fitting method that is basically the same as in our previous studies (Honda et al. 2015; Notsu et al. 2017). This is originally based on the method described in Takeda & Kawanomoto (2005). We summarize the method in the following.

Zoom In Zoom Out Reset image size

In the process of calculating A(Li), we used the MPFIT program contained in the SPTOOL software package also used above. We assumed local thermodynamic equilibrium (LTE) and derived Li abundances using the synthesis spectrum with interpolated model atmospheres taken from Kurucz (1993). We also assumed ⁶Li/⁷Li = 0 throughout this study. The line data around the Li i λ6708 region adapted here are the same as those used in Section 2.2 of Takeda & Kawanomoto (2005). The estimated A(Li) values of the target stars are listed in Table 4. For the stars where Li features are absent (e.g., below the detectable limit), we estimated the upper-limit values of A(Li) by applying the method that we described in Honda et al. (2015) and Notsu et al. (2017).

In Honda et al. (2015) and Notsu et al. (2017) we discussed the typical errors of Li abundances by considering errors arising from multiple causes (errors linked to atmospheric parameters, uncertainties arising from profile-fitting errors, and non-LTE effects). Since we conduct basically the same analyses as the above two papers, we here roughly assume that errors of A(Li) are ∼0.15 dex, as in Notsu et al. (2017).