# Starspot mapping with adaptive parallel tempering. II. Application to TESS
data for M-dwarf flare stars, AU Microscopii, YZ Canis Minoris, and EV
Lacertae

Kai Ikuta Current Affiliation: Department of Multidisciplinary Sciences, The
University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan  Astronomical
Observatory, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo, Kyoto
606-8502, Japan  Kosuke Namekata ALMA Project, National Astronomical
Observatory of Japan, NINS, Mitaka, Tokyo 181-8588, Japan  Yuta Notsu
Department of Earth and Planetary Sciences, Tokyo Institute of Technology,
2-12-1 Ookayama, Meguro, Tokyo 152-8551, Japan  Laboratory for Atmospheric and
Space Physics, University of Colorado Boulder, 3665 Discovery Drive, Boulder,
Colorado 80303, USA  National Solar Observatory, 3665 Discovery Drive,
Boulder, Colorado 80303, USA  Hiroyuki Maehara Astronomical Observatory, Kyoto
University, Kitashirakawa-Oiwake-cho, Sakyo, Kyoto 606-8502, Japan  Okayama
Branch Office, Subaru Telescope, National Astronomical Observatory of Japan,
NINS, Kamogata, Asakuchi, Okayama 719-0232, Japan  Soshi Okamoto  Department
of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo, Kyoto
606-8502, Japan  Japan Meteorological Agency, 1-3-4 Otemachi, Chiyoda, Tokyo
100-8122, Japan  Satoshi Honda  Nishi-Harima Astronomical Observatory,
University of Hyogo, 407-2 Nishigaichi, Sayo-cho, Sayo, Hyogo 679-5313, Japan
Daisaku Nogami  Department of Astronomy, Kyoto University, Kitashirakawa-
Oiwake-cho, Sakyo, Kyoto 606-8502, Japan  Kazunari Shibata  Kwasan
Observatory, Kyoto University, 17 Ohmine-cho, Kita-Kazan, Yamashina, Kyoto
607-8471, Japan  Department of Environmental Systems Science, Doshisha
University, 1-3 Tataramiyakodani, Kyotanabe, Kyoto 610-0394, Japan

(Received July 28, 2022; Revised January 23, 2023; Accepted February 16, 2023)

###### Abstract

Starspots and stellar flares are indicators of stellar magnetic activity. The
magnetic energy stored around spots is thought to be the origin of flares, but
the connection is not completely understood. To investigate the relation
between spot locations deduced from the light curves and occurrence of flares
therein, we perform starspot modeling for TESS light curves of three M-dwarf
flare stars, AU Mic, YZ CMi, and EV Lac, using the code implemented in Paper
I. The code enables to deduce multiple stellar/spot parameters by the adaptive
parallel tempering algorithm efficiently. We found that flare occurrence
frequency is not necessarily correlated with the rotation phases of the light
curve for each star. The result of starspot modeling shows that either spot is
always visible to the line of sight in all phases, and we suggest that this
can be one of the reasons that there is no or less correlation between
rotation phases and flare frequency. In addition, the amplitude and shape of
the light curve for AU Mic and YZ CMi have varied in two years between
different TESS Cycles. The result of starspot modeling suggests that this can
be explained by the variations of spot size and latitude.

Starspots (1572), Stellar flares (1603), M dwarf stars (982), Markov chain
Monte Carlo (1889), Importance sampling (1892), Model selection (1912),
Astrostatistics (1882)

††journal: ApJ††software: stella (Feinstein et al., 2020), corner (Foreman-
Mackey, 2016)

##  1 Introduction

Starspots are apparent manifestations of magnetic activity on the stellar
surface, and can be ubiquitously observed on various type of stars (for
reviews, Berdyugina, 2005; Strassmeier, 2009; Kochukhov, 2021). For active
young stars, cool stars (M, K, G-dwarfs), and RS CVn-type stars, starspots
have been extensively studied through ground-based observations of photometry
(e.g., Henry et al., 1995; Messina & Guinan, 2002) and spectroscopy (e.g.,
Strassmeier & Bopp, 1992; Strassmeier et al., 1999). Stellar flares are
intense explosions in the stellar atmosphere by releasing magnetic energy
around starspots, and have extensively been studied in close association with
spots (e.g., Pettersen, 1980; Hawley & Pettersen, 1991; Gershberg, 2005; Benz
& Güdel, 2010). In particular, starspot and stellar flare studies have been
enabled with high precision and long-term photometry by the advent of Kepler
space telescope (Koch et al., 2010) and Transiting Exoplanet Survey Satellite
(TESS; Ricker et al., 2015). Many superflares are reported on a large number
of cool stars for Kepler data (Maehara et al., 2012; Hawley et al., 2014;
Davenport et al., 2014; Karoff et al., 2016; Notsu et al., 2019; Okamoto et
al., 2021) and TESS data (e.g., Günther et al., 2020; Tu et al., 2020;
Feinstein et al., 2020; Tu et al., 2021; Namekata et al., 2022a, b). The spot
size estimated from the amplitude of the light curve is determined by the
upper limit of the flare energy for Kepler data of solar-type stars, and stars
with larger spots cause larger flares as in the case of solar flares (Notsu et
al., 2013; Shibata et al., 2013; Maehara et al., 2017; Okamoto et al., 2021).
Hawley et al. (2014) and Silverberg et al. (2016) investigate the relation
between the rotation phase of Kepler light curve and occurrence of flares
therein for an active dMe star GJ1243 and found that the rotation phase and
flare frequency are not correlated. They suggest that this could be because
there are large spots on the visible pole or small spots distribute in the
whole stellar disk. The ensemble studies also suggest that the flare frequency
is not necessarily correlated with the phase of the light curve for Kepler/K2
data (e.g., Doyle et al., 2018) and TESS data (Doyle et al., 2020; Feinstein
et al., 2020).

It is unknown where the spots are actually on the stellar surface because the
surface cannot be spatially resolved with a simple analysis from the light
curve. Therefore, it is important to investigate the relation between the spot
locations and flare occurrence by mapping the surface with sophisticated
analysis, hereafter referred to as starspot modeling. Davenport et al. (2015)
conducted starspot modeling for the Kepler light curve of GJ1243 with two-spot
model and showed that the light curve can be explained by a spot on the
stellar equator and the other spot at the high latitude. Their starspot
modeling supports the suggestion by Hawley et al. (2014) that flare frequency
is not correlated with the spot location because of a always visible large
spot on the pole. This kind of sophisticated analysis can be important to
understand the origin of flares, but a small number of stars including GJ1243
have been investigated from this point. Then, it is also necessary for
modeling of more number of light curves to understand the relation between
spot locations and flare frequency on flare stars.

AU Microscopii (AU Mic), YZ Canis Minoris (YZ CMi), and EV Lacertae (EV Lac)
are bright and magnetically active M-dwarf stars. These have been observed for
monitoring their flares and to investigate their spots by (Zeeman) Doppler
imaging through ground-based observations. Recently, these targets have been
spectroscopically observed at the same time as the TESS observations to delve
into the flare mechanisms (e.g., Maehara et al., 2021). The spot locations can
be one of indicators where flares occur on the stellar disk and which
direction the plasma associated with flares erupts to (Namekata et al.,
2022a). AU Mic is a young dM1e dwarf with a spatially resolved debris disk
(Kalas et al., 2004), and was revealed to be a young planetary system
harboring two warm Neptunes (Plavchan et al., 2020; Martioli et al., 2020;
Gilbert et al., 2022) from the transit search (Hebb et al., 2007). Since AU
Mic is magnetically active and frequently causes large flares (Kunkel, 1970),
it has been spectroscopically observed from radio (Kundu et al., 1987) to
X-ray (Mitra-Kraev et al., 2005). Its surface is largely filled with strong
magnetic fields (Klein et al., 2021, 2022), and the light curve exhibits
rotational periodic modulations ascribed to starspots on the surface (Torres
et al., 1972; Rodono et al., 1986; Butler et al., 1987; Plavchan et al.,
2020). YZ CMi is a magnetically active dM4.5e dwarf. Flares are reported from
radio (Spangler et al., 1974) through optical (Lacy et al., 1976) to X-ray
(Grindlay & Heise, 1975), and has been spectroscopically observed for the
flare properties (Kowalski et al., 2010; Maehara et al., 2021). The modulation
of light curve is reported through ground-based observations (Pettersen et
al., 1983; Zboril, 2003). There are spots on the pole by Zeeman Doppler
imaging (Morin et al., 2008). EV Lac is a magnetically active dM3.5e dwarf and
has been investigated for their flare properties (Osten et al., 2005; Honda et
al., 2018; Muheki et al., 2020; Paudel et al., 2021). The modulation of light
curve is reported through ground-based observations (Pettersen, 1980;
Pettersen et al., 1983), and there are also spots on the pole by (Zeeman)
Doppler imaging (Morin et al., 2008; Jeffers et al., 2022).

In this study, for the purpose of investigating the spot properties and their
relation with flare occurrences, we conduct starspot modeling for TESS light
curves of M-dwarf flare stars, AU Mic, YZ CMi, and EV Lac, using a code
implemented in Ikuta et al. (2020) (hereafter, referred to as Paper I). In
Paper I, the code is implemented for starspot modeling with an adaptive
parallel tempering (PT) algorithm (e.g., Hukushima, 1999; Vousden et al.,
2016) and an importance sampling algorithm to deduce stellar/spot parameters
specified in macula (Kipping, 2012) and compare the number of spots in the
Bayesian framework. This paper is organized as follows. In Section 2, we
describe the TESS data, flare detection, calculation of flare energy, and the
numerical setup for starspot modeling of the data. In Section 3, we discuss
the results of starspot modeling in terms of comparing with previous studies
by photometry and spectroscopy, the relation between spot locations and flare
frequency, the variation of the amplitude and shape of the light curve for AU
Mic and YZ CMi in two years, and comparing the number of spots in the Bayesian
framework in addition to refering the validity of starspot modeling. In
Section 4, we conclude this paper and describe future prospects. In Appendix
A, we describe the result with different assumptions for the light curve of EV
Lac Cycle 2 to evaluate the uncertainties of parameters. In Appendix B, we
exhibit supplementary figures delineating the joint posterior distributions
for each of the model.

##  2 Method

###  2.1 Preprocessing TESS data and flare detection

AU Mic (TIC 441420236), YZ CMi (TIC 266744225), and EV Lac (TIC 154101678)
were observed with 2-min cadence by TESS in Sector 1 (July 25 to August 22 in
2018), 7 (January 7 to February 2 in 2019), and 16 (September 11 to October 7
in 2019) through the TESS prime mission (Cycle 1 and 2), respectively. After
that, AU Mic and YZ CMi were also observed in Sector 27 (July 4 to 30 in 2020)
and 34 (January 13 to February 9 in 2021) through the TESS extended mission
(Cycle 3), respectively. We retrieve the PDC-SAP (Pre-search Data Conditioning
Simple Aperture Photometry) flux from the MAST (Mikulski Archive for Space
Telescopes) Portal
webcite.111https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html To
confirm effects by the TESS systematics due to scattered light and spacecraft
drift, we derive the rotation period from the SAP light curves using the TESS
Systematics-insensitive Periodogram (Hedges et al., 2020). As a result, we
obtain the periods of $4.872$, $4.825$, $2.778$, $2.778$, and $4.364$ days,
for AU Mic Cycle 1 and 3, YZ CMi Cycle 1 and 3, and EV Lac Cycle 2,
respectively. These values are consistent with the results deduced from the
PDC-SAP light curves, and the detrended light curves correspond to that of the
PDC-SAP light curves. Thus, we adopt the PDC-SAP light curves for the starspot
modeling. For the case of EV Lac Cycle 2 as an example case, we exhibit the
result by the two-spot model with the systematic noise model by Gaussian
process (Appendix A.3).

Figure 1 exhibits TESS light curves in terms of the relative flux (black), the
PDC-SAP flux normalized by its average excluding missing values. Vertical
lines in each figure also represent the detected flares at each time (red) as
following procedures. First, we detect flare candidates in the light curves by
using stella (Feinstein et al., 2020). The code detects flares from TESS light
curves with convolutional neural networks, and we adopt their pretrained model
with detected flares in Günther et al. (2020). Then, we adopt the flare
probability of 0.1 and an amplitude of $10^{-3}$ for all trained hundred
models as fewer false negatives as possible. We confirm the flare candidates
as flares by eye because many false positives are included in the detected
candidates (Feinstein et al., 2022). Second, after extracting confirmed flares
from the light curves, we also extract outliers of the modulations ascribed to
spots from the remaining light curves with Bessel filter (Maehara et al.,
2021). Although a small number of outliers remain in the light curves, they do
not affect the result of the starspot modeling since the number of data points
is much less than that of total data points. In this study, since the
modulation timescale ascribed to spots is much longer than the cadence, we use
data every three points for computational efficiency. As a result of the flare
detection, the number of detected flares are summed up to 55, 70, 124, 137,
and 76, for AU Mic Cycle 1 and 3, YZ CMi Cycle 1 and 3, and EV Lac Cycle 2,
respectively (Figure 1).222Each flare properties is listed on
https://github.com/KaiIkuta/StarspotMapping.git: the peak time (BJD-2745000),
amplitude, equivalent duration (sec), and e-folding time (day). The
detectability of flares may depends on the rotation phase of the light curve
because the amplitude of the light curve for these stars are large, and flares
in the training set of the pretrained model are detected by reducing the
periodic signal of the light curve (Günther et al., 2020). Thus, we consider
the equivalent duration = 1 (sec) as the sufficient threshold to discuss the
correlation between spot locations and flare frequency. Among detected flares,
the number of flares with larger equivalent duration than 1 (sec) are 48, 50,
109, 106, and 71, respectively.

According to Shibayama et al. (2013), a bolometric flare energy is formulated
as

$\displaystyle E_{\rm flare}$ $\displaystyle=\int L_{\rm
flare}(t)dt=\sigma_{\rm SB}T_{\rm flare}^{4}\times 2\pi R_{\rm
star}^{2}\times\frac{\int R(\lambda)B(\lambda,T_{\rm eff})d\lambda}{\int
R(\lambda)B(\lambda,T_{\rm flare})d\lambda}\times\int\frac{\Delta F}{F}(t)dt,$
(1)

where $L_{\rm flare}(t)$ and $\Delta F/F(t)$ are the bolometric luminosity and
relative flux at the time $t$, and $\sigma_{\rm SB}=5.67\times 10^{-1}$ (erg
K-4 m-2 sec-1) is the Stefan–Boltzmann constant. The flare is assumed to be
the blackbody radiation with the flare temperature $T_{\rm flare}=10000$ (K).
The last time-integrated factor corresponds to the equivalent duration, which
is calculated by the integral of the relative flux with time from the start to
end of a flare (Hawley et al., 2014). The coefficient specific to a star is
the flare energy with the equivalent duration $=1$ (sec). Those for AU Mic, YZ
CMi, and EV Lac are calculated to be $3.68\times 10^{32}$, $3.56\times
10^{31}$, $5.25\times 10^{31}$ (erg), respectively. We note that the
calculation error up to a factor is inherent in the assumed flare model
(Maehara et al., 2021). If assuming the flare temperature $T_{\rm flare}=7000$
(K), the energy reduces by a factor of two. Figure 2, 3, and 4 accumulates
flares as their flare frequency for AU Mic, YZ CMi, and EV Lac, respectively
(cf. Maehara et al., 2021). We discuss the relation between spot locations and
flare frequency in Section 3.2 and 3.3.

###  2.2 Numerical setup

Table 2.2 lists stellar parameters, references, deduced spot temperature and
intensity, and limb-darkening coefficients for AU Mic, YZ CMi, and EV Lac. The
spot temperature $T_{\rm spot}$ is deduced from the quadratic formula in the
stellar effective temperature $T_{\rm eff}$ (Equation 4 in Herbst et al.,
2021):

$T_{\rm spot}=-3.58\times 10^{-5}T_{\rm eff}^{2}+0.801T_{\rm eff}+666.5.$ (2)

Thereby, the spot relative intensity $f_{\rm spot}$ for TESS band is
calculated by the stellar effective temperature $T_{\rm eff}$ and the spot
temperature $T_{\rm spot}$:

$f_{\rm spot}=\frac{\int R(\lambda)B(\lambda,T_{\rm spot})d\lambda}{\int
R(\lambda)B(\lambda,T_{\rm eff})d\lambda},$ (3)

where $R(\lambda)$ and $B(\lambda,T)$ are the TESS-band response function and
the Plank function with each wavelength $\lambda$, respectively. The stellar
limb-darkening law is adopted as the following formula for the cosine of the
angle to the line of sight $\mu$:

$I(\mu)/I(1)=1-\Sigma_{k=1}^{4}c_{k}(1-\mu^{k/2}),$ (4)

and the coefficients $c_{k}$ based on the stellar parameters are deduced from
Claret (2018). Then, we adopt the nonlinear limb-darkening law. For the case
of EV Lac Cycle 2 as an example case, we exhibit the result by the quadratic
limb-darkening law (Appendix A.1).

In Paper I, we implement a code for starspot modeling using an adaptive PT
algorithm. The model evidence for each model is also computed with an
importance sampling algorithm together with the PT transition. We adopt an
analytical spotted model (macula; Kipping, 2012) for the spot model, and the
model is specified in the stellar and spot parameters: the stellar inclination
angle $i$ (deg); equatorial rotation period $P_{\rm eq}$ (day); degree of
differential rotation $\kappa$; spot relative intensity $f_{\rm spot}$;
latitude $\Phi$ (deg); initial longitude $\Lambda$ (deg); reference time
$t_{\rm ref}$ (day) (the time at the midpoint of the interval over which the
spot has its maximum radius); maximum radius $\alpha_{\rm max}$ (deg);
emergence duration ${\cal I}$ (day); decay duration ${\cal E}$ (day); stable
duration ${\cal L}$ (day). The rotation period at the latitude $\Phi$
(Equation 16 in Paper I) is described by

$P(\Phi)=\frac{P_{\rm eq}}{1-\kappa\sin^{2}\Phi}.$ (5)

In this study, it is unnecessary to adopt the spot parameters of the reference
time $t_{p}$, emergence duration ${\cal I}$, decay duration ${\cal E}$, and
stable duration ${\cal L}$ in the model because the amplitudes of the light
curves are approximately constant in the observation period of a TESS Sector.
Therefore, the light curves are specified with fewer stellar and spot
parameters than those used in Paper I for simplicity (Table 2, 3, 4, 5, and
6). In fact, we ascertain that the posterior distribution does not converge to
an unimodal distribution due to these spot parameters by exploiting the test
runs. Each of the spot is discerned by the range of latitude $\Phi$ from the
southern to northern pole, unlike the range of the reference time $t_{p}$ in
Paper I.

Light curves of AU Mic Cycle 1 and 3, YZ CMi Cycle 3, and EV Lac exhibit two
local minimum per one equatorial rotation period, and we optimize them by two-
and three-spot models (Table 2, 3, 5, and 6). The light curve of YZ CMi Cycle
1 exhibits one local minimum per one equatorial rotation period, and we
optimize it by two- and one-spot models (Table 4). As for the one-spot model,
the differential rotation parameter $\kappa$ is fixed to 0. We note that the
light curves of YZ CMi have one local minimum per one equatorial rotation
period possibly due to its low inclination angle (e.g., Basri & Shah, 2020).
We ascertained the posterior distributions of the inclination angle and spot
intensity would converge to unphysical modes (the edges of the parameter
range) by exploiting test runs because they have respective degeneracies with
the spot latitude and radius (Paper I). Therefore, the values of the
inclination angle and spot intensity are fixed to the values in Table 2.2, and
we abbreviate the name followed by Paper I for each of the model: for
instance, “two-spot model fixed $\sin i$ and $f_{\rm spot}$” in the context of
Paper I is simply renamed “two-spot model” in this study. For the case of EV
Lac Cycle 2 as an example case, we exhibit the result that the spot intensity
is included in the parameters under an uniform and normal prior distributions
(Appendix A.2).

Table 1: Stellar parameters

Stellar parameter | AU Mic | YZ CMi | EV Lac
---|---|---|---
Effective temperature $T_{{\rm eff}}$ (K) |  $3700\pm 100$ c |  $3100\pm 50$ e |  $3400\pm 18$ h
Rotation period $P_{{\rm rot}}$ (day) |  $4.863\pm 0.010$ c |  $2.776\pm 0.010$ e |  $4.359$ j
Stellar radius $R_{\rm star}$ ($R_{\rm sun}$) |  $0.75\pm 0.03$ c |  $0.37^{+0.03}_{-0.06}$ e |  $0.35\pm 0.02$ j
Surface gravity $\log g$aa The surface gravity is calculated from the stellar radius $R_{\rm star}$ and mass $M_{\rm star}$. The spot temperature is formulated by the stellar effective temperature $T_{{\rm eff}}$ in Equation 2, and the spot relative intensity $f_{\rm spot}$ for TESS band is deduced from Equation 3.  | $4.39\pm 0.10$ | $4.87\pm 0.33$ | $4.89\pm 0.13$
Inclination angle $i$ (deg) |  $75$ d |  $36^{+17}_{-14}$ e |  $60$ k
Spot temperature $T_{\rm spot}$ (K)aa The surface gravity is calculated from the stellar radius $R_{\rm star}$ and mass $M_{\rm star}$. The spot temperature is formulated by the stellar effective temperature $T_{{\rm eff}}$ in Equation 2, and the spot relative intensity $f_{\rm spot}$ for TESS band is deduced from Equation 3.  | $3140\pm 54$ | $2806\pm 29$ | $2976\pm 10$
Spot relative intensity $f_{\rm spot}$aa The surface gravity is calculated from the stellar radius $R_{\rm star}$ and mass $M_{\rm star}$. The spot temperature is formulated by the stellar effective temperature $T_{{\rm eff}}$ in Equation 2, and the spot relative intensity $f_{\rm spot}$ for TESS band is deduced from Equation 3.  | $0.43\pm 0.04$ | $0.56\pm 0.03$ | $0.48\pm 0.01$
Limb-darkening coefficients ($c_{1}$, $c_{2}$, $c_{3}$, $c_{4}$)bb The stellar limb-darkening coefficient for TESS data is characterized by the effective temperature $T_{\rm eff}$ and surface gravity $\log g$ under the solar metalicity (Claret, 2018). We adopt the coefficient values of ($T_{\rm eff}$, $\log g$)=($3700$, $4.5$), ($3100$, $5.0$), and ($3400$, $5.0$), for AU Mic, YZ CMi, and EV Lac, respectively.  | ($2.87$, $-4.35$, $3.92$, $-1.32$) | ($3.15$, $-4.62$, $4.01$, $-1.34$) | ($3.00$,$-4.54$, $4.01$, $-1.35$)

c Plavchan et al. (2020);d Wisniewski et al. (2019); e Baroch et al. (2020); f
Cifuentes et al. (2020); h Rojas-Ayala et al. (2012); j Paudel et al. (2021);
k Morin et al. (2008);

References. —

##  3 Results and discussion

We optimize the light curves of AU Mic Cycle 1 and 3, YZ CMi Cycle 3, and EV
Lac Cycle 2, by the two- and three-spot models, and the light curve of YZ CMi
Cycle 1 by the two- and one-spot models. The number of spots is set as the
number of local minimum in the light curve and one more number (Section 2.2).
For every cases, unimodal posterior distributions are deduced regardless of
the number of spots, and joint ones are delineated in Appendix B. Table 2, 3,
4, 5, and 6 show the modes of the deduced posterior distributions, their
credible regions, reproduced average flux, and the logarithm of model evidence
for each model, together with their prior distributions for each of the
parameters, for AU Mic Cycle 1 and 3, YZ CMi Cycle 1 and 3, and EV Lac Cycle
2, respectively. The preferable models and their values of the logarithm of
model evidence are marked in bold. Figure 1 exhibits the reproduced light
curves with each mode of the posterior distributions in each of the two-spot
model (red). All light curves can be folded in phase with the rotation period
since the deduced degree of differential rotation $\kappa$ is approximately 0
or the absolute value of spot latitudes is almost the same value, and the
amplitude of the light curve is approximately constant in the observation
period of a TESS Sector. Figure 5, 6, 7, 8, and 9 exhibit visualized surface
reproduced with each mode of posterior distributions; (a, f) phase-folded
light curve (black) and reproduced one with the model (red); (b, g) the
temporal variation of visible projected area of each spot (red, blue, and
green) and the total (black) relative to the disk; (c, h) that of visible area
relative to the photosphere; (d) flare frequency in each bin for all flares
(red) and flares with larger equivalent duration than 1 (sec) (blue); and (e)
flare energy in each bin for those. We discuss the result in comparison with
other studies by photometry and spectroscopy in Section 3.1, correlation
between spot locations and flare frequency in Section 3.2, two-year variation
of the light curves for AU Mic and YZ CMi in Section 3.3, and model selection
and the validity for starspot modeling in Section 3.4.

###  3.1  Deduced parameters in comparison with other studies by photometry
and spectroscopy

AU Mic

In the two-spot models (Table 2 and 3), the separation of spot longitude
extends from $129.54\pm 0.04$ (deg) for Cycle 1 to $138.93\pm 0.05$ (deg) for
Cycle 3. This slight change in the separation can be attributed to the
differential rotation if assuming that the same spots have existed for two
years and the star has significant differential rotation. In contrast, the
deduced degree of differential rotation $\kappa$ for Cycle 1 and 3 are
approximately 0, respectively. This value of the differential rotation can not
cause the change in the separation of spots. This may mean that the change in
the relative longitude of two spots can be caused for different reasons rather
than the differential rotation, although we need to pay attention to this
small value because the observation period of a TESS Sector ($\sim 27$ day) is
not enough to deduce the differential rotation accurately.

Wisniewski et al. (2019) also conducted starspot modeling for the same light
curve with two spots using STSP code (Davenport et al., 2015). Their result of
deduced spot latitudes $\simeq$ $9.6$ and $44.8$ (deg) are different from our
results $=-13.5^{+0.1}_{-0.2}$ and $39.7^{+0.1}_{-0.1}$ (deg) probably due to
the following three reasons: (i) the high inclination angle $i=75$ (deg) can
result in large uncertainty in deducing spot latitudes, (ii) we sets the spot
relative intensity $f_{\rm spot}=0.43$ (Table 2.2) by calculating the
intensity in TESS band based on the spot temperature deduced from the formula
(Equation 2 and 3), whereas Wisniewski et al. (2019) assumes 0.7 (Davenport et
al., 2015), which results in the difference due to degeneracies between the
spot relative intensity and latitudes ( Paper I, ), and (iii) it is also
difficult to discern whether a spot is on northern or southern hemisphere as
the inclination angle increases (Walkowicz et al., 2013). In addition, the
difference of deduced spot longitude in our result $=129.54\pm 0.04$ (deg) is
consistent with their result that the spots are separated by 131 (deg). This
consistency in the longitude is simply understood since the local minima of
the light curve only correspond to the separation of spot longitude in the
case of the two-spot model.

The surface map has been investigated through Doppler imaging and Zeeman
Doppler imaging between Cycle 1 and 3 (September to November 2019) (Klein et
al., 2021) and after Cycle 3 (November 2020 to September 2021) (Klein et al.,
2022). As a result of Doppler imaging, Klein et al. (2021) reported there were
two dark features from the equator and 40 (deg) in addition to bright one
around the equator in opposite to the dark one in 2019. Klein et al. (2022)
also reported there were also dark features on the equator and 60 (deg) in
2020. This results are different from our results of deduced latitudes for
Cycle 1 and 3 in the two-spot model. It should be noted that the result of
Wisniewski et al. (2019) is almost consistent with the results of Doppler
imaging. One reason for this difference between our result and their result by
Doppler imaging could be that there are degeneracies between the inclination
angle and latitude, especially the high inclination angle $i=75$ (deg)
(Walkowicz et al., 2013), in the starspot modeling. We need to be careful to
interpret these results obtained by both starspot modeling and Doppler imaging
because these results largely depend on the mapping method and assumptions,
such as prior distribution (Roettenbacher et al., 2017; Luger et al., 2021c).

In the three-spot models (Table 2 and 3), the latitudes in each model are
almost the same absolute values. Rotation period of each spot is almost the
same value because of the almost the same latitude (Equation 5), and therefore
the degree of differential rotation is not able to be measured accurately,
although large absolute values are deduced for Cycle 1 and 3.

YZ CMi

In the two-spot models (Table 4 and 5), the separation of spot longitude
extends from $96.32\pm 0.28$ (deg) for Cycle 1 to $193.95\pm 0.53$ (deg) for
Cycle 3. As the same as AU Mic, this slight change in the separation can be
attributed to the differential rotation if assuming that the same spots have
existed for two years and the star has significant differential rotation. In
contrast, the deduced degree of differential rotation $\kappa$ are
approximately 0, respectively. This may mean that the change in the relative
longitude of two spots can be caused for different reasons rather than the
differential rotation, although we need to pay attention to this small value
because the observation period of a TESS Sector ($\sim 27$ day) is not enough
to deduce the differential rotation accurately, as discussed for AU Mic.

Bicz et al. (2022) also conducted starspot modeling for the same light curves
with three and four spots for Cycle 1 and 3 using BASSMAN code, respectively.
Their model is different in the point that they assume the inclination angle
$=60$ (deg) and deduce the spot temperature (i.e, relative intensity) for each
spot, although we set the angle $i=36$ (deg) and the spot temperature $T_{\rm
spot}=2806$ (K) (relative intensity $f_{\rm spot}=0.66$) for all spots. We do
not respectively adopt three- and four-spot models for Cycle 1 and 3, and our
result of deduced spot parameter is very different from their result. The
difference is due to the models and assumptions, such as the number of spots
(see Section 3.4).

The surface map has been investigated through Doppler imaging before Cycle 1
(September 2016 to May 2017) (Baroch et al., 2020). The result suggests there
are spots on the latitude from 75 to 81 (deg), and this result is consistent
with our result of deduced latitudes $=79.0^{+0.01}_{-0.02}$ (deg) for Cycle 1
in the one-spot model.

EV Lac

The surface map has been investigated through Doppler imaging before Cycle 2
(January 2016 to December 2017) (Jeffers et al., 2022). The result suggests
that there are spots near the north pole, and is different from our result of
deduced latitudes in both two- and three-spot models (Table 6).

###  3.2 Spot location vs flare frequency

We investigate the correlation between spot locations and flare frequency in
each rotation phase of the light curve. In Figure 5, 6, 7, 8, and 9, for each
flares in the phase-folded light curve, we assume a null hypothesis that the
flares are evenly distributed across ten bins of the rotation phase and
perform a chi-squared test with nine degrees of freedom. For AU Mic Cycle 1
and 3, YZ CMi Cycle 1 and 3, EV Lac Cycle 2, we obtain a reduced chi-square
$=$ 0.899 (0.824), 0.381 (0.533), 0.505 (0.376), 0.439 (0.444), and 0.825
(0.890), and $p$-value $=$ 0.525 (0.594), 0.945 (0.851), 0.872 (0.947), 0.915
(0.911), and, 0.593 (0.533), for all flares (flares with larger equivalent
duration than 1 sec), respectively. Thus, the hypothesis that flares are
evenly distributed across ten bins is not rejected, and we can not conclude
that there is a correlation between the rotation phase of the light curve and
occurrence of flares. This result is consistent with the result in previous
studies for different stars (Hawley et al., 2014; Doyle et al., 2018, 2020).

In contrast, Bicz et al. (2022) suggested that the hypothesis that flares are
evenly distributed in the rotation phase can not be rejected for 80 flares in
Cycle 1 and can be rejected for 80 flares in Cycle 3. Our result is different
from their result because more number of flares are detected for 124 in Cycle
1 and 137 in Cycle 3, respectively (Section 2.1).

As a result of the starspot modeling, we can interpret the correlation between
rotation phases and flare frequency as follows. Our result shows that either
spot is always visible in all phases (Figure 5, 6, 7, 8, and 9) due to
multiple spots or moderate inclination angle for each star. This suggests that
flare frequency can be less correlated with the rotation phase if assuming
that flares occur around spot visible to the line of sight at that time.
Therefore, a polar spot is not necessarily needed to explain that the flares
distribute in all rotation phases, and no or less correlation between rotation
phases and flare frequency do not mean no relation between spot locations and
flare occurrence.

###  3.3 Two-year variation of light curves

AU Mic and YZ CMi were observed in TESS Cycle 1 and 3 with the interval of
approximately two years. It is shown that the amplitude of both light curves
decreases (Figure 1), and the flare frequency is almost constant within an
error (Figure 2 and 3)333We note that, when the flare frequency is compared,
the flare frequency is need be adopted rather than the cumulative one because
the cumulative one highly depends on large flares. In particular, the light
curve of YZ CMi changes from one to two local minima per one equatorial
rotation period. Therefore, we compare the starspot modeling of Cycle 1 with
that of Cycle 3 for AU Mic and YZ CMi in terms of spot latitude and size to
explain the variations of light curves ascribed to spots in two years between
different TESS Cycles. Then, we focus on the result of the two-spot model
because two-spot models are adopted for Cycle 1 and 3.

Comparing each of the two-spot models for AU Mic (Table 2 and 3), deduced
latitudes and radii are quite different. Although the latitude and radius of
the first spot near the equator respectively correspond within an error of 14
and 3 (deg), the latitude of the second spot changes from mid-latitude to the
pole. The radii of the second spot on the pole is not determined due to the
unspotted level of the light curve. We describe the unspotted level on the
starspot modeling below.

Comparing each of the two-spot models for YZ CMi (Table 4 and 5), deduced
latitudes and radii are slightly different. The latitude of the first spot
near the equator is different within an error of 25 (deg), and that of the
second one at mid-latitude within an error of 5 (deg) because the latitude is
determined by the visible duration of the spot (Walkowicz et al., 2013). It is
possible that the second spot decayed and the relative longitude changed due
to the stellar differential rotation or the spot in Cycle 1 decayed and a new
spot emerged in the interval of Cycle 1 and 3.

The total visible area relative to the stellar hemisphere ($2\pi R^{2}_{\rm
star}$) (Figure 5, 6, 7, and 8) changed from 0.076 to 0.142 for AU Mic and
from 0.063 to 0.035 for YZ CMi. We found that the flare frequency did not
change from Cycle 1 to Cycle 3 for both AU Mic and YZ CMi (Figure 2 and 3). We
suggest that this may be because the spot area deduced from the starspot
modeling changes by a factor of two in two years, but this could not affect to
change the flare frequency within the error.

We note that it is difficult to determine the area of always visible spot on
high latitude and the existence of a polar spot because such spot reduce the
average flux and does not significantly contributes the modulation of the
light curve. Thus, the unspotted level is not well determined in our result.
The unspotted level of flux is of importance to estimate spot area and other
spot parameters in the starspot modeling. For each of the two-spot (three- or
one-spot) models (Table 2, 3, 4, 5, and 6), the reproduced average flux
$F_{\rm ave}$ in the model are 0.831 (0.830), 0.816 (0.823), 0.839 (0.778),
0.847 (0.848), and 0.850 (0.842), respectively. Then, the fluxes of Cycle 3
relative to those of 1 for AU Mic and YZ CMi are 0.982 and 1.010,
respectively. This reflects that the spot area has become smaller or a always
visible polar spot that does not contribute to the brightness change has
appeared, and the amplitude of the light curve has become smaller from Cycle 1
to 3. The reason that the radius of the polar spots reach the upper limit of
the parameter range is that the unspotted level is not determined, and the
same light curve can be reproduced by reducing the average flux. Therefore, it
can be validated that there is only a polar spot in the two- and three-spot
models for AU Mic Cycle 3. The possible large spot on the high latitude could
affect the above conclusion (i.e., the change of the total spot area in two
years), and more careful analysis may be necessary.

###  3.4 Model selection and the validity for starspot modeling

More spots are indicated to be on the surface than the number of local minimum
per one equatorial rotation period by analyses of spot occultations (Morris et
al., 2017; Namekata et al., 2020a). In other words, the light curve of spotted
stars always have two or one local minimum per one equatorial rotation period
even if produced with any number of spots, but not two or one spot is actually
present on the surface (Basri & Nguyen, 2018; Basri, 2018; Luger et al.,
2021a, b). Then, as well as Paper I, we compare the number of spots based on
the model selection by comparing the model evidence (Kass & Raftery, 1995).
The values of the logarithm of model evidence $\log{\cal Z}$ for AU Mic Cycle
1 and 3, YZ CMi Cycle 1 and 3, and EV Lac Cycle 2 are listed in Table 2, 3, 4,
5, and 6, respectively.

In cases of AU Mic Cycle 1 and 3, YZ CMi Cycle 3, and EV Lac Cycle 2, the
three-spot model is much more decisive than the two-spot model by orders of
magnitude: the evidences of the three-spot model relative to that of the two-
spot model are $\Delta\log{\cal Z}=22568.926$, $53509.442$, $380.686$, and
$3162.217$, respectively. In the case of YZ CMi Cycle 1, the two-spot model is
much more decisive than the one-spot model by orders of magnitude: the
evidence of the two-spot model relative to that of the one-spot model is
$\Delta\log{\cal Z}=10356.959$. In all cases, it is suggested that each of the
model with one more additional spot is preferable in the Bayesian framework.
While we can also conduct starspot modeling with an additional spot, models
with less than three spots are sufficient to investigate the correlation
between the spot location and flare frequency (Section 3.2). We note that the
model evidence is affected by the choice of prior distributions (Gelman et
al., 1995). To constrain the number of spots precisely, it is also essential
to scan the surface with spot occultations by multiple transits of exoplanets
(Morris et al., 2017; Namekata et al., 2020a; Bruno & Deleuil, 2021).

The peak-to-trough amplitude of light curves for AU Mic 1 and 3, YZ CMi 1 and
3, and EV Lac Cycle 2 are 0.046 (0.049), 0.032 (0.035), 0.035 (0.036), 0.014
(0.015), and 0.027 (0.028), for each of the two-spot model (three- or one-spot
models), respectively. According to Figure 4 in Rackham et al. (2018), the
peak-to-trough amplitude can be only explained by giant spots on the surface
when the filling factors are approximately 0.22, 0.10, 0.12, 0.02, and 0.08,
respectively. Since large spot groups are necessary to generate the large
amplitude of light curves, it is sufficient to assume a few large spots for
M-dwarf stars.

In addition, Roettenbacher et al. (2017) showed that the comparison between
light curve inversion and Doppler imaging for a RS CVn star $\sigma$ Geminorum
reveals qualitatively different surface map is obtained due to the various
degeneracies inherent in each method. Therefore, as future studies, multi-
color and multi-spectrum mappings of the stellar surface are important to
constrain the surface map more robustly.

##  4 conclusion and future prospects

We conducted starspot modeling for TESS light curves of three M-dwarf flare
stars, AU Mic, YZ CMi, and EV Lac, using the code implemented in Paper I, for
the purpose of investigating the spot properties and their relation with flare
occurrences. The code enables to deduce multiple stellar/spot parameters by
the adaptive parallel tempering algorithm efficiently. The light curves of AU
Mic Cycle 1 and 3, YZ CMi Cycle 3, and EV Lac Cycle 2 are optimized by the
two- and three-spot models, and that of YZ CMi Cycle 1 by the two- and one-
spot models. For each model, unimodal posterior distributions are deduced
(Table 2, 3, 4, 5, and 6) regardless of the number of spots, and joint ones
are delineated in Appendix B. Flares are also detected (Figure 1) and
accumulated as the flare frequency (Figure 2, 3, and 4). We discuss the
results in terms of comparing with previous studies by photometry and
spectroscopy, the relation between spot locations and flare frequency, the
variation of the amplitude and shape of the light curve structures for AU Mic
and YZ CMi in two years, the model selection for each numbers of model, and
the validity of starspot modeling.

In Section 3, we described the main results summarized as follows:

  * (i)

The flare frequency is not necessarily correlated with rotation phase of the
light curve for each star. We suggest that this can be because either spot is
always visible in all phases for any model due to multiple spots or moderate
inclination angle for each star. As the result, flares can be distributed in
all rotation phases.

  * (ii)

For AU Mic and YZ CMi, the amplitude and shape of the light curve have varied
in two years, whereas the flare frequency is almost constant within an error.
The result of starspot modeling suggests that the variation of the light curve
can be explained by the variations of spot size and latitude. We also found
that the total spot area does not change significantly in two years enough to
change the flare frequency within the error, and this is why the flare
frequency did not vary in two years.

In particular, EV Lac is to be observed in Sector 56 and 57 (September 1 to
October 29 in 2022) through TESS second extended mission (Cycle 5). Then,
three-year variation of the light curve can be investigated as in Section 3.3.
In future studies, it is necessary to conduct surface mapping in multi-band
wavelengths (e.g., Morris et al., 2018; Namekata et al., 2020b; Toriumi et
al., 2020), taking into account the effect of bright spots (faculae) and
constraint on the spot temperature (e.g., Johnson et al., 2021). It is also
helpful to scan the stellar surface with spot occultations by multiple
transits of exoplanets (e.g., Morris et al., 2017; Namekata et al., 2020a;
Bruno & Deleuil, 2021). Surface mapping on the stellar surface plays one of
important roles in quantifying the stellar effects on exoplanet
characterizations (Rackham et al., 2018, 2022).

We sincerely appreciate the referee for providing feedback that helped to
improve both content of this manuscript and the clarity. This study is based
on public available data obtained by TESS mission. Funding for the TESS
mission is provided by NASA Science Mission Directorate. K.I. is grateful to
Taichi Kato and Naoto Kojiguchi for their instructions of data analysis.
Numerical computations were carried out on Cray XC40, Yukawa-21 (the Yukawa
Institute for Theoretical Physics, Kyoto University), and Cray XC50 (Center
for Computational Astrophysics, National Astronomical Observatory of Japan).
Our study was also supported by JSPS KAKENHI Grant Numbers: JP21H01131(all
authors), JP18J20048, JP21J00316(K.N.), JP21J00106(Y.N.), JP20K04032,
JP20H05643, and JP20H00173(H.M.). K.I. is supported by JST CREST Grant Number
JPMJCR1761. Y.N. was supported by JSPS Overseas Research Fellowship Program
and NASA ADAP award program number 80NSSC21K0632 (PI: Adam Kowalski). We also
acknowledges the International Space Science Institute and the supported
International Team 510: ”Solar Extreme Events: Setting Up a Paradigm”.

## References

  * Basri & Nguyen (2018) Basri, G. & Nguyen, H. T. 2018, ApJ, 863, 190. doi:10.3847/1538-4357/aad3b6
  * Basri (2018) Basri, G. 2018, ApJ, 865, 142. doi:10.3847/1538-4357/aade45
  * Basri & Shah (2020) Basri, G. & Shah, R. 2020, ApJ, 901, 14. doi:10.3847/1538-4357/abae5d
  * Baroch et al. (2020) Baroch, D., Morales, J. C., Ribas, I., et al. 2020, A&A, 641, A69. doi:10.1051/0004-6361/202038213
  * Benz & Güdel (2010) Benz, A. O. & Güdel, M. 2010, ARA&A, 48, 241. doi:10.1146/annurev-astro-082708-101757
  * Berdyugina (2005) Berdyugina, S. V. 2005, Living Reviews in Solar Physics, 2, 8. doi:10.12942/lrsp-2005-8
  * Bicz et al. (2022) Bicz, K., Falewicz, R., Pietras, M., et al. 2022, ApJ, 935, 102. doi:10.3847/1538-4357/ac7ab3
  * Bruno & Deleuil (2021) Bruno, G. & Deleuil, M. 2021, arXiv:2104.06173
  * Butler et al. (1987) Butler, C. J., Doyle, J. G., Andrews, A. D., et al. 1987, A&A, 174, 139
  * Cifuentes et al. (2020) Cifuentes, C., Caballero, J. A., Cortés-Contreras, M., et al. 2020, A&A, 642, A115. doi:10.1051/0004-6361/202038295
  * Claret (2018) Claret, A. 2018, A&A, 618, A20. doi:10.1051/0004-6361/201833060
  * Davenport et al. (2014) Davenport, J. R. A., Hawley, S. L., Hebb, L., et al. 2014, ApJ, 797, 122. doi:10.1088/0004-637X/797/2/122
  * Davenport et al. (2015) Davenport, J. R. A., Hebb, L., & Hawley, S. L. 2015, ApJ, 806, 212. doi:10.1088/0004-637X/806/2/212
  * Doyle et al. (2018) Doyle, L., Ramsay, G., Doyle, J. G., et al. 2018, MNRAS, 480, 2153. doi:10.1093/mnras/sty1963
  * Doyle et al. (2020) Doyle, L., Ramsay, G., & Doyle, J. G. 2020, MNRAS, 494, 3596. doi:10.1093/mnras/staa923
  * Feinstein et al. (2020) Feinstein, A., Montet, B., & Ansdell, M. 2020, The Journal of Open Source Software, 5, 2347. doi:10.21105/joss.02347
  * Feinstein et al. (2020) Feinstein, A. D., Montet, B. T., Ansdell, M., et al. 2020, AJ, 160, 219. doi:10.3847/1538-3881/abac0a
  * Feinstein et al. (2022) Feinstein, A. D., Seligman, D. Z., Günther, M. N., et al. 2022, ApJ, 925, L9. doi:10.3847/2041-8213/ac4b5e
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24. doi:10.21105/joss.00024
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., et al. 2017, AJ, 154, 220. doi:10.3847/1538-3881/aa9332
  * Gelman et al. (1995) Gelman, A., Carlin, J. B., Stern, H. S., Rubin, D. B. 1995, Bayesian Data Analysis, Chapman & Hall, New York
  * Gershberg (2005) Gershberg, R. E. 2005, Solar-Type Activity in Main-Sequence Stars: , Astronomy and Astrophysics Library. ISBN 978-3-540-21244-7. Springer Berlin Heidelberg, 2005. doi:10.1007/3-540-28243-2
  * Gilbert et al. (2022) Gilbert, E. A., Barclay, T., Quintana, E. V., et al. 2022, AJ, 163, 147. doi:10.3847/1538-3881/ac23ca
  * Günther et al. (2020) Günther, M. N., Zhan, Z., Seager, S., et al. 2020, AJ, 159, 60. doi:10.3847/1538-3881/ab5d3a
  * Günther & Daylan (2021) Günther, M. N. & Daylan, T. 2021, ApJS, 254, 13. doi:10.3847/1538-4365/abe70e
  * Grindlay & Heise (1975) Grindlay, J. & Heise, J. 1975, International Cosmic Ray Conference, 1, 154
  * Hawley & Pettersen (1991) Hawley, S. L. & Pettersen, B. R. 1991, ApJ, 378, 725. doi:10.1086/170474
  * Hawley et al. (2014) Hawley, S. L., Davenport, J. R. A., Kowalski, A. F., et al. 2014, ApJ, 797, 121. doi:10.1088/0004-637X/797/2/121
  * Hedges et al. (2020) Hedges, C., Angus, R., Barentsen, G., et al. 2020, Research Notes of the American Astronomical Society, 4, 220. doi:10.3847/2515-5172/abd106
  * Hebb et al. (2007) Hebb, L., Petro, L., Ford, H. C., et al. 2007, MNRAS, 379, 63. doi:10.1111/j.1365-2966.2007.11904.x
  * Henry et al. (1995) Henry, G. W., Eaton, J. A., Hamer, J., et al. 1995, ApJS, 97, 513. doi:10.1086/192149
  * Herbst et al. (2021) Herbst, K., Papaioannou, A., Airapetian, V. S., et al. 2021, ApJ, 907, 89. doi:10.3847/1538-4357/abcc04
  * Honda et al. (2018) Honda, S., Notsu, Y., Namekata, K., et al. 2018, PASJ, 70, 62. doi:10.1093/pasj/psy055
  * Hukushima (1999) Hukushima, K. 1999, Phys. Rev. E, 60, 3606. doi:10.1103/PhysRevE.60.3606
  * Ikuta et al. (2020) Ikuta, K., Maehara, H., Notsu, Y., et al. 2020, ApJ, 902, 73. doi:10.3847/1538-4357/abae5f (Paper I)
  * Jeffers et al. (2022) Jeffers, S. V., Barnes, J. R., Schöfer, P., et al. 2022, A&A, 663, A27. doi:10.1051/0004-6361/202141880
  * Johnson et al. (2021) Johnson, L. J., Norris, C. M., Unruh, Y. C., et al. 2021, MNRAS, 504, 4751. doi:10.1093/mnras/stab1190
  * Kalas et al. (2004) Kalas, P., Liu, M. C., & Matthews, B. C. 2004, Science, 303, 1990. doi:10.1126/science.1093420
  * Kass & Raftery (1995) Kass R. E., & Raftery A. E., 1995, J. Am. Stat. Assoc., 90, 773
  * Karoff et al. (2016) Karoff, C., Knudsen, M. F., De Cat, P., et al. 2016, Nature Communications, 7, 11058. doi:10.1038/ncomms11058
  * Kipping (2012) Kipping, D. M. 2012, MNRAS, 427, 2487. doi:10.1111/j.1365-2966.2012.22124.x
  * Klein et al. (2021) Klein, B., Donati, J.-F., Moutou, C., et al. 2021, MNRAS, 502, 188. doi:10.1093/mnras/staa3702
  * Klein et al. (2022) Klein, B., Zicher, N., Kavanagh, R. D., et al. 2022, MNRAS, 512, 5067. doi:10.1093/mnras/stac761
  * Koch et al. (2010) Koch, D. G., Borucki, W. J., Basri, G., et al. 2010, ApJ, 713, L79. doi:10.1088/2041-8205/713/2/L79
  * Kowalski et al. (2010) Kowalski, A. F., Hawley, S. L., Holtzman, J. A., et al. 2010, ApJ, 714, L98. doi:10.1088/2041-8205/714/1/L98
  * Kochukhov (2021) Kochukhov, O. 2021, A&A Rev., 29, 1. doi:10.1007/s00159-020-00130-3
  * Kundu et al. (1987) Kundu, M. R., Jackson, P. D., White, S. M., et al. 1987, ApJ, 312, 822. doi:10.1086/164928
  * Kunkel (1970) Kunkel, W. E. 1970, PASP, 82, 1341. doi:10.1086/129036
  * Lacy et al. (1976) Lacy, C. H., Moffett, T. J., & Evans, D. S. 1976, ApJS, 30, 85. doi:10.1086/190358
  * Luger et al. (2021a) Luger, R., Foreman-Mackey, D., Hedges, C., et al. 2021, AJ, 162, 123. doi:10.3847/1538-3881/abfdb8
  * Luger et al. (2021b) Luger, R., Foreman-Mackey, D., & Hedges, C. 2021, AJ, 162, 124. doi:10.3847/1538-3881/abfdb9
  * Luger et al. (2021c) Luger, R., Bedell, M., Foreman-Mackey, D., et al. 2021, arXiv:2110.06271
  * Maehara et al. (2012) Maehara, H., Shibayama, T., Notsu, S., et al. 2012, Nature, 485, 478
  * Maehara et al. (2017) Maehara, H., Notsu, Y., Notsu, S., et al. 2017, PASJ, 69, 41. doi:10.1093/pasj/psx013
  * Maehara et al. (2021) Maehara, H., Notsu, Y., Namekata, K., et al. 2021, PASJ, 73, 44. doi:10.1093/pasj/psaa098
  * Martioli et al. (2020) Martioli, E., Hébrard, G., Moutou, C., et al. 2020, A&A, 641, L1. doi:10.1051/0004-6361/202038695
  * Messina & Guinan (2002) Messina, S. & Guinan, E. F. 2002, A&A, 393, 225. doi:10.1051/0004-6361:20021000
  * Mitra-Kraev et al. (2005) Mitra-Kraev, U., Harra, L. K., Güdel, M., et al. 2005, A&A, 431, 679. doi:10.1051/0004-6361:20041201
  * Morin et al. (2008) Morin, J., Donati, J.-F., Petit, P., et al. 2008, MNRAS, 390, 567. doi:10.1111/j.1365-2966.2008.13809.x
  * Morris et al. (2017) Morris, B. M., Hebb, L., Davenport, J. R. A., et al. 2017, ApJ, 846, 99. doi:10.3847/1538-4357/aa8555
  * Morris et al. (2018) Morris, B. M., Agol, E., Davenport, J. R. A., et al. 2018, ApJ, 857, 39. doi:10.3847/1538-4357/aab6a5
  * Muheki et al. (2020) Muheki, P., Guenther, E. W., Mutabazi, T., et al. 2020, MNRAS, 499, 5047. doi:10.1093/mnras/staa3152
  * Namekata et al. (2020a) Namekata, K., Davenport, J. R. A., Morris, B. M., et al. 2020, ApJ, 891, 103. doi:10.3847/1538-4357/ab7384
  * Namekata et al. (2020b) Namekata, K., Maehara, H., Sasaki, R., et al. 2020, PASJ, 72, 68. doi:10.1093/pasj/psaa051
  * Namekata et al. (2022a) Namekata, K., Maehara, H., Honda, S., et al. 2021, Nature Astronomy, 6, 241. doi:10.1038/s41550-021-01532-8
  * Namekata et al. (2022b) Namekata, K., Maehara, H., Honda, S., et al. 2022, ApJ, 926, L5. doi:10.3847/2041-8213/ac4df0
  * Notsu et al. (2013) Notsu, Y., Shibayama, T., Maehara, H., et al. 2013, ApJ, 771, 127. doi:10.1088/0004-637X/771/2/127
  * Notsu et al. (2019) Notsu, Y., Maehara, H., Honda, S., et al. 2019, ApJ, 876, 58. doi:10.3847/1538-4357/ab14e6
  * Okamoto et al. (2021) Okamoto, S., Notsu, Y., Maehara, H., et al. 2021, ApJ, 906, 72. doi:10.3847/1538-4357/abc8f5
  * Osten et al. (2005) Osten, R. A., Hawley, S. L., Allred, J. C., et al. 2005, ApJ, 621, 398. doi:10.1086/427275
  * Paudel et al. (2021) Paudel, R. R., Barclay, T., Schlieder, J. E., et al. 2021, ApJ, 922, 31. doi:10.3847/1538-4357/ac1946
  * Pettersen (1980) Pettersen, B. R. 1980, AJ, 85, 871. doi:10.1086/112750
  * Pettersen et al. (1983) Pettersen, B. R., Kern, G. A., & Evans, D. S. 1983, A&A, 123, 184
  * Plavchan et al. (2020) Plavchan, P., Barclay, T., Gagné, J., et al. 2020, Nature, 582, 497. doi:10.1038/s41586-020-2400-z
  * Rackham et al. (2018) Rackham, B. V., Apai, D., & Giampapa, M. S. 2018, ApJ, 853, 122. doi:10.3847/1538-4357/aaa08c
  * Rackham et al. (2022) Rackham, B. V., Espinoza, N., Berdyugina, S. V., et al. 2022, arXiv:2201.09905
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003. doi:10.1117/1.JATIS.1.1.014003
  * Rodono et al. (1986) Rodono, M., Cutispoto, G., Pazzani, V., et al. 1986, A&A, 165, 135
  * Roettenbacher et al. (2017) Roettenbacher, R. M., Monnier, J. D., Korhonen, H., et al. 2017, ApJ, 849, 120. doi:10.3847/1538-4357/aa8ef7
  * Rojas-Ayala et al. (2012) Rojas-Ayala, B., Covey, K. R., Muirhead, P. S., et al. 2012, ApJ, 748, 93. doi:10.1088/0004-637X/748/2/93
  * Shibata et al. (2013) Shibata, K., Isobe, H., Hillier, A., et al. 2013, PASJ, 65, 49. doi:10.1093/pasj/65.3.49
  * Shibayama et al. (2013) Shibayama, T., Maehara, H., Notsu, S., et al. 2013, ApJS, 209, 5. doi:10.1088/0067-0049/209/1/5
  * Silverberg et al. (2016) Silverberg, S. M., Kowalski, A. F., Davenport, J. R. A., et al. 2016, ApJ, 829, 129. doi:10.3847/0004-637X/829/2/129
  * Spangler et al. (1974) Spangler, S. R., Shawhan, S. D., & Rankin, J. M. 1974, ApJ, 190, L129. doi:10.1086/181523
  * Stassun et al. (2018) Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018, AJ, 156, 102. doi:10.3847/1538-3881/aad050
  * Strassmeier & Bopp (1992) Strassmeier, K. G. & Bopp, B. W. 1992, A&A, 259, 183
  * Strassmeier et al. (1999) Strassmeier, K. G., Lupinek, S., Dempsey, R. C., et al. 1999, A&A, 347, 212
  * Strassmeier (2009) Strassmeier, K. G. 2009, A&A Rev., 17, 251. doi:10.1007/s00159-009-0020-6
  * Toriumi et al. (2020) Toriumi, S., Airapetian, V. S., Hudson, H. S., et al. 2020, ApJ, 902, 36. doi:10.3847/1538-4357/abadf9
  * Torres et al. (1972) Torres, C. A. O., Ferraz Mello, S., & Quast, G. R. 1972, Astrophys. Lett., 11, 13
  * Tu et al. (2020) Tu, Z.-L., Yang, M., Zhang, Z. J., et al. 2020, ApJ, 890, 46. doi:10.3847/1538-4357/ab6606
  * Tu et al. (2021) Tu, Z.-L., Yang, M., Wang, H.-F., et al. 2021, ApJS, 253, 35. doi:10.3847/1538-4365/abda3c
  * Vousden et al. (2016) Vousden, W. D., Farr, W. M., & Mandel, I. 2016, MNRAS, 455, 1919. doi:10.1093/mnras/stv2422
  * Walkowicz et al. (2013) Walkowicz, L. M., Basri, G., & Valenti, J. A. 2013, ApJS, 205, 17. doi:10.1088/0067-0049/205/2/17
  * Wisniewski et al. (2019) Wisniewski, J. P., Kowalski, A. F., Davenport, J. R. A., et al. 2019, ApJ, 883, L8. doi:10.3847/2041-8213/ab40bf
  * Zboril (2003) Zboril, M. 2003, Astronomische Nachrichten, 324, 527. doi:10.1002/asna.200310162

##  Appendix A Another models with different assumptions

To evaluate the uncertainties of parameters (Section 2.2), we also perform
starspot modeling with the two-spot model for the light curve of EV Lac Cycle
2 with different assumptions: the quadratic limb-darkening law in Section A.1,
parameters including the spot relative intensity under uniform and normal
prior distributions in Section A.2, and the two-spot model with systematic
noise model by Gaussian process in Section A.3. The deduced parameters are
listed in Table 7.

###  A.1 Quadratic limb-darkening law

To compare the result by the quadratic limb-darkening law with that by the
nonlinear law (Table 2.2), we also fix the coefficient values $c_{1}=c_{3}=0$,
$c_{2}=u_{1}+2u_{2}=1.04$, and $c_{4}=-u_{2}=-0.44$ in Equation 4, where
$u_{1}=0.16$ and $u_{2}=0.44$ (Equation 17 in Paper I) from Claret (2018). As
a result, the uncertainties of spot parameters are the same degree with those
in the two-spot model (Table 6). The uncertainties of the stellar parameters
increase by an order of magnitude, but remain small. The deduced latitude of
the second spot is lower because there is degeneracies between the spot
latitude and size due to the high inclination angle $=60$ (deg) (Paper I).

###  A.2 Spot relative intensity with uniform and normal prior distributions

To take account of the degeneracies between spot latitude and spot relative
intensity, we also include it as parameters under uniform and normal prior
distributions extended to the range from 0 to 2 for bright spots as reported
in Morris et al. (2018). As a result, the uncertainties of stellar parameters
are also the same degree with those in the two-spot model (Table 6). The
uncertainties of spot parameters increase by a factor, but remain small. The
deduced latitude of the second spot is lower for the same reason in Section
A.1, and the deduced radii are larger because the deduced relative intensity
is larger (Section 2.2). The deduced spot intensities $f_{\rm spot}=0.833$ and
$0.811$ for the uniform and normal prior distributions correspond to the spot
temperature $T_{\rm spot}=3284$ and $3268$ (K), respectively (Equation 3).

###  A.3 Two-spot and systematic noise models

To evaluate systematic noise in the residual from the two-spot model, and we
simultaneously perform modeling of the noise with Gaussian process (GP;
Foreman-Mackey et al., 2017; Günther & Daylan, 2021). Then, we employ a
multivariate normal likelihood function with the Matérn-3/2 kernel
characterized by hyperparameters $\sigma_{\rm sys}$ and $\rho_{\rm sys}$:

$p({\cal
D}|\theta)=\frac{1}{(\sqrt{2\pi})^{N}\sqrt{|{\boldsymbol{K}}+\sigma_{i}^{2}{\boldsymbol{I}}|}}\exp{\biggl{[}-\frac{1}{2}\bigl{(}{\boldsymbol{F}}_{{\rm
obs}}-{\boldsymbol{F}}_{{\rm mod}}({\theta})\bigr{)}^{\rm
T}({\boldsymbol{K}}+\sigma_{i}^{2}{\boldsymbol{I}})^{-1}\bigl{(}{\boldsymbol{F}}_{{\rm
obs}}-{\boldsymbol{F}}_{{\rm mod}}({\theta})\bigr{)}\biggr{]}}$ (A1)

and

$K(t_{i},t_{j})=\sigma_{\rm
sys}^{2}\bigl{(}1+\frac{\sqrt{3}|t_{i}-t_{j}|}{\rho_{\rm
sys}}\bigr{)}\exp\bigl{(}-\frac{\sqrt{3}|t_{i}-t_{j}|}{\rho_{\rm
sys}}\bigr{)},$ (A2)

where $\theta$, ${\cal D}$, $\sigma_{i}$, ${\boldsymbol{F}}_{{\rm obs}}$, and
${\boldsymbol{F}}_{{\rm mod}}$ denote parameters, observed data, photometric
error, the vector of the observed flux, and that of the model flux,
respectively. As a result, the uncertainties of spot parameters are slightly
larger than those in the two-spot model (Table 6) up to an order of magnitude
because the systematic noise model enables to optimize the residual from the
spot model flexibly. The uncertainties of stellar parameters increase by two
order of magnitude, but remain small. The deduced latitude of the second spot
is also lower for the same reason in Section A.1.

Figure 1: (a, b, c, d, e) TESS light curves of AU Mic Cycle 1 (Sector 1, July 25 to August 22 in 2018) and 3 (Sector 27, July 4 to 30 in 2020), YZ CMi Cycle 1 (Sector 7, January 7 to February 2 in 2019) and 3 (Sector 34, January 13 to February 9), and EV Lac Cycle 2 (Sector 16, September 11 to October 7 in 2019) (black), respectively. The relative flux is calculated by the PDC-SAP flux normalized by its average excluding missing values. Detected flares are represented in each peak time (red bar) (for details, Section 2.1), and the large flares are marked as outliers (red triangle). The reproduced light curves with the two-spot models are also exhibited under the TESS light curves (red).  Figure 2: Flare frequency for AU Mic Cycle 1 (red) and 3 (blue). The number of flares per observation period of a TESS Sector (Sector 1 and 27) are 55/25.131 = 2.189 and 70/23.307 = 3.003 (day-1), respectively. (a) Cumulative flare frequency distribution versus bolometric flare energy: the number of flares per day with an energy larger than the energy value. (b) Flare frequency distribution versus bolometric flare energy: the number of flares per day and energy in each energy bin. The error bar equals an square root of (the number of flares in each energy bin + 1) per observation period and energy in each energy bin.  Figure 3: Same as Figure 2 but for YZ CMi Cycle 1 and 3. The number of flares per observation period of a TESS Sector (Sector 7 and 34) are 124/22.725 = 5.457 and 137/23.398 = 5.855 (day-1), respectively. Figure 4: Same as Figure 2 but for EV Lac Cycle 2. The number of flares per observation period of a TESS Sector (Sector 16) is 76/19.650 = 3.868 (day-1). Table 2: AU Mic Cycle 1 case  Deduced Parameters | Two-spot Model | Three-spot Model | Prior Distributionaa ${\cal U_{\text{log}}}(a,b)=1/(\theta\log(b/a))$ and ${\cal U}(a,b)=1/(b-a)$ represent bounded log uniform distribution (Jeffreys prior) and bounded uniform distribution defined in $a\leq\theta\leq b$, respectively.
---|---|---|---
(Stellar parameters) |  |  |
1\. Equatorial period $P_{{\scriptsize\textrm{eq}}}$ (day) | $4.8624^{+0.0001}_{-0.0001}$ | $4.8089^{+0.0010}_{-0.0006}$ | ${\cal U}_{\log}(4.5000,5.0000)$
2\. Degree of differential rotation $\kappa$ | $-0.0001^{+0.0001}_{-0.0001}$ | $0.1998^{+0.0002}_{-0.0022}$ | ${\cal U}(-0.2000,0.2000)$
(Spot parameters) |  |  |
(1st spot) |  |  |
3\. Latitude $\Phi_{1}$ (deg) | $-13.54^{+0.12}_{-0.11}$ | $-14.72^{+0.11}_{-0.07}$ |  ${\cal U}(-90.00,\Phi_{2})$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
4\. Initial longitude $\Lambda_{1}$ (deg) | $-22.95^{+0.02}_{-0.03}$ | $-134.06^{+0.05}_{-0.06}$ | ${\cal U}(-180.00,180.00)$
5\. Maximum radius $\alpha_{{\scriptsize\textrm{max,1}}}$ (deg) | $15.61^{+0.01}_{-0.02}$ | $12.86^{+0.02}_{-0.01}$ | ${\cal U}(0.01,30.00)$
(2nd spot) |  |  |
6\. Latitude $\Phi_{2}$ (deg) | $39.67^{+0.08}_{-0.07}$ | $12.08^{+0.13}_{-0.09}$ |  ${\cal U}(\Phi_{1},\Phi_{3})$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
7\. Initial longitude $\Lambda_{2}$ (deg) | $-152.49^{+0.02}_{-0.02}$ | $179.39^{+0.06}_{-0.12}$ | ${\cal U}(-180.00,180.00)$
8\. Maximum radius $\alpha_{{\scriptsize\textrm{max,2}}}$ (deg) | $16.02^{+0.01}_{-0.01}$ | $10.90^{+0.01}_{-0.02}$ | ${\cal U}(0.01,30.00)$
(3rd spot) |  |  |
9\. Latitude $\Phi_{3}$ (deg) | - | $13.92^{+0.10}_{-0.10}$ |  ${\cal U}(\Phi_{2},90.00)$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
10\. Initial longitude $\Lambda_{3}$ (deg) | - | $-24.31^{+0.03}_{-0.02}$ | ${\cal U}(-180.00,180.00)$
11\. Maximum radius $\alpha_{{\scriptsize\textrm{max,3}}}$ (deg) | - | $14.47^{+0.01}_{-0.01}$ | ${\cal U}(0.01,30.00)$
Average flux $F_{\rm ave}$ | 0.831 | 0.830 |
Logarithm of model evidence $\log{\cal Z}$ | $-63950.957$ | $\boldsymbol{-41382.031}$ |
Table 3: AU Mic Cycle 3 case  Deduced Parameters | Two-spot Model | Three-spot Model | Prior Distributionaa ${\cal U_{\text{log}}}(a,b)=1/(\theta\log(b/a))$ and ${\cal U}(a,b)=1/(b-a)$ represent bounded log uniform distribution (Jeffreys prior) and bounded uniform distribution defined in $a\leq\theta\leq b$, respectively.
---|---|---|---
(Stellar parameters) |  |  |
1\. Equatorial period $P_{{\scriptsize\textrm{eq}}}$ (day) | $4.8289^{+0.0001}_{-0.0002}$ | $5.2004^{+0.0112}_{-0.0061}$ | ${\cal U}_{\log}(4.5000,5.0000)$
2\. Degree of differential rotation $\kappa$ | $0.0094^{+0.0001}_{-0.0001}$ | $-0.1072^{+0.0018}_{-0.0032}$ | ${\cal U}(-0.2000,0.2000)$
(Spot parameters) |  |  |
(1st spot) |  |  |
3\. Latitude $\Phi_{1}$ (deg) | $-0.55^{+0.18}_{-0.15}$ | $-59.04^{+0.06}_{-0.06}$ |  ${\cal U}(-90.00,\Phi_{2})$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
4\. Initial longitude $\Lambda_{1}$ (deg) | $-14.36^{+0.03}_{-0.02}$ | $-1.32^{+0.06}_{-0.07}$ | ${\cal U}(-180.00,180.00)$
5\. Maximum radius $\alpha_{{\scriptsize\textrm{max,1}}}$ (deg) | $12.80^{+0.01}_{-0.01}$ | $30.00^{+0.00}_{-0.01}$ | ${\cal U}(0.01,30.00)$
(2nd spot) |  |  |
6\. Latitude $\Phi_{2}$ (deg) | $80.22^{+0.01}_{-0.01}$ | $52.57^{+0.12}_{-0.10}$ |  ${\cal U}(\Phi_{1},\Phi_{3})$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
7\. Initial longitude $\Lambda_{2}$ (deg) | $-153.29^{+0.04}_{-0.03}$ | $164.28^{+0.13}_{-0.14}$ | ${\cal U}(-180.00,180.00)$
8\. Maximum radius $\alpha_{{\scriptsize\textrm{max,2}}}$ (deg) | $30.00^{+0.00}_{-0.01}$ | $16.42^{+0.02}_{-0.04}$ | ${\cal U}(0.01,30.00)$
(3rd spot) |  |  |
9\. Latitude $\Phi_{3}$ (deg) | - | $56.28^{+0.05}_{-0.05}$ |  ${\cal U}(\Phi_{2},90.00)$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
10\. Initial longitude $\Lambda_{3}$ (deg) | - | $-96.59^{+0.09}_{-0.11}$ | ${\cal U}(-180.00,180.00)$
11\. Maximum radius $\alpha_{{\scriptsize\textrm{max,3}}}$ (deg) | - | $16.13^{+0.03}_{-0.03}$ | ${\cal U}(0.01,30.00)$
Average flux $F_{\rm ave}$ | 0.816 | 0.823 |
Logarithm of model evidence $\log{\cal Z}$ | $-139256.404$ | $\boldsymbol{-85746.962}$ |
Table 4: YZ CMi Cycle 1 case  Deduced Parameters | Two-spot Model | One-spot Model | Prior Distributionaa ${\cal U_{\text{log}}}(a,b)=1/(\theta\log(b/a))$ and ${\cal U}(a,b)=1/(b-a)$ represent bounded log uniform distribution (Jeffreys prior) and bounded uniform distribution defined in $a\leq\theta\leq b$, respectively.
---|---|---|---
(Stellar parameters) |  |  |
1\. Equatorial period $P_{{\scriptsize\textrm{eq}}}$ (day) | $2.7753^{+0.0002}_{-0.0003}$ | $2.7737^{+0.0001}_{-0.0001}$ | ${\cal U}_{\log}(2.0000,3.0000)$
2\. Degree of differential rotation $\kappa$ | $-0.0008^{+0.0002}_{-0.0002}$ |  $0.0000$ (fixed) | ${\cal U}(-0.2000,0.2000)$
(Spot parameters) |  |  |
(1st spot) |  |  |
3\. Latitude $\Phi_{1}$ (deg) | $-6.60^{+0.24}_{-0.19}$ | $79.03^{+0.01}_{-0.02}$ |  ${\cal U}(-90.00,\Phi_{2})$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the one-spot model, we set $\Phi_{2}=90.0$ (the upper limit of the latitude).
4\. Initial longitude $\Lambda_{1}$ (deg) | $113.68^{+0.20}_{-0.22}$ | $31.35^{+0.06}_{-0.06}$ | ${\cal U}(-180.00,180.00)$
5\. Maximum radius $\alpha_{{\scriptsize\textrm{max,1}}}$ (deg) | $12.73^{+0.07}_{-0.08}$ | $29.99^{+0.01}_{-0.02}$ | ${\cal U}(0.01,30.00)$
(2nd spot) |  |  |
6\. Latitude $\Phi_{2}$ (deg) | $41.17^{+0.30}_{-0.28}$ | - |  ${\cal U}(\Phi_{1},90.00)$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the one-spot model, we set $\Phi_{2}=90.0$ (the upper limit of the latitude).
7\. Initial longitude $\Lambda_{2}$ (deg) | $17.36^{+0.12}_{-0.18}$ | - | ${\cal U}(-180.00,180.00)$
8\. Maximum radius $\alpha_{{\scriptsize\textrm{max,2}}}$ (deg) | $15.92^{+0.01}_{-0.02}$ | - | ${\cal U}(0.01,30.00)$
Average flux $F_{\rm ave}$ | 0.839 | 0.778 |
Logarithm of model evidence $\log{\cal Z}$ | $\boldsymbol{14982.598}$ | $4625.639$ |
Table 5: YZ CMi Cycle 3 case  Deduced Parameters | Two-spot Model | Three-spot Model | Prior Distributionaa ${\cal U_{\text{log}}}(a,b)=1/(\theta\log(b/a))$ and ${\cal U}(a,b)=1/(b-a)$ represent bounded log uniform distribution (Jeffreys prior) and bounded uniform distribution defined in $a\leq\theta\leq b$, respectively.
---|---|---|---
(Stellar parameters) |  |  |
1\. Equatorial period $P_{{\scriptsize\textrm{eq}}}$ (day) | $2.7705^{+0.0003}_{-0.0002}$ | $2.7707^{+0.0003}_{-0.0003}$ | ${\cal U}_{\log}(2.0000,3.0000)$
2\. Degree of differential rotation $\kappa$ | $0.0082^{+0.0007}_{-0.0008}$ | $0.0088^{+0.0016}_{-0.0011}$ | ${\cal U}(-0.2000,0.2000)$
(Spot parameters) |  |  |
(1st spot) |  |  |
3\. Latitude $\Phi_{1}$ (deg) | $18.48^{+0.28}_{-0.30}$ | $-10.92^{+0.50}_{-0.45}$ |  ${\cal U}(-90.00,\Phi_{2})$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
4\. Initial longitude $\Lambda_{1}$ (deg) | $-23.68^{+0.17}_{-0.15}$ | $133.41^{+0.81}_{-0.66}$ | ${\cal U}(-180.00,180.00)$
5\. Maximum radius $\alpha_{{\scriptsize\textrm{max,1}}}$ (deg) | $12.91^{+0.06}_{-0.05}$ | $7.91^{+0.14}_{-0.19}$ | ${\cal U}(0.01,30.00)$
(2nd spot) |  |  |
6\. Latitude $\Phi_{2}$ (deg) | $38.82^{+1.21}_{-0.81}$ | $16.36^{+0.36}_{-0.29}$ |  ${\cal U}(\Phi_{1},\Phi_{3})$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
7\. Initial longitude $\Lambda_{2}$ (deg) | $170.27^{+0.50}_{-0.50}$ | $-17.50^{+0.35}_{-0.52}$ | ${\cal U}(-180.00,180.00)$
8\. Maximum radius $\alpha_{{\scriptsize\textrm{max,2}}}$ (deg) | $7.77^{+0.12}_{-0.10}$ | $12.34^{+0.05}_{-0.05}$ | ${\cal U}(0.01,30.00)$
(3rd spot) |  |  |
9\. Latitude $\Phi_{3}$ (deg) | - | $49.48^{+5.76}_{-2.69}$ |  ${\cal U}(\Phi_{2},90.00)$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
10\. Initial longitude $\Lambda_{3}$ (deg) | - | $-146.64^{+2.46}_{-2.57}$ | ${\cal U}(-180.00,180.00)$
11\. Maximum radius $\alpha_{{\scriptsize\textrm{max,3}}}$ (deg) | - | $6.59^{+0.41}_{-0.18}$ | ${\cal U}(0.01,30.00)$
Average flux $F_{\rm ave}$ | 0.847 | 0.848 |
Logarithm of model evidence $\log{\cal Z}$ | $13111.858$ | $\boldsymbol{13492.544}$ |
Table 6: EV Lac Cycle 2 case  Deduced Parameters | Two-spot Model | Three-spot Model | Prior Distributionaa ${\cal U_{\text{log}}}(a,b)=1/(\theta\log(b/a))$ and ${\cal U}(a,b)=1/(b-a)$ represent bounded log uniform distribution (Jeffreys prior) and bounded uniform distribution defined in $a\leq\theta\leq b$, respectively.
---|---|---|---
(Stellar parameters) |  |  |
1\. Equatorial period $P_{{\scriptsize\textrm{eq}}}$ (day) | $4.3573^{+0.0002}_{-0.0003}$ | $4.3651^{+0.0003}_{-0.0002}$ | ${\cal U}_{\log}(4.0000,5.0000)$
2\. Degree of differential rotation $\kappa$ | $0.0020^{+0.0001}_{-0.0001}$ | $-0.0033^{+0.0001}_{-0.0001}$ | ${\cal U}(-0.2000,0.2000)$
(Spot parameters) |  |  |
(1st spot) |  |  |
3\. Latitude $\Phi_{1}$ (deg) | $24.89^{+0.28}_{-0.21}$ | $9.91^{+0.12}_{-0.18}$ |  ${\cal U}(-90.00,\Phi_{2})$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
4\. Initial longitude $\Lambda_{1}$ (deg) | $60.35^{+0.11}_{-0.13}$ | $-113.64^{+0.20}_{-0.20}$ | ${\cal U}(-180.00,180.00)$
5\. Maximum radius $\alpha_{{\scriptsize\textrm{max,1}}}$ (deg) | $13.44^{+0.02}_{-0.02}$ | $14.67^{+0.04}_{-0.06}$ | ${\cal U}(0.01,30.00)$
(2nd spot) |  |  |
6\. Latitude $\Phi_{2}$ (deg) | $55.89^{+0.31}_{-0.26}$ | $9.95^{+0.16}_{-0.13}$ |  ${\cal U}(\Phi_{1},\Phi_{3})$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
7\. Initial longitude $\Lambda_{2}$ (deg) | $-84.11^{+0.11}_{-0.11}$ | $111.85^{+0.32}_{-0.35}$ | ${\cal U}(-180.00,180.00)$
8\. Maximum radius $\alpha_{{\scriptsize\textrm{max,2}}}$ (deg) | $14.36^{+0.06}_{-0.06}$ | $10.56^{+0.07}_{-0.09}$ | ${\cal U}(0.01,30.00)$
(3rd spot) |  |  |
9\. Latitude $\Phi_{3}$ (deg) | - | $54.15^{+0.53}_{-0.47}$ |  ${\cal U}(\Phi_{2},90.00)$bb We discern each spot by its latitude $\Phi_{k}$, not by its reference time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT. In case of the two-spot model, we set $\Phi_{3}=90.0$ (the upper limit of the latitude).
10\. Initial longitude $\Lambda_{3}$ (deg) | - | $13.71^{+0.36}_{-0.25}$ | ${\cal U}(-180.00,180.00)$
11\. Maximum radius $\alpha_{{\scriptsize\textrm{max,3}}}$ (deg) | - | $15.99^{+0.06}_{-0.04}$ | ${\cal U}(0.01,30.00)$
Average flux $F_{\rm ave}$ | 0.850 | 0.842 |
Logarithm of model evidence $\log{\cal Z}$ | $-7297.802$ | $\boldsymbol{-4135.585}$ |
Figure 5: (Left) Visualized surface reproduced with each mode of posterior
distributions for the rotation phases = 0.1, 0.3, 0.5, 0.7, and 0.9 (vertical
dashed line). The light curve is folded with the period = 4.862 (day) from the
day at the first peak = 1329.1 (BJD-2457000): (a) the phase-folded light curve
of AU Mic Cycle 1 (black) and phase-folded reproduced one with the two-spot
model (red); (b) the temporal variation of visible projected area of each spot
(red, blue, and green) and the total area (black) relative to the stellar disk
or solar hemisphere; (c) same as (b) but for the visible area relative to the
stellar or solar hemispheres; (d) flare frequency per day in each bin for all
flares (red) and flares with larger equivalent duration than $1.0$ (sec)
(blue). The error bar equals a square root of (the number of flares in each
bin + 1) per observation period; and (e) flare energy and the corresponding
equivalent duration in each bin. (Right) (f, g, h) Same as the left but for
the three-spot model.  Figure 6: Same as Figure 5 but for AU Mic Cycle 3 with
the period = 4.829 (day) from the day at the first peak = 2040.1
(BJD-2457000).  Figure 7: (Left) Visualized surface reproduced with each mode
of posterior distributions for the rotation phases = 0.1, 0.3, 0.5, 0.7, and
0.9 (vertical dashed line). The light curve is folded with the period = 2.774
(day) from the day at the first peak = 1498.3 (BJD-2457000): (a) the phase-
folded light curve of YZ CMi Cycle 1 (black) with large flares marked as
outliers (red triangle) and phase-folded reproduced one with the two-spot
model (red); (b) the temporal variation of visible projected area of each spot
(red, blue, and green) and the total area (black) relative to the stellar disk
or solar hemisphere; (c) same as (b) but for the visible area relative to the
stellar or solar hemispheres; (d) flare frequency per day in each bin for all
flares (red) and flares with larger equivalent duration than $1.0$ (sec)
(blue). The error bar equals a square root of (the number of flares in each
bin + 1) per observation period; and (e) flare energy and the corresponding
equivalent duration in each bin. (Right) Same as the left but for the one-spot
model.  Figure 8: Same as Figure 7 but for YZ CMi Cycle 3 with the period =
2.773 (day) from the day at the first peak = 2231.1 (BJD-2457000). Figure 9:
(Left) Visualized surface reproduced with each mode of posterior distributions
for the rotation phases = 0.1, 0.3, 0.5, 0.7, and 0.9 (vertical dashed line).
The light curve is folded with the period = 4.359 (day) from the day at the
first peak = 1741.1 (BJD-2457000): (a) the phase-folded light curve of EV Lac
Cycle 2 (black) with large flares marked as outliers (red triangle) and Phase-
folded reproduced one with the two-spot model (red); (b) the temporal
variation of visible projected area of each spot (red, blue, and green) and
the total area (black) relative to the stellar disk or solar hemisphere; (c)
Same as (b) but for the visible area relative to the stellar or solar
hemispheres; (d) flare frequency per day in each bin for all flares (red) and
flares with larger equivalent duration than $1.0$ (sec) (blue). The error bar
equals a square root of (the number of flares in each bin + 1) per observation
period; and (e) flare energy and the corresponding equivalent duration in each
bin. (Right) Same as the left but for the three-spot model.

Table 7: EV Lac Cycle 2 additional case  Deduced Parameters | Quadratic limb-darkening | Uniform | Normal | Gaussian process | Prior Distributionaa ${\cal U_{\text{log}}}(a,b)=1/(\theta\log(b/a))$, ${\cal U}(a,b)=1/(b-a)$, and ${\cal TN}(\mu,\sigma^{2},a,b)$ represent bounded log uniform distribution (Jeffreys prior), bounded uniform distribution, and truncated normal distribution, defined in $a\leq\theta\leq b$, respectively.
---|---|---|---|---|---
(Stellar parameters) |  |  |  |  |
1\. Equatorial period $P_{{\scriptsize\textrm{eq}}}$ (day) | $4.3482^{+0.0007}_{-0.0010}$ | $4.3554^{+0.0002}_{-0.0004}$ | $4.3559^{+0.0003}_{-0.0003}$ | $4.3518^{+0.0019}_{-0.0106}$ | ${\cal U}_{\log}(4.0000,5.0000)$
2\. Degree of differential rotation $\kappa$ | $0.0147^{+0.0011}_{-0.0010}$ | $0.0037^{+0.0003}_{-0.0002}$ | $0.0032^{+0.0002}_{-0.0001}$ | $0.0120^{+0.0129}_{-0.0026}$ | ${\cal U}(-0.2000,0.2000)$
(Spot parameters) |  |  |  |  |
(1st spot) |  |  |  |  |
3\. Latitude $\Phi_{1}$ (deg) | $22.59^{+0.24}_{-0.20}$ | $25.70^{+0.31}_{-0.39}$ | $25.41^{+0.38}_{-0.33}$ | $21.11^{+1.92}_{-1.52}$ |  ${\cal U}(-90.00,\Phi_{2})$bb We discern each spot by its latitude $\Phi_{k}$, not by its refereznce time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT.
4\. Initial longitude $\Lambda_{1}$ (deg) | $50.30^{+0.13}_{-0.10}$ | $56.56^{+0.27}_{-0.20}$ | $57.28^{+0.22}_{-0.20}$ | $50.30^{+0.73}_{-0.71}$ | ${\cal U}(-180.00,180.00)$
5\. Maximum radius $\alpha_{{\scriptsize\textrm{max,1}}}$ (deg) | $12.41^{+0.01}_{-0.01}$ | $24.35^{+0.21}_{-0.17}$ | $22.78^{+0.20}_{-0.19}$ | $12.46^{+0.05}_{-0.04}$ | ${\cal U}(0.01,30.00)$
(2nd spot) |  |  |  |  |
6\. Latitude $\Phi_{2}$ (deg) | $30.29^{+0.33}_{-0.26}$ | $46.14^{+0.65}_{-0.60}$ | $48.06^{+0.67}_{-0.47}$ | $29.16^{+1.64}_{-1.89}$ |  ${\cal U}(\Phi_{1},90.00)$bb We discern each spot by its latitude $\Phi_{k}$, not by its refereznce time $t_{k}$ in Paper I, to improve the sampling efficiency of the PT.
7\. Initial longitude $\Lambda_{2}$ (deg) | $-85.66^{+0.09}_{-0.12}$ | $-86.97^{+0.17}_{-0.14}$ | $-86.44^{+0.17}_{-0.14}$ | $-85.56^{+0.79}_{-0.72}$ | ${\cal U}(-180.00,180.00)$
8\. Maximum radius $\alpha_{{\scriptsize\textrm{max,2}}}$ (deg) | $11.77^{+0.01}_{-0.01}$ | $24.48^{+0.21}_{-0.21}$ | $23.14^{+0.22}_{-0.22}$ | $11.78^{+0.06}_{-0.04}$ | ${\cal U}(0.01,30.00)$
(Additional parameters) |  |  |  |  |
9\. Relative intensity $f_{\rm spot}$ |  $0.48$ (fixed) | $0.833^{+0.002}_{-0.002}$ | $0.811^{+0.003}_{-0.003}$ |  $0.48$ (fixed) |  ${\cal U}(0.00,2.00)$ / ${\cal TN}(0.48,0.01,0.00,2.00)$
10\. GP amplitude $\ln\sigma_{\rm sys}$ | - | - | - | $-6.5994^{+0.0181}_{-0.0130}$ | ${\cal U}(-10.0,0.0)$
11\. GP exponential length $\ln\rho_{\rm sys}$ | - | - | - | $-5.1250^{+0.0235}_{-0.0198}$ | ${\cal U}(-10.0,0.0)$
Average flux $F_{\rm ave}$ | $0.854$ | $0.851$ | $0.851$ | $0.856$ |
Logarithm of model evidence $\log{\cal Z}$ | $-5995.666$ | $-5654.106$ | $-5753.989$ | $21813.240$ |

##  Appendix B Joint posterior distributions

The joint posterior distributions of the equatorial period $P_{\rm eq}$,
degree of differential rotation $\kappa$, latitude $\Phi_{k}$, initial
longitude $\Lambda_{k}$, and maximum radius
$\alpha_{{\scriptsize\textrm{max},k}}$, are delineated in Figure 10, 11, 12,
13, 14, 15, 16, 17, 18, and 19 for AU Mic Cycle 1 and 3, YZ CMi Cycle 1 and 3,
and EV Lac Cycle 2, generated with corner (Foreman-Mackey, 2016). All figures
show each of the distribution is an unimodal distribution.

Figure 10: The joint posterior distribution of parameters for the light curve
of AU Mic Cycle 1 with the two-spot model. Each column represents the
equatorial period $P_{\rm eq}$, degree of differential rotation $\kappa$,
latitude $\Phi_{k}$, initial longitude $\Lambda_{k}$, and maximum radius
$\alpha_{{\scriptsize\textrm{max},k}}$, respectively.  Figure 11: Same as
Figure 10 but for the three-spot model.  Figure 12: Same as Figure 10 but for
AU Mic Cycle 3.  Figure 13: Same as Figure 12 but for the three-spot model.
Figure 14: The joint posterior distribution of parameters for the light curve
of YZ CMi Cycle 1 with the two-spot model. Each column represents the
equatorial period $P_{\rm eq}$, degree of differential rotation $\kappa$,
latitude $\Phi_{k}$, initial longitude $\Lambda_{k}$, and maximum radius
$\alpha_{{\scriptsize\textrm{max},k}}$, respectively.  Figure 15: Same as
Figure 14 but for the one-spot model.  Figure 16: Same as Figure 14 but for YZ
CMi Cycle 3.  Figure 17: Same as Figure 16 but for the three-spot model.
Figure 18: The joint posterior distribution of parameters for the light curve
of EV Lac Cycle 2 with the two-spot model. Each column represents the
equatorial period $P_{\rm eq}$, degree of differential rotation $\kappa$,
latitude $\Phi_{k}$, initial longitude $\Lambda_{k}$, and maximum radius
$\alpha_{{\scriptsize\textrm{max},k}}$, respectively.  Figure 19: Same as
Figure 18 but for the three-spot model.