Temporal Evolution of Spatially-Resolved Individual Star Spots on a Planet-Hosting Solar-type Star: Kepler 17

The first method is a local minima tracing method which is firstly proposed by Hall & Henry (1994) to measure the temporal evolutions of star spots. In this method, each star spot can be identified by the repetition of the local minima over the rotational phases. The light curve of a rotating star with star spots shows several local minima when the spots are on the visible side (Notsu et al., 2013; Maehara et al., 2017). By tracing the local minima with the rotational period of the star, we can estimate the lifetimes of the large-scale spot patterns. We can estimate the temporal evolution of star spot coverage by measuring the local depth of the local minima from the nearby local maxima as an indicator of spot area (hereafter local-minima spot; see Namekata et al., 2019, for the detailed method). We first obtained the smoothed light curve by using locally weighted polynomial regression fitting (LOWESSFIT; Cleveland, 1979) to remove flare signature and noise. In the LOWESSFIT algorithm, a low degree polynomial is fitted to the data subset by using weighted least squares, where more weight is given to nearby points. We used the $lowess$ function incorporated in the statsmodels $python$ package. We detected the local minima as downward convex points of the smoothed light curve; i.e., the smoothed stellar fluxes $F(t)$ satisfy $F(t_{(\rm n-m)})$ $\rm<$ $F(t_{(\rm n-m-1)})$ and $F(t_{(\rm n+m)})$ $\rm<$ $F(t_{(\rm n+m+1)})$. Here, $m$ takes a value of [0, 1, 2], $t$ is time, and $n$ is time step. When we estimate the spot area from the local depth of the light curve, we assume that the spot contrast is 0.3. The advantage of this method is its low computational cost, so it is suitable for statistical analyses, as discussed in Namekata et al. (2019). This star is a good candidate to evaluate how the simple method can estimate the temporal evolution of star spot areas by comparing with the other methods.

Light curve modeling methods for a rotating, spotted star have been developed by several authors (e.g., Croll, 2006; Fröhlich et al., 2012; Lanza et al., 2014). Here, we also developed a light curve modeling method under some simple assumptions (see Ikuta et al., 2019, in prep. for details). The spot contrast is assumed to be constant (=0.3) because it is highly degenerate with spot area (Walkowicz et al., 2013). We are interested in the spot evolution, so each spot area is assumed to emerge and decay linearly for simplicity. The real spot evolution may be more complex, but it will not give a significant effect on our results because we will focus on the relations of spot maximum size and mean variation rate as presented by Namekata et al. (2019). A more realistic spot evolution model should be done in our future works.

In the left panel of Figure 1, modeled stellar light curves are plotted in several cases of spot size and latitude. In the model, we use the stellar surface model separated with grids by assuming linear-limb darkening (see Notsu et al., 2013). As one can see, the stellar light curve of a low-latitude small spot is highly degenerate with those of high-latitude large spots within the Kepler photometric errors of 0.1% (see also Walkowicz et al., 2013). By considering this, we did not consider the latitude information, so there should be uncertainty in area estimations caused by the projection effect at higher latitudes.

There are analytical models which reproduce stellar light curves from the spot parameters (e.g., Kipping, 2012). In this study, we adopt an approach to use a Gaussian function to approximate stellar rotational light curves. As in Figure 1, all of the light curve in the different cases show similar light curves in Kepler 17, which can be fitted by a Gaussian function as indicated in black line, with the Kepler’s photometric errors ($\sim$ 0.1 %). We used the following function as a standard light curve:

where $A$ is spot area, $\sigma$ is non-dimensional factor (=0.110) derived by the Gaussian fitting, $t$ is time, $t_{0}$ is the standard time when the spot are on the visible side, and $P_{\rm rot}$ is the rotational period of a given spots. The right of Figure 1 is an example of the application of this Gaussian light curve to the multi-spot case. In this case, the model light curve $\Delta F(t)$ can be obtained as

Here, $A(t)$ is a spot area as a function of time $t$, and $n$ is the spot number. Note that our original code described by Ikuta et al. (2019) can generate a more realistic light curve considering the inclination and spot latitude.

Here, the following six parameters are necessary for each spot: maximum area, emergence rate, decay rate, maximum timing, standard longitude ($t_{0}$), and rotational period ($P_{\rm rot}$). The total number of parameters is six times the number of spots (here, the number of spots is set to be five). We carried out a parameter search to estimate the most-likely parameters well reproducing the observed light curve with Markov chain Monte Carlo (MCMC). MCMC methods have become an important algorithm in not only astronomical (Sharma, 2017) but also various scientific fields (e.g., Liu, 2001). It can generate samples that follow a posterior distribution by selecting the sampling way based on the likelihood function $L$ in the Bayesian theorem. Here, we take the likelihood function as a product of Gaussian function on each time:

We adopt the traditional Metropolis-Hastings algorithm (Metropolis et al., 1953; Hasting, 1970). We also adopt the Gaussian function as the proposal distribution, and the step length (proposal variance) of the Gaussian was adaptively-tuned for each step by considering the acceptance ratio of MCMC sampling converged to be 0.2 (Gilks et al., 1998; Araki & Ikeda, 2013; Yamada et al., 2017, so-called adaptive MCMC algorithm, e.g.,). A uniform distribution is adopted as a prior distribution.

In our case, when sampling from a multi-modal distribution with a large number of parameters, the chains theoretically can converge to the posterior distribution, but practically seapking, the chains may not converge because of the limited sampling times. They can be trapped in a local mode for a very long time. In order to avoid this, we apply the Parallel Tempering (PT) algorithm to our MCMC estimations (e.g. Hukushima & Nemoto, 1996; Araki & Ikeda, 2013). The PT algorithm introduces auxiliary distribution by using the so-called inverse temperature ($\beta$), runs multiple MCMC samplings (hereafter “replica”) simultaneously for each inverse temperature, and sometimes, exchange replica with a certain percentage. The tempering procedure and exchange processes help samples to escape from a local maximum. In practice, each Markov chain is controlled by the inverse temperature $\beta_{i}$ ($1=\beta_{0}>\beta_{1}>\beta_{2}...>\beta_{n}>0$), and the auxiliary likelihood for each replica is expressed as $L_{i}^{\beta_{i}}$. The higher-order replica has, therefore, the smaller valley to another local maximum than the lower-order one, and the parameters are easy to be searched across the wide range of the parameter space without being trapped in local maxima. By repeating the exchange between the orders, the highest-inverse-temperature replica (= the target samples) can sample around the maximum likelihood within finite computation times.

We run the MCMC sampler for 500,000 steps after the exchanges no longer occur, and set the burn-in period as the first 100,000 steps. As noted above, the replica change can occur theoretically even when sampling around the highest likelihood, but actually did not. We finally take the parameter showing the highest likelihood value in the 400,000 steps, and the error bars are estimated from the posterior distributions. For the parameter range, we take [0.001, 0.05] for maximum spot fraction ($A_{\rm max}$), [-5, -1] for emergence rate ($\Delta A$/day) in log scale, [-5, -1] for decay rate ($\Delta A$/day) in log scale, [start time, end time] for maximum timing, [0^∘, 360^∘] for standard Carrington longitude, and [11.63, 12.86] for rotational period ($P_{\rm rot}$). In this study, we adopt a four-spot model by considering the model evidence and its output. We also carried out five spot model and no significant new information was obtained. Initial input values are randomly selected in the parameter range and independently selected for each replica. The convergence was checked by visually checking the parameter changing for each step and marginalized posterior distribution with one parameter or two parameters. The MCMC method can estimate the parameters to best reproduce the light curve, but the reproduced spots do not necessarily have much to do with the spots that are actually present as we describe in the following paragraph.

We modeled the spots occulted during the exoplanet transit (hereafter in-transit spots) by using STSP code²²2The source code of STSP can be seen here: https://github.com/lesliehebb/stsp (Davenport, 2015; Morris et al., 2017). The spot modeling approach applied in the present study is the same as that already described in Chap. 5 of Davenport (2015) and L. Hebb et al. in preparation, to which we refer the reader for details. In brief, STSP is a C code for quickly modeling the variations in the flux of a star due to circular spots, both in and out of the path of a transiting exoplanet. Static star spots are currently assumed in the code, and evolution of spot activity is typically analyzed by modeling many windows (e.g. individual transits) independently. STSP also allows the user to fit light curves with arbitrary sampling for the maximum-likelihood spot positions and size parameters using MCMC, and has been used for modeling systems with no transiting exoplanets (Davenport et al., 2015), as well as transiting systems with a range of geometries (e.g. Morris et al., 2017). Note that this application of MCMC algorithm is different from that described in the section 3.2 (see, Morris et al., 2017). The MCMC technique allows us to properly explore the degeneracies between star spot positions and sizes. The modeled data is mostly taken from the work by Davenport (2015). Please note that the model by Davenport (2015) was originally carried out by assuming by the stellar radius of 1.1 $R_{\odot}$ and planet radius of 0.10 $R_{\rm star}$, which is slightly different from reported from the other study (Désert et al., 2011; Bonomo et al., 2012; Valio et al., 2017). This can lead to 10 % errors in the estimated area and variation rates of area, but these do not significantly affect our results, because we will discuss much trends that are much larger than this difference. For the STSP outputs, we remove solutions of spots which are located on the stellar limb (i.e. distance from the disk center $>$ 1 - 2 $R_{\rm planet}$) and spots with the too large size (i.e., radius $>$ 0.25 $R_{\rm star}$) because the solutions in the ranges are likely to be biased. We also calculated the error bars of the estimated spot parameters based on the posterior distributions of MCMC samplings.

Figure 3 shows the result of local-minima spots and a comparison with in-transit spots. The middle panel shows the temporal evolution of local-minima spots. As described by Lanza et al. (2019), there are two active longitudes in this star as indicated with the different colors, and spot group Group A (red) is dominant compared to the spot group B (blue) in this period. We use this notation because local light curve minima on the Sun are sometimes due to spot groups or “active longitudes” (e.g., Berdyugina, 2005; Lanza et al., 2019; Namekata et al., 2019). In general these terms may be misleading because the source of light curve minima might not be either of these phenomena. Minima can easily arise (even on the Sun) from a scattered configuration where more spots are on one hemisphere than the other. We have chosen to use conventional terminology, but one can more accurately read “spot group A” as “side A”. The temporal evolution of the area of spot group A shows a clear emergence and decay phase, while spot group B does not.

The bottom panel of Figure 3 is a comparison of spot locations between local-minima spots and in-transit spots. Two local minima are dominant in one rotational phase, while there are 5-6 spots are visible in the transit path. The longitude of the local minima well matches those where in-transit spots are crowded, especially in the first 200 days (BJD 350 to 550+245833). On the other hand, after BJD 550+245833, the red and blue circles lose track of most of the transit spots, and the longitude of local minima becomes delayed compared to the in-transit spots (i.e., the slower rotation period than that of the equator). It could be due to the complex spot distribution changing, but could be caused by the new emergence of star spot on the higher latitude where spots are rotating slower because of the solar-like differential rotation of the star (see, Lanza et al., 2019).

Figure 4 indicates the comparison of the area of local-minima spots and in-transit spots. The upper left panel shows a comparison between Kepler light curve and the reconstructed one from the in-transit spots. Note that the reconstructed light curve in Figure 4 is subtracted by its mean value, and the zero level for the calculation of the reconstructed light curves is around ($F$-$F_{\rm av}$)/$F_{\rm av}$ $\sim$ 0.28. The differential amplitude of the observed light curve is comparable to or a little smaller than the one reconstructed from the spots in-transit.

The lower left panel of Figure 4 shows the comparison of temporal evolutions of the local-minima spot area (black), maximum visible area of in-transit spot (green), and total spot coverage in transit path on half the star (orange). The total spot coverage in the transit chord was calculated for the closest transits to the local minima. The right panel of Figure 4 is the comparison for each data point of the lower left panel of Figure 4. The area of local-minima spots is consistent with the maximum spot area in-transit. The area of local-minima spots also have a positive correlation with the total spot area in transit, although the former is smaller by a factor of two. This tendency can be also seen in the case of the Sun (Basri, 2018).

Figure 5 shows the results of light-curve-modeling spots in comparison with in-transit spots. The upper panel shows a comparison between Kepler light curve and the modeled one. The Kepler light curve is reproduced by our simple Gaussian model. That is because there is only hemispheric information in the light curve, which is straightforward to model (but not so accurate). The middle panel shows the estimated areal evolution of the light-curve-modeling spots. The estimated parameters are listed on Table 1. The bottom panel shows that the estimated location of the light-curve-modeling spots matches that of the in-transit spots, especially for the red and purple spots. This would be because the purple and red spots are located near the equator. On the other hand, green and blue spots are not rotating with the same rotational period as the equator, so there is no corresponding spot occulted by transit. If we compare Figure 3 and 5, we can see that the locations of the local minima and the estimated spot are quite similar to each other.

We also tried a five-spot model and no significant new information was obtained. Here, we note that long-term spot modeling, covering over a quarter ($\sim$ 90 days), should be done carefully because Kepler light curves have an inevitable long-term instrumental trend in the light curve, and the absolute values may be unreliable.

We estimated the location and area of star spots occulted by the planet for all Kepler short time cadence data (16 quarters in total). The estimated result is plotted in Figure 6. There are many long-lived recurrent spots that are located on the same longitude for a long time. To pick up candidates of long-lived spots, we apply DBSCAN, a commonly-used data clustering algorithm, in the Python package scipy (Martin, 1996; Davenport, 2015). In brief, this algorithm finds the core points which have more than $N$ satellite points within the length of $\epsilon$, and find clusters by connecting the core points with each other. The detected clusters are plotted with the colored symbols in Figure 6. We interpret these clusters as the long-lived spots and measure the temporal evolution of the spot area in each cluster. Most of the spots show very complex areal evolutions probably due to the consecutive flux emergence events in the same place, while some of them show clear emergence and decay phases as plotted in Figure 7.

We then compared the evolution of star spot areas detected by (a) local minima tracing method, (b) light curve modeling method, and (c) transit method in Figure 8. We focus on the spot group A in Figure 3 to see how the temporal evolution of star spot areas computed by the different methods compare. In panel (a) the spot group A indicated in Figure 3 is plotted. In panel (b) the corresponding spot evolution estimated by light-curve modeling is plotted on the basis of the location in Figure 8. As we expected, both of them show very similar temporal evolutions. This is very natural because both are obtained from the same data, but important for the validation of the local minima tracing method. In panel (c) we plot the temporal evolution of the selected spot area estimated by the transit method which is located on the longitude between 0^∘ and 100^∘ in the bottom panel of Figure 3. As one can see, the spot group A actually consists of at least four spots (two spots at the same time), and the temporal evolution of in-transit spots does not match with that of spot group A. The red circle shows the most dominant spot showing a clear emergence/decay phase, but the peak time is different from those seen in panel (a) and (b). This means that the temporal evolutions of the individual spots are different from those of spot groups (active longitudes), especially in this period.

We estimated the temporal evolution of individual star spot areas occulted by the planet and compared with those derived by other methods. We estimated the evolution of star spot area shown in Figure 7 with linear fitting (see red lines in Figure 7), and calculated the maximum area (flux), emergence rate, decay rates, and lifetime in the same manner as Namekata et al. (2019). We assume the mean magnetic field is about 2000 G when we compare them with the solar values. The estimated emergence rate is $1.1\times 10^{21}$ Mx$\cdot$h^-1 on average, and the decay rate is $-7.8\times 10^{20}$ Mx$\cdot$h^-1, for the spot area of $1.3\times 10^{24}$ Mx. On the other hand, the emergence/decay rates of spot group derived from local minima is $3.8\times 10^{20}$ Mx$\cdot$h^-1 and $-5.6\times 10^{10}$ Mx$\cdot$h^-1 for the spot area of $2.0\times 10^{24}$ Mx, respectively. Likewise, the emergence/decay rates of spot group derived from light curve modeling is $1.1\times 10^{20}$ Mx$\cdot$h^-1 and $-3.9\times 10^{20}$ Mx$\cdot$h^-1 for the spot area of $1.4\times 10^{24}$ Mx, respectively. The spot area is quite similar, but the flux emergence/decay rates derived from rotational modulations are smaller by one order of magnitude than those derived from the transit method. Here, the spot group A has an equatorial rotational period and expected to be on round the equator, so the latitudinal effects are negligible. Each data are plotted in Figure 9 in comparison with those of sunspots and star spot in our previous work. As you can see, however, the emergence/decay rates of star spot occulted by transit look consistent with sunspots and star spots in our previous work, although they scatter by an order of magnitude. The emergence rates can be explained by the solar scaling relation (d$\Phi$/d$t\propto\Phi^{0.3-0.5}$, Otsuji et al., 2011; Norton et al., 2017), and the decay rates can be also explained by the solar relation (d$\Phi$/d$t\propto\Phi^{\sim 0.5}$, Petrovay & van Driel-Gesztelyi, 1997).

We also compared the lifetime-area relation in Figure 10. As a result, the lifetimes of star spots occulted by transit are consistent with previous star spot studies (Bradshaw & Hartigan, 2014; Davenport, 2015; Giles et al., 2017; Namekata et al., 2019), and the lifetimes of star spot are much shorter than extrapolated from the solar empirical relations ($T\propto A$, so-called Gnevyshev-Waldmeier law; Gnevyshev, 1938; Waldmeier, 1955). On the other hand, the area-lifetime relation is roughly derived to be $T$ [d] $\sim$ C ($A$ [MSH])^0.5, where C is $\sim 0.8$, from the above dependence of emergence/decay rates on the total flux. This empirical relation would be more consistent with the stellar observations than the Gnevyshev-Waldmeier law as discussed in the following section.

In this section, we summarize and discuss the locations and spot areas on Kepler 17. In Sect. 4.1, we showed that one local minimum actually consists of several dominant spots in the case of Kepler 17, which has been already indicated by Lanza et al. (2019). This can be easily understood in analogy with sunspot distributions where we can see several active regions at the same time during the maximum activity. The locations of local minima sometimes nicely match those of (nearly equatorial) transited dominant spots; in those cases we can pin much of the light curve modulations to those spots (see, Figure 4).

We also showed the amplitudes of local minima values also have a positive correlation with the visible maximum spot area in transit and the total projected spot area in transit in this observational period. This positive correlation may imply that (1) the unseen spots are randomly distributed and so have a relatively small effect on the brightness variation and the dominant spots (group) in transit mainly create the rotational modulation (e.g. Eker, 1994) or (2) the unseen spots follow the same general hemispheric pattern as the transited spots. In the case of (2), the unseen spots are perhaps part of the same general active areas, but not always, because sometimes the areas trend away from each other.

We also showed that the amplitudes of local minima match the visible maximum spot area in transit, while they are smaller by a factor of 2 than the total projected spot area in transit. This factor difference would be caused because spots are widely distributed in longitude (see Figure 6), which decreases the brightness variation amplitude as well demonstrated by Eker (1994).

The consistency between the amplitudes of local minima and the visible maximum spot area in transit indicates that the spot area derived from local minima (or amplitude of the brightness variations) can be still a good indicator of that of the largest spot on the disk, but do not always correspond to the total filling factor. The similar properties can be seen in the case of the Sun as Basri (2018) showed.

In this section, we discuss the temporal evolution of star spot area in comparison with those of sunspots and our previous study (Namekata et al., 2019). In Figure 9, it is clear that the emergence/decay rates of spatially resolved star spot (resolution is still $\sim$ 10-20 deg) are consistent with those of sunspots within one order of magnitude of error bars of solar data. Also, this result is roughly consistent with our previous study within one order of magnitude. (red and blue circles, Namekata et al., 2019), which measured the variations rates of the favorable individual spots (i.e., isolated local-minima series showing clear emergence and decay). This possibly supports that the spot emergence/decay can be explained by the same mechanism, and imply a possibility to apply solar physics to star spot emergence/decay. However, the solar ($P_{\rm rot}\sim$ 25 d) and Kepler-17 ($P_{\rm rot}\sim$ 12 d) data is more consistent with the rapid rotators ($P_{\rm rot}<$ 7 d; blue circles) and larger than the slowly rotators ($P_{\rm rot}>$ 7 d; red circles). We speculate that the discrepancy between spatially resolved and non-resolved spots on the slowly rotating stars can be a result of the superposition effect of the several spot evolutions, as showed in Section 4.3 and discussed in the next paragraph. As for the rotational period dependence, Namekata et al. (2019) have discussed some possible mechanisms of the dependence of the rotational period on the spot evolution (e.g. decay due to the differential rotations of the stars), although it is still debated.

How about the comparison between the in-transit (spatially-resolved) and rotational-modulation (spatially-unresolved) spots? In the Kepler-17 system, it is reported that the two active longitudes are prominent, whose lifetimes are over the $Kepler$ observational period ($\sim$ 1400 days), while the brightness variation amplitude is varying every moment (Lanza et al., 2019). This is partly due to intermittent flux emergence and decay on the stellar surface. In fact, in Sect. 4.3, we also showed that the temporal evolution of individual star spot area in transit looks different from those derived from the local minima tracing method and light curve modeling method. These results would imply that the temporal evolution of star spot derived from out-of-transit light curves (i.e., rotational modulations) can be actually a superposition of the several dominant spots existing at the same active longitude. As a result of the superposition effect, as in Figure 9, the emergence/decay rates of the star spots in transits ($\sim\pm 1\times 10^{21}$ Mx$\cdot$h^-1) are larger by more than a factor of two (up to one order of magnitude) than those derived from out-of-transit light curves ($\sim\pm 1$-$6\times 10^{20}$ Mx$\cdot$h^-1), although the maximum spot area is quite similar ($\sim 10^{24}$ Mx). This is a feature discovered only on Kepler 17, but can be applicable to the other stars. So far, most of the star spot evolutions are estimated based on the $Kepler$/$COROT$ rotational modulations (Bradshaw & Hartigan, 2014; Giles et al., 2017; Namekata et al., 2019), our results propose that there is a certainty that the superposition effect changes the lifetime and variation rates by some factor, and the lifetimes may not be those of spot groups but those of active longitudes. Even though Namekata et al. (2019) carefully chose the favorable targets, we therefore have to be careful on the qualitative discussions.

In fact, most of the in-transit spots do not show clear emergence/decay phase. This may mean that there are continuous flux emergences and as a result they do not show clear emergence/decay phase, but can mean that the we may pick up only spots having rapid emergence/decay and the variation rates can actually have much larger diversity (more than one order of magnitude).

We also comment on the difference between the light curve modeling method and the local minima tracing method. In Sect. 4.3, the temporal evolution of star-spot area derived by light curve modeling are consistent with those derived by the local minima tracing method, and they have similar values of the variation rates ($\sim\pm 1$-$6\times 10^{20}$ Mx$\cdot$h^-1) and the maximum spot sizes ($\sim 10^{24}$ Mx). The reasons for the small difference in variation rates (a factor of $\sim$ 5) are as follows. First, the light curve modeling method improves the results of the local minima tracing method by considering the contamination of different spots. Second, the $Kepler$ light curves have inevitable long-term trends, so the relative values between different observational quarters are not reliable. In our light curve, the mean value is set to be zero for each quarter. Because of this, the light curve modeling method can generate a pseudo-long-lived spots, which result in the above difference (see, Basri, 2018, for more detais).

Finally, we comment on the star spot lifetime. Figure 10 shows that the lifetimes of detected star spots are smaller than those expected from solar empirical relation (T = A/D; D=10), and the result is consistent with the other studies. This is not surprising because the positive correlations of the variation rates and spot area (d$A$/d$t\propto A^{\alpha}$, $\alpha$=0.3-0.5) result in more small power-law relations $T\propto A^{1-\alpha}$ (the detail about the star spot lifetimes are discussed by Namekata et al. (2019)). Please note that Figure 6 that all of the individual star spots seem to appear for only several hundred days (see also Davenport, 2015; Lanza et al., 2019). This fact can be a strong restriction on the spot-evolution physics because our previous study Namekata et al. (2019) could not exclude the existence of individual spots with lifetimes of more than 1,000 days as predicted from the solar Gnevyshev-Waldmeier relation.

In the Kepler 17, we cannot find any stellar flares in the light curve, although the stellar brightness variations indicate the existence of large star spots that have a potential to produce a superflare ($>10^{34}$ erg). Also, Maehara et al. (2017) reported solar-type stars that do not show any superflares but have large star spots. The reason why such stars with large star spots do not show any superflares is an open question. Statistically, stars with such large star spots $>10^{-1.5}$ cause superflares 0.1-1 times per one year (Maehara et al., 2017), so such stars without superflares during $Kepler$’s observational period can be classified to have relatively less flare-productive spots. One possibility is that the large spots without any flares, like spots on Kepler 17, can have a simple polarity shape (e.g., $\alpha$-type or $\beta$-type spots), which is known to rarely produce extreme flares in the case of the sunspots (e.g., Toriumi & Wang, 2019). On sunspots, complex spots show relatively fast decay, so that lifetimes of star spots can be an indicator of spot complexity and flare productivity. The comparison of star-spot lifetimes between flare-productive and less flare-productive stars would be important for why and how stellar superflares occur.

Roettenbacher & Vida (2018) reported that the stellar superflares do not appear to be correlated with the rotational phase on solar-type stars. They tried to explain this result by the existence of large polar spots, visible large flares over the limb, or the flares between the active longitudes. Our transit model reveals that there are, more or less, visible spots on the disk regardless of its rotational phase, and the local maxima are not always the unspotted brightness level. These observations imply that, even if we assume that the superflares are accompanied with the solar-like low-latitude spots, superflares can apparently occur regardless of its rotational phase because large spots are always visible to some extent.

Finally, large star spots having a potential to produce superflares are found to survive more than 100 days (up to 1 year). This means that the surrounding exoplanets can be exposed to the threat of stellar superflares such a long period once large star spots appear, which can be critical for the exoplanet habitability (e.g., Segura et al., 2010; Airapetian et al., 2016; Yamashiki et al., 2019).