Characterization of starspots on a young M-dwarf K2-25: multi-band observations of stellar photometric variability and planetary transits

Detailed atmospheric characterization of exoplanets by transmission spectroscopy requires careful consideration of stellar surface inhomogeneities induced by starspots. This effect is particularly problematic for planetary systems around M-dwarfs, and their spot properties are not fully understood. We investigated the stellar activity of the young M-dwarf K2-25 and its effect on transit observations of the sub-Neptune K2-25b. From multi-band monitoring observations of stellar brightness variability using ground-based telescopes and TESS, we found that the temperature difference between the spots and photosphere is<190 K and the spot covering fraction is<61% (2$\sigma$). We also investigated the effect of starspot activity using multi-epoch, multi-band transit observations. We rule out cases with extremely low spot temperatures and large spot covering fractions. The results suggest that spots could distort the transmission spectrum of K2-25b by as much as $\sim$100 ppm amplitude, corresponding to the precision of JWST/NIRSPEC of the target. Our study demonstrates that simultaneous multi-band observations with current instruments can constrain the spot properties of M-dwarfs with good enough precision to support atmospheric studies of young M-dwarf planets via transmission spectroscopy.


INTRODUCTION
The number of discovered exoplanets now exceeds 5000, largely due to space-based transit surveys, such as the Kepler space telescope (Borucki et al. 2010) and the Transiting Exoplanets Survey Satellite (TESS; Ricker et al. 2015).Among the various planetary groups, planets feasible for transmission spectroscopy are one of the most interesting targets to characterize their atmospheres.The James Webb Space Telescope (JWST; Gardner et al. 2006) and future space telescopes such as Ariel (Tinetti et al. 2018(Tinetti et al. , 2021)), with their extremely ★ E-mail: mayukomori.519@gmail.comhigh precision, enable us to detect atmospheric signals of less than 100 ppm, which allows more detailed analysis of the planetary atmospheric compositions.
Recently, the transit light source effect (TLSE) has been identified as a significant factor affecting planetary transmission spectra (Rackham et al. 2018(Rackham et al. , 2023)).In brief, unocculted spots or faculae during planetary transits make the observed transit depth shallower or deeper than the expected transit depth without spots or faculae.This effect is particularly complicated for the planets around active M-dwarfs: the wavelength dependence of this effect is confounded by the absorption of molecules present in both stellar photosphere and spots, resulting in false positive signals of planetary atmospheres.face models assuming a spot temperature of 2800 K or 3000 K on a stellar surface of 3200 K.They found that the transmission spectrum of K2-25 b prefers a flat model to a clear solar-abundance atmosphere model, as a large spot covering fraction (22% or 36%) was needed to fit a solar-abundance atmosphere model, and such a large spot covering fraction contradicts the observed amplitude of brightness variations.Therefore, they concluded that, although TLSE could be at play to some extent, the atmospheric spectrum of K2-25 b is flat, which indicates a cloudy/hazy atmosphere and/or a high mean molecular weight atmosphere.Gaidos et al. (2020) and Stefansson et al. (2020) measured the stellar obliquity and planetary mass via spectroscopic observations.Their radial velocity observations yielded a planetary mass higher than that expected from a mass-radius relation, making the amplitude of the atmospheric signal in the transmission spectrum smaller without necessitating the presence of clouds or hazes.
No detailed space-based transmission spectroscopy has ever been performed for K2-25, and such observations by JWST and similar instruments are needed to investigate whether the atmospheric spectrum of K2-25 b is flat.On the other hand, high-precision spectroscopy will need to pay more attention to TLSE, as it will be sensitive to small differences in spot properties, which can complicate the interpretation of atmospheric models.Thao et al. (2020) made empirical assumptions about spot temperatures, but more realistic stellar surface models can be constrained by observations.Interest in the K2-25 system itself and the need to understand its current spot activity are the reasons we chose this system as our target, in addition to the fact that it was a good test case for our methods.
The remainder of the paper is organized as follows: Section 2 summarizes the observational data used in our analyses; Section 3 describes our analysis methods for constraining the spot characteristics from multi-band monitoring and transit light curves; Section 4 describes K2-25's spot characteristics and the effect of the spots on the transit observations, and compares our results with previous studies.Section 5 summarizes the results and discussions, with the future prospects of studying spot properties of M-dwarfs.

OBSERVATIONS & DATA REDUCTIONS
The observations are divided into two categories: monitoring photometry to investigate stellar brightness variations using TESS, Sinistro, and ZTF (Section 2.1, 2.2, and 2.3); transit photometry to evaluate the effect of spots on transit depth measurements by TESS, MuSCAT2, and MuSCAT3 (Section 2.1 and 2.4).These observations were conducted at optical wavelengths, where the spot effect is expected to be prominent.The data are summarized in Table 2. K2-25 (TIC 434226736, TOI-5095) was observed by TESS in Sector 44 (2021 October 12 to November 6) during TESS Cycle 4 with a 2-min cadence in the TESS bandpass (600-1000 nm).We used the PDCSAP (Pre-search Data Conditioning Simple Aperture Photometery) light curves produced by the Science Processing Operations Center (SPOC) photometry pipeline (Jenkins 2002;Jenkins et al. 2010Jenkins et al. , 2016) ) for both the monitoring and transit analyses.We also checked the SAP (Simple Aperture Photometry-) light curves, as sometimes stellar modulation is removed from the PDCSAP light curves.We confirmed that the brightness variation on K2-25 was almost identical in the SAP and PDCSAP light curves, both in shape and normalized amplitude, with instrumental systematics and a long-term trend corrected in PDCSAP light curve.In addition, while nearby stars blending into the TESS aperture may affect photometric accuracy, we confirmed that their effects were negligible (Appendix A1).

Monitoring and transit photometry -TESS
To analyse the brightness variability, we removed outliers with more than 3 away from the median of the normalized light curve.In the transit analyses, for the six transits observed by TESS, nondetrended light curves were cut out for a range of 0.1 days before and after the predicted transit-centre time.

Monitoring photometry -LCO 1m / Sinistro
Sinistro is a single-band optical camera mounted on five 1m telescopes worldwide, operated by Las Cumbres Observatory (LCO, Brown et al. 2013).Our observations were performed specifically with the telescopes at McDonald Observatory in Texas and at Teide Observatory in Tenerife, Canary islands.Sinistro has a 26 ′ .5 × 26 ′ .5 field of view with a pixel scale of 0.389 ′′ pix −1 .
We observed K2-25 from 2021 October 02 to November 05 with Sinistro in -band (400-550 nm), and from 2021 October 13 to 2022 January 31 in -band (820-920 nm).With the exposure time of 120 seconds in -band and 40 seconds in   -band (denoted -band hereafter for simplicity), 5-15 frames were taken every 4-5 hours during good observing conditions with airmass < 1.6, moon distance > 30 deg, and no weather or instrumental issues.Data reduction was carried out through the LCO BANZAI pipeline (McCully et al. 2018).The photometry was obtained using the package astropy-Photutils (Bradley et al. 2021).We found that circular apertures with a radius of 8 pixels were optimal to reduce the scatter in the light curves.We also subtracted the sky background noise calculated using an annulus with a radius of 50 to 80 pixels.We adopted the differential photometry method to account for the telluric variability.In practice, the telluric effect on stellar brightness variations should depend on the colour of the star, given that atmospheric extinction occurs to different degrees for stars of different colours.To minimise this effect, we only used stars with colours similar to the target as comparison stars.We calculated the  −  colour of each star in the field of view from the catalogue of Pan-STARRS Data Release 2 (Chambers et al. 2016;Flewelling et al. 2020).As the  −  value of K2-25 is 3.48, we selected stars with  −  > 2 for use in the calculation.
We excluded frames with the following three conditions as outliers: (i) the number of detected stars was smaller than 50; (ii) the median stellar magnitude was > 0.7 mag fainter than that of the reference frame; (iii) the obtained target brightness was more than 3 different from the mean model obtained by fitting a third-order polynomial.

Monitoring photometry -ZTF
The Zwicky Transient Facility (ZTF) is a time-domain survey that conducts a full sweep of the northern visible sky every two days (Bellm et al. 2019).The survey is performed by a camera with a 47 square degree field of view mounted on the Samuel Oschin 1.2m telescope at Palomar Observatory in California (Bellm et al. 2019).
We used the archival photometry of K2-25 in   and   -bands (367.6-561.4 and 549.8-739.4 nm), which is similar to SDSS  and -band (Masci et al. 2019), from 2018 March 27 until 2022 March 2. For the ZTF light curves, we removed the data points with the quality warning flags (expressed as 'catflags' in the ZTF catalogue), and data points where zero-magnitude is 3 away from the median value of zero-magnitude.In addition, the magnitude was corrected for the second-order extinction effect by a linear function of colour coefficients ('clrcoeff' in the ZTF catalogue).Finally, the magnitude was converted to flux, from which we produced a normalised light curve.
We found that the shape of brightness variations of K2-25 had changed over a timescale of ∼ 200 days possibly due to the variation of spot distribution (see Appendix A3).Therefore, in order to maintain simultaneous observations with the other dataset, we only used the data after MJD =59400 (2021 July 5) in the analyses.
We observed four transits of K2-25 b with MuSCAT2 on 2021 December 8, 22, 29, and 2022 January 5, and two transits with MuS-CAT3 on 2022 January 8 and 29.The exposure times and number of data for each observing night are summarised in Table 2.A transit number  T was introduced to indicate how many transits have occurred with the transit at  0 = 2458515.64206(BJD) (Stefansson et al. 2020) as the starting point.We note that the -band camera was not available on 2021 December 8 and 22, because of technical difficulties with the instrument.
These observations were not exactly simultaneous with the TESS observations, and there was a gap of 35 days (∼ 10 planetary orbital periods) between the last transit with TESS and the first transit with MuSCAT2.
The reduction for MuSCAT2 data and the aperture photometry of both MuSCAT2/3 datasets were conducted by the custom pipeline described in Fukui et al. (2011), and the reduction for MuSCAT3 data was carried out by the BANZAI pipeline.For each night and each band, the optimal aperture radius was found by selecting the value between 8 and 14 pixels that minimized the scatter in the light curve.Between one and three stars brighter than the target were selected as comparison stars for each night and each band.
The derived raw light curves are shown in Figure A.9. Apart from the transit signal, some rapid flux increases were detected especially in shorter wavelengths, which are likely to be flares.In particular, we visually detected a significant flare-like slope in  and -band light curves for transit  T = 307.We also detected a flux bump during the transit for light curve  T = 314, which might be a flare during the transit or a spot crossing.
The removal of outliers was conducted as follows.First, we removed a flux jump caused by a flare-like signal in  T = 307.As the highest peak of the signal was at BJD= 2459585.36556, the data points 0.01 day before and after this timestamp were removed.Although the flare-like signal was not visually obvious in  and band, we removed the points in all bands to reduce the bias from arbitrary manipulation.Then, we fit each light curve at each band and each date one by one, using mean models as the transit model combined with the baseline model.The transit model was calculated by the package batman (Kreidberg 2015) from the transit parameters described in Section 3.6.The baseline model was defined by a third-order polynomial of time.The models were optimized using MCMC package emcee (Foreman-Mackey et al. 2013).We removed outliers with more than 3 away from the mean model.We repeated this operation three times (Figure A.9).

ANALYSES & RESULTS
The analyses described in this section are divided into two main parts.First is the analysis of stellar brightness variability to constrain the spot characteristics and to calculate the degree of TLSE (Section 3.1 to 3.5).Second is the multi-band transit depth derivation to evaluate the TLSE (Section 3.6).Additional analyses and results are described in the Appendix.

Monitoring light curve analyses -rotation period
We obtained five light curves: the TESS light curve, Sinistro light curves in -and -band, and ZTF light curves in   -and   -band.These light curves allow us to constrain the temperature and distribution of spots on K2-25.
We derived the periodic signals in the light curves using the Generalized Lomb-Scargle Periodogram (GLS; Zechmeister & Kürster 2009) in the Python package PyPeriod (Czesla et al. 2019).We searched for the period of the light curves between 1 and 20 days and obtained the strong GLS power peak at ∼ 1.88 days in all light curves (Figure 1).False alarm probability (FAP) values below 10 −3 indicate robust detection for those periods.The one-day alias of that signal at ∼ 2.14 days also shows a strong peak in the power, especially in ZTF light curves, because the ZTF data is collected every few nights.
The mean of the derived periods weighted by the errors was calculated to be 1.87708 ± 0.00066 day, and we adopted this value as the stellar rotation period  rot .This value is consistent with the previously reported value  rot = 1.88 ± 0.02 day from the K2 light curve (Mann et al. 2016).

Monitoring light curve analyses -modulation amplitudes
We folded the light curves by the period  rot = 1.87708 days, assuming that the brightness variability was stable during the observation window and there was no red noise (time-correlated noise) in the light curves in all bands.It is possible that instrumental systematics behave as red noise in the light curves, but they are unlikely to be correlated with the stellar rotation phase.The data points were binned every 0.02 day in all bands in the following analyses.
We assumed that all the phase-folded light curves have have local maxima and minima at the same phases and differ only in amplitude, and the light curves can be simply modeled as a summation of sine functions.The derivation of the amplitudes had two steps.First, we optimized the TESS light curve by a summation of the sine functions as the reference model.We used the TESS light curve here as it has the most clear modulational shape.Second, we modelled the other light curves by changing the amplitude and offset of the reference model.
In the first step, the assumed reference model  TESS ( ) along with the stellar rotational phase  is described as where   ,   , and   are the parameters to be optimized.The optimization was performed using scipy.optimize, to make the chi-squared value smallest for the TESS light curve.The optimum parameter values are summarized in Table 3.The number of sine-curves   was selected to be 2 because the Bayesian information criterion (BIC) became smallest when the   value was tested from 1 to 4. In step two, the flux   at each band  ={Sinistro-, ZTF-, ZTF-, and Sinistro-} is represented by the  TESS ( ) with the optimum parameters in Table 3, as in the following equation, where   is the amplitude and (.) is the offset at each band.
For each band independently, we derived the value that maximises the log-likelihood with two parameters,   and (.) .Python package emcee was used for the parameter estimation via the Markov chain Monte Carlo method (MCMC).In this process, we found that the white noise in the light curves are larger than the photometric errors, which may be due to instrumental systematics or observing conditions.In order to correctly account for this noise, we amplified the errors so that the reduced chi-squared value approached 1, and we re-fit the model to the data with amplified errors.Figure 2 summarizes the resulting   values with their uncertainties, and Figure 3  the derived best-fit models with 1 uncertainties.It is clearly seen that the amplitude becomes smaller to longer wavelengths.As the Sinistro- and ZTF- have almost identical wavelength ranges, the weighted-mean value is calculated and used in the following analyses.We also tried another method using a Gaussian process (GP) model, which is more flexible for modelling the light curves, to independently derive the modulation amplitudes (see Appendix A2) and confirmed that the results are consistent.

Derivation of spot temperature and covering fraction
The differences in the amplitude at each wavelength reflect differences in the spot contrast, defined as the spot intensity relative to the photosphere.Based on that, we derive the spot temperature and covering fraction from the modulation amplitudes in this section.
We note that spots visible regardless of the stellar rotation phase (referred to as "always-visible spots" hereafter) are not considered in this analysis.It has been noted that even with multi-band light curves, constraining the covering fraction of always-visible spots is challenging (Rosich et al. 2020).We also confirmed that the estimation of the spot temperature and spot covering fraction, which contributes to the brightness variation, remains unaffected whether the presence of always-visible spots is considered or not.On the other hand, while always-visible spots do not contribute to the brightness variation, they do cause TLSE during planetary transits.We account for them in the TLSE estimation in Section 3.5 by assuming various values of spot covering fraction.
Assuming that spots are scarcely visible during the peak brightness of a star and the spots with temperature  spot and covering fraction  spot during the dimmest phase, the amplitude of stellar brightness variation   at each band  = {, , TESS, } can be expressed by where   are the spot contrast (e.g., Equation 8in Notsu et al. 2013).
The spot contrast   is calculated from the following formula (e.g., Morris et al. 2018;Ikuta et al. 2023): where T  is the filter transmission function and F ,spot/phot are the stellar flux from the spots or the photosphere at the temperature of  spot and  phot , respectively.We calculated   through Equation 4using the filter profiles from the SVO Filter Profile Service (Rodrigo et al. 2012;Rodrigo & Solano 2020) and theoretical stellar spectra from the SVO Theoretical Model Services.For stellar spectra models, we used the BT-Settl model with the surface gravity log  = 5 (in cgs), and the metallicity [Fe/H] = 0 dex.The spot temperature  spot ranged from 2700 K to 3200 K in each 100 K step for the calculation and then interpolated.The photosphere temperature  phot was fixed to be 3200 K, but we confirmed that the ∼ 100 K difference of  phot did not cause significant change to the relation of  spot and   .
Figure 4 shows the relation of calculated   and Δ, which is defined by  phot −  spot , for all observed bands.The   values show banddependent variation when Δ is around 500 K, but the variation gets smaller when Δ is close to 0 K or more than 1000 K, in the case of  phot = 3200 K.The combination of Equations 3 and 4 enables us to estimate the two unknown values,  spot and  spot from the estimated amplitude   in each band.As there are four bands in our analysis,  spot and  spot can be obtained simultaneously by formulating four of these equations with each error.The parameters are set to Δ and  spot , and the best parameter values were searched to describe the observed amplitudes for the four bands by calculating chi-squared values.We used MCMC package emcee for the parameter estimation.We put prior for Δ to take positive values and  spot to take values between 0 to 1.
Figure 5 shows the posterior distribution of the parameters using 0 500 1000 1500 corner (Foreman-Mackey 2016).Although they are highly correlated, we derived the 2 upper limit of Δ is 190 K and  spot is 61%.Note that we did not consider the possibility of spots with temperature higher than the stellar photosphere (referred to as "bright spots").We discuss the case for bright spots in Appendix A4.

Investigation of the spot distribution
To evaluate the analysis of estimating spot temperature and covering fraction only from the amplitudes in multi-band light curves (Section 3.3), we further explored possible spot distributions using a spot mapping method (Ikuta et al. 2020(Ikuta et al. , 2023)).Given that there are known degeneracy in various spot distributions as a solution to the starspot mapping from the light curves (Basri & Shah 2020;Ikuta et al. 2020;Luger et al. 2021a), this study adopted the method of trying out several patterns of spot distribution and exploring their trends rather than searching for a single solution.
We used an analytical spotted model macula (Kipping 2012) as in Ikuta et al. (2023).It calculates the flux relative to the spotless photosphere assuming a particular stellar inclination angle, for a specified number of circular spots characterized by four parameters: spot contrast, radius, latitude, and longitude (for more details, see Ikuta et al. 2023).We calculated the light curves in , , TESS, and -band under different spot contrast and stellar limb-darkening coefficients.The spot temperature was set to be identical for all spots.The spot contrast at each band was calculated by Equation 4 under the assumed temperature.The limb-darkening coefficients were calculated by the package PyLDTk (Parviainen & Aigrain 2015;Husser et al. 2013), assuming the stellar effective temperature  eff = 3207 ± 58 K, metallicity [Fe/H] = 0.15 ± 0.03 dex, and stellar surface gravity log  = 4.944 ± 0.031 (in cgs).The stellar inclination and the degree of differential rotation were set to 90 deg (Stefansson et al. 2020) and 0. The rotation period was set to 1 because the light curves were phase-folded.
We reproduced various patterns of stellar surface models, with the number of spots set to two, three, four, and five.In each spot configuration with the number of spots of   , the light curve model was generated with 3  +1 parameters: radius, latitude, and longitude for each spot and Δ.Since a larger number of parameters require a higher computational cost to obtain the optimal values, we provided some additional settings to simplify the calculations.First, we fixed the  phot at 3200 K. Also, the latitude of any spot was restricted to taking only positive values.It is because the stellar inclination of K2-25 is ∼ 90 deg so it is impossible to determine whether a spot is on the northern or southern hemisphere of the star.For spot size, we set uniform prior so that the maximum spot radius was 30 • .We limited Δ to be a positive value so that faculae (or any bright spots) were not considered.
At each number of spots   , we derived the posterior distributions using emcee to fit the reproduced light curves with the observed ones by calculating the chi-squared values.In the calculation, we adopted the amplified uncertainties of the observed light curves calculated in Section 3.2.
In addition, to compare the results by spot mapping with those from the amplitude, we calculated the spot covering fraction  spot from the posterior distribution, by the summation of the projected area of each spot (Equation 11 in Ikuta et al. 2020) at the minimum of the light curve.
Figure 6 shows the posterior distributions of Δ and  spot at each spot number case using corner (Foreman-Mackey 2016).The 2 upper limit of Δ was 142, 86, 78, 70 K and the upper limit of  spot was 18, 31, 35, 41% for two, three, four and five spots cases, respectively.
While the two-spot case showed slightly different values, the upper limit values were similar for more than three spots.The model evidence was calculated via the importance sampling algorithm (e.g., Díaz et al. 2014) as well as Ikuta et al. (2020).The logarithm of the model evidence was calculated to be −711.34,−707.81,−706.60, and −705.78 for two, three, four, and five-spot models, respectively.The model evidence is almost equal for the three, four, and five-spot cases within a factor of 8 (∼ exp(2.02)),while the two-spot case showed relatively lower model evidence by three orders of magnitude (∼ exp(5.56))than the five-spot case.Then, the model with more than three spots is preferred in the model comparison (Kass & Raftery 1995).The joint posterior for the three-spot case is shown in Figure A.11 using corner (Foreman-Mackey 2016).
Figure 7 visualises the best-fit spot distribution at phase = 0 in the Mollweide projection for the two, three, four and five-spot cases.In the case of two spots, spots appear in the higher latitude (> 60 deg for both spots in the best-fit values) whereas in the case of three or more spots, they are found to be distributed similarly at about similar latitudes around 38 deg, while each spot latitude value has large (> 15 deg) uncertainties.In any case, the optimal solution preferred a pattern where the spot is barely visible when the star is brightest.
We note that our calculation does not take spot overlaps into account (i.e. the overlapping spots are calculated to be twice as dark).The spot latitude is restricted to take only positive values (northern hemisphere) because of the degeneracy with negative values (southern hemisphere) under the stellar inclination of ∼ 90 deg.This overlap could be eliminated if one of the overlapped spots is placed alternately on the southern hemispheres.
The resulting upper limit of spot temperature and spot covering fraction (Figure 6) by spot mapping was consistent with the result of amplitude-spot temperature relations (Figure 5).We will further discuss the results in Section 4.1.

Estimation of transit light source effect
From the derived spot characteristics, we estimated how much the observed transit depths were affected by spots.The degree of TLSE to the transit depth at each wavelength is quantified by the contamination factor   (Rackham et al. 2017) as follows: where  ,obs is the observed transit depth at each wavelength and  ,true is the true transit depth originating from the planetary atmosphere.From Rackham et al. (2018), the   value is determined by the ratio of the flux from the photosphere to the flux from the transit chord.Simply, it is calculated by the equation For each band, the equation can be represented as using the spot covering fraction  spot and spot contrast   at the band , obtained from equations ( 3) and ( 4 assume that there are no spots in the transit chord.This should be a reasonable assumption based on the fact that our observations show no signs of spot crossings (see Appendix A6).
Figure 8 shows the contamination factor (Equation 6) for each wavelength, assuming spots with a temperature of 3100 K on the stellar surface with a temperature of 3200 K as a reasonable assumption from our estimation of Δ from the stellar brightness variability (Sections 3.3 and 3.4).As the  spot can be larger if there are alwaysvisible spots, we calculated the TLSE with different spot covering fractions from 5 % to 30 %.It is clearly seen that TLSE is larger in the optical wavelength range than in the infrared.In the near-infrared wavelength range that is important for the atmosphere characterization, the value of   varies by up to 1%; as the transit depth of K2-25 is ∼ 1%, this is expected to result in a transit depth variation of 100 ppm.
Separately, TLSE can also be calculated from the mapping results.Since Equation 6 represents the ratio of the flux coming from the spotless surface to the current stellar surface, we can divide the maximum value of the flux by the flux coming from the star at each rotation phase derived from the spot mapping.In this approach, the information regarding the spot positions and stellar limb darkening, which are known to potentially influence the estimation of TLSE (Thompson et al. 2024), was incorporated, as it directly utilizes the results of spot mapping.We confirmed that the typical TLSE value derived by this method was comparable to the value calculated by the method described above.This approach also enables the prediction of TLSE variations across different stellar rotational phases.Further details are given in the Appendix A5.

Transit light curve analyses
In this section, we focus on the transit light curves taken by the MuS-CAT2, MuSCAT3, and TESS (summarized in Table 2).MuSCAT2 and MuSCAT3 observed the transits of K2-25 b six times in , , and  band, and four times in  band.TESS observed the six transits.
Therefore, we have 28 transit light curves in total (Figure 10).We analysed these transit light curves to derive the transit depths at each band to evaluate the degree of TLSE.
We jointly fit all 28 light curves using the mean models defined by the transit parameters and noise models defined as Gaussian processes (GP).The parameters, priors, and derived values from the following procedures are summarized in Table 4.
For the mean model of the planetary transit, we used the package batman (Kreidberg 2015), which is based on the analytic formula in Mandel & Agol (2002).As transit parameters, we adopt the transitcentre time  0 , the orbital period  orb , the orbital semi-major axis over the stellar radius / ★ , the impact parameter , and the ratio of the planetary radius to the stellar radius   / ★ .Only   / ★ values were set as different parameters for each band.We assumed a linear ephemeris, as we did not detect any transit timing variations from our light curves.As the previous studies have shown K2-25 b should have an eccentric orbit (Thao et al. 2020;Stefansson et al. 2020), we added two more parameters √  cos  and √  sin , where  and  are the eccentricity and longitude of periastron, respectively.In addition, modified quadratic limb-darkening coefficients ( 1 ,  2 ) were set as parameters for each band (Kipping 2013).We used normal priors (Gaussian priors) for  0 , ,  based on the values in Stefansson et al. (2020), and normal priors for the limb-darkening coefficients based on the theoretical calculation by the package PyLDTk assuming the stellar effective temperature  eff = 3207±58 K, metallicity [Fe/H] = 0.15 ± 0.03 dex, and stellar surface gravity log  = 4.944 ± 0.031 (in cgs).In addition, we set the normal prior for stellar density  ★ to the value derived from Stefansson et al. (2020).Uniform priors were set for the other parameters.
To account for baseline trends in light curves, we used the GP with an approximated Matérn-3/2 kernel: which is implemented in celerite (Foreman-Mackey et al. 2017; Foreman-Mackey 2018).Here,  is the absolute time difference of each observing point, meaning that a data point is correlated with another data point.The hyperparameters  and  determine the correlation of the data points:  corresponds to the amplitude of correlation strength, and  is the timescale of the exponential decay.
was fixed to 0.01 by default.Assuming that the baseline variation is caused by the brightness variation of the star itself or the wavelengthdependent instrument noise, we set hyperparameters for each band.As we do not know any prior information on these hyperparameter values, we set log-uniform priors for the hyperparameters.
As shown in Figures 10 and 11, we jointly optimized 28 transit light curves in total with the mean and noise models.The total number of parameters and hyperparameters were summed up to 31.The parameter estimation was performed using emcee, and the derived values are listed in Table 4.The derived transit depths for each band are also shown in Figure 13.Considering the large uncertainties in  and TESS-band, we do not see significant transit depth dependence on the colour.

Constraints on spot temperature and covering fraction from transit depths
As shown in Figure 8, the TLSE by dark spots makes the observed transit depths larger in shorter-wavelength bands in the optical.The derived transit depths did not show such a trend, although their uncertainties are large.Independently of the method obtained from the monitoring light curves, we investigated the possible spot characteristics only from the transit depths observed by MuSCAT2/3 and TESS.
Assuming the true transit depth is constant for all observed bands (i.e., flat transmission spectra model), The observed transit depths at each band   can be described as where   is a contamination factor for each band and  true is a true transit depth (constant value).  for each band is calculated through the Equation 7. We modelled the   for band  = {, , , , TESS} with three parameters, Δ,  spot , and  true .We set uniform prior U (0, 1500) for Δ and U (0, 1) for  spot . true was restricted to take only positive values.The parameter estimation was done using emcee.
Figure 9 shows the derived posterior distribution using corner (Foreman-Mackey 2016).The 2 upper limit of Δ and  spot values were 1427 K and 0.56, respectively.This means that little restriction could be placed on the spot characteristics from only the wavelength dependence of the transit depths.The posterior of Δ shows a weak bimodal distribution, which is attributed to the fact that the obtained transit depths are consistent with flat considering the uncertainties and does not show a strong wavelength dependence.This is because the wavelength dependence of the spot contrast values is small at both low and high Δ values, as shown in be required, which is unlikely in light of the other observed facts.Nevertheless, models with low-temperature spots spread over a wide area were ruled out.These results were consistent with the results from the monitoring light curves (Figure 5, 6).We also tried to derive Δ,  spot , and  true by jointly fitting the transit depths and modulation amplitudes (Section 3.3).The derived posterior distribution of the parameters was almost identical to the results from the fitting of modulation amplitudes (Figure 5), meaning that the spot characteristics are constrained much better from the multi-band modulation amplitudes than the transit depths, considering the quality of our datasets.

Characteristics of spots on K2-25
From the amplitudes of stellar brightness variations (Sections 3.2 and 3.3), we derived the spot temperature  spot to be Δ < 190 K lower than the assumed photosphere temperature of  phot = 3200 K, and the spot covering fraction  spot to be < 61% with 2 confidence level (Figure 5).While the upper limit of spot covering fraction  spot < 61% is not a strong constraint, the upper limit of spot temperature, Δ spot < 190K, provides sufficient information to scrutinize spot properties and assess the extent of the TLSE, as detailed in Section 4.4.
Furthermore, we explored the potential spot distribution through spot mapping (Section 3.4).Despite limiting the number of spots up to five due to parameter degeneracies, we determined that cases involving three or more spots exhibited a general trend in spot distribution.A further increase in the number of spots could elevate the value of  spot due to certain spots not contributing to brightness variations (i.e., always visible spots).Nonetheless, Δ remained rel-atively stable as the number of spots increased, indicating minimal variation in Δ prediction with different numbers of spots.
Comparing the results from the amplitude (Figure 5) and spot mapping (Figure 6), the spot mapping results showed more strict constraints on both the spot temperature and spot covering fraction in each spot number case.This is because the spot mapping takes into account not only the amplitude of the brightness variation but also the shape of the light curves.Possible sizes and positions of the spots were determined in order to represent the smooth brightness variations of K2-25 monitoring light curves, which resulted in a stronger constraint on the spot temperatures.We note that spot mapping entails various assumptions (e.g., assuming a perfectly circular spot, a small number of spots, etc.).The results may vary if a more flexible spot configuration is allowed (e.g., Luger et al. 2021b).Therefore, the more conservative results obtained from the amplitudes are used in the discussion part of this paper.
In addition, from the wavelength dependence of the transit depths (Section 3.7), we also derived the 2 upper limit of both spot temperature difference (<1184 K) and covering fraction (<0.43) (Figure 9).Despite weak constraints due to transit depth uncertainties, we can rule out the extreme cases with large spot temperature differences and large spot covering fractions.The results are consistent with the results derived by stellar brightness monitoring.Consequently, the resulting spot characteristics were in agreement with the overall features from the modulation amplitudes (Section 3.3), spot mapping (Section 3.4) and from transit depths (Section 3.7).

Effectiveness of simultaneous multi-band observations
In our study, simultaneous multi-band monitoring observations allowed us to achieve exceptional precision to constrain the spot characteristics, especially spot temperature difference Δ.This is mainly because, thanks to the multi-band light curves, degeneracies between spot latitude, spot size and spot temperature are broken to some extent.When the spot is located near the poles, the wavelength dependence of the dimming due to the spot is smaller because of limbdarkening.In this study, the wavelength dependence of the brightness variability amplitudes resulted in taking a low likelihood value of the case with spots near the pole, thereby better defining the spot size and temperature.
We confirmed that our method using multi-band light curves from TESS and ground-based telescopes resulted in better constraint especially on spot temperature estimation, comparing the spot mapping result using five light curves and using only the TESS light curve (details are in Appendix A8).The result of spot mapping effectively depends on the TESS light curve with higher time cadence and photometric precision, but multi-band data from ground-based telescopes are useful in putting further constraints on the solutions.
There are a few studies of such multi-band (but not simultaneous) monitoring observations conducted for M-dwarfs.Morris et al. (2018) analysed the monitoring light curves of an M8 star TRAPPIST-1, by Kepler and Spitzer, and constrained the spot temperature to be extremely hot ( spot ≳ 5300 ± 200 K).Miyakawa et al. (2021) analysed the monitoring light curves of four M dwarfs in the Hyades cluster in the Kepler-band and , , -band, and derived their spot temperatures within 1 uncertainties of ≲ 700 K.
Our observations constrained the  spot values more precisely than these studies.This is because our observations obtained multi-band light curves in optical, which are more sensitive to spot temperature than in infrared.Also, simultaneity was important, as the spot characteristics may change in a timescale of the spot emergence and decay (Appendix A5).Our observations showed that simultaneous multi-band monitoring using a combination of ground-based and space-based telescopes will be effective in determining the "current" spot characteristics of the host star, for example, in future planetary atmosphere research.

Indication to the spot-photosphere temperature relations
Figure 12 shows the relation of Δ spot =  phot − spot and  phot values from this work in comparison with the ones from recent literature for cool dwarfs derived by various methods (Berdyugina 2005;Frasca et al. 2009;Miyakawa et al. 2021;Almenara et al. 2022;Waalkes et al. 2023).
There were several studies to derive the empirical relations of Δ and  phot with polynomial functions (Berdyugina 2005;Maehara et al. 2017;Herbst et al. 2021).Berdyugina (2005) reports that Δ increases as the value of  phot increases for both the dwarfs and giants.Herbst et al. (2021) updates the empirical relations by incorporating the latest samples, mainly for stars with effective temperatures of  eff > 4500 K.After Berdyugina (2005), a few observations slightly increased the number of samples with  eff < 4000 K (Frasca et al. 2009;Miyakawa et al. 2021;Almenara et al. 2022;Waalkes et al. 2023).The results of those studies and our result suggest that the overall trend almost corresponds to the empirical relationship, but a trend towards smaller Δ with lower photosphere temperature, except for a young pre-main sequence M dwarf AU Mic (Waalkes et al. 2023).
Whereas the number of samples is still limited, these results suggest that the relations might no longer be suitable for M-dwarfs, especially with  eff ≲ 3500 K.It is important to note that in the prediction of the TLSE for each planetary system,  spot values are often assumed based on the empirical formula mentioned above (e.g., Rackham et al. 2018).In such cases, there is a possibility of overestimating Δ values (i.e., underestimating  spot values).It is possible to update the relationship for M-dwarfs by adding samples with high-precision measurements of Δ as future issues.In addition, the accumulation of the measurements on the  spot for low-temperature stars would be  useful in theoretical studies of the starspot properties.In Panja et al. (2020), the temperature dependence of the H − opacity explains the trend of Δ spot and  phot , and this reproduces the trend derived by Berdyugina (2005).However, the calculation has not yet been done for stars with temperatures lower than 3700 K (i.e., M0-type stars).It is important to increase observational data for the evaluation of these theoretical explanations.

Indication to the K2-25b atmosphere retrievals
In Figure 13, we show the derived transit depths of K2-25 b in five bands (, , , , and TESS-band) along the wavelength.We overplotted them with the previous transit depths derived by Thao et al. (2020), which used the data taken with Kepler, Sinistro -band, MEarth, Spitzer channel 1, and Spitzer channel 2. We found that the transit depths from our observations have smaller overall values than the ones from Thao et al. (2020), despite the wavelength ranges being similar.However, given the large errors especially in the , , and TESS bands, all transit depths are consistent within a 2 range.We did not ascertain the trend of deeper transit in the visible compared to the infrared, which was reported in Thao et al. (2020).This result might be due to some systematics from different observations, different analytical methods, and different stellar activity levels.
We also overlaid the calculated transmission spectrum models assuming TLSE and planetary atmosphere compositions.The atmosphere model with 100× solar metallicity is calculated using the python package PLATON (Zhang et al. 2019) assuming the chemical equilibrium, from  ★ = 0.2932 ⊙ (Thao et al. 2020),   /  = 0.108 (this work),   = 24.5⊕ and the equilibrium temperature  eq = 494 K (Stefansson et al. 2020).The TLSE was calculated assuming the spots with  spot = 3100 K and  spot = 15 %, which are reasonable values from our monitoring observations.While changing the y-axis offset as a parameter, we fit the models to data and calculated the chi-squared and reduced chi-squared values for four models: a flat atmosphere spectrum with and without TLSE, and a 100×-Solar-metalicity atmosphere spectrum with and without  5.We do not see a significant difference in chi-squared values for all models, in both cases using the transit depths only from this work and from this work combined with Thao et al. (2020).
While Thao et al. (2020) rejected the clear 100× solar metalicity model, our analysis shows that this possibility cannot be dismissed.This is not due to differences in stellar surface models, but mainly because the atmospheric model features were reduced due to the larger planetary mass estimates, with reference to the results of Stefansson et al. (2020).For the combined dataset from our work and Thao et al. (2020), the chi-squared value is slightly smaller when the spot effects are taken into account, for both the flat and 100×-Solar-metalicity atmosphere cases, but this is also not a statistically significant difference.
Also, in Thao et al. (2020), they pointed out that spot covering fractions of 22% to 36% are too large and unlikely given that the stellar variability is 0.5% to 2%, but our analysis does not rule out such a large spot covering fraction, considering the existence of spots that do not contribute to the stellar brightness variability.
It should be noted that the spectrum from Thao et al. ( 2020) is stacked, i.e. consisting of individual datasets taken between 2015 and 2017, while our spectrum was taken in 2021.Given that our observations show that the spot distribution in K2-25 varies on scales of hundreds of days (see Appendix A3), it would be not good practice to use all datasets to compare with models.We therefore stop here for a detailed comparison using both results.To make a real comparison for the atmosphere constraint of K2-25b, we need to simultaneously obtain spectra over a wide range of wavelengths with high-precision spectroscopy, such as by JWST.
As in Figure 8, the TLSE is not very significant in the JWST/NIRSPEC wavelength range (600 -5000 nm) compared to the optical.However, for example, when comparing the contamination factor at wavelength ∼ 1500 nm and 4000 nm, there are some cases where the   differs by ∼ 1 %.This means that even if the "true" transit depth is the same (i.e., flat transmission spectra), the observed depth differs by 100 ppm relative to each other.
Given the photometric precision of JWST's observations, this effect cannot be ignored.We have confirmed from the simulation of JWST spectra using software PandExo (Batalha et al. 2017), that when K2-25 is observed with PRISM at JWST/NIRSPEC, a precision of < 100 ppm can be achieved in a single transit, although this depends on the wavelength binning.This indicates that stellar surface models need to be examined when inferring atmospheric models for K2-25 b from transmission spectroscopy by JWST.In that case, our knowledge of the spot temperature and the covering fraction of K2-25 from the monitoring observations is a good prior to investigating the stellar surface models.Of course, ideally, the monitoring observations should be conducted at the same time as the transmission spectroscopy by JWST.

Other stellar activity features on K2-25: indication for the future atmosphere research
We also investigated some other stellar activity signals on K2-25 from the analyses of multi-band light curves (see Appendix A3, A5 to A7 for details).We found there is a long-term variability of the modulation amplitude in a timescale of a few ×100 days.This might be explained by the timescale of the change in the spot distribution.This means it is important to obtain the concurrent status of stellar activity when the transmission spectroscopy is conducted.
From the analysis of flare frequency, we confirmed three flares in K2-25 with the energy at ∼ 10 33 erg, in a 21.76 days-long TESS  light curve.This is a typical flare frequency for a young mid-M star (Günther et al. 2020;Murray et al. 2022), although the number of flare detections was limited because of the detection limit.Smaller flares should happen more frequently, which could be the dominant source of baseline fluctuation in transit light curves (especially in the bluecolour band).When performing high-precision transit observations using space telescopes, at least multiple observations should be made, since flares may occur during the observation and cause the baseline modulation in the transit light curves.

CONCLUSION
In this study, we examined the stellar brightness variability of the young mid-M-dwarf K2-25 and analyzed the transits of its sub-Neptune companion K2-25 b using data from TESS and groundbased telescopes.Our aim was to investigate the characteristics of starspots within the system and their impact on transit observations.The significance of our study stems from the utilization of (i,ii) "multi-band" and (iii) "simultaneous" observations, which enabled precise constraints on spot characteristics.
(i) Multi-band monitoring of stellar brightness provided robust constraints on spot temperature.Using monitoring light curves in four different bands (utilizing five different instruments/filters), we revealed the wavelength dependence of the amplitudes of stellar brightness modulation.From there, we found the spot temperature difference from the photosphere temperature of K2-25 to be < 190 K and the spot covering fraction to be < 61% with 2 confidence level.These results were consistent with the one derived from spot mapping, assuming two, three, four or five spot cases (Section 4.1).The derived spot temperature difference was smaller than the extrapolation of the empirical relation with the stellar photosphere temperature.This suggests that such an empirical relation is no longer suitable for M-dwarfs (Section 4.3).(ii) Multi-band transit photometry was also obtained to constrain stellar surface models through the calculation of the TLSE.Although the derived transit depths by MuSCAT series and TESS had uncertainties too large to precisely constrain the spot temperature and covering fraction, we confirmed that the results were consistent with the spot characteristics constrained by the monitoring observations (Section 3.7).We calculated that the expected TLSE would distort the transmission spectrum of the system with the amplitude of ∼ 100 ppm in the JWST/NIRSPEC wavelength range, which is detectable with JWST (Section 4.4).
(iii) This study was the first trial of taking stellar monitoring and transit observations simultaneously (Section 4.2).From the longbaseline light curves, we found that the spot distributions should change on timescales of hundreds of days.For more precise transit observations, it is important to know the spot distribution at each observing epoch.
There is room for further improvements to our methods, such as improving observation precision and refining stellar surface models for more precise derivation of stellar parameters from light curves.For example, faculae were not considered in our analysis because of the difficulty in quantifying their contributions to the light curves: faculae are known to become brighter towards the star limb (Norris et al. 2023).Although faculae are expected to contribute less to the brightness variability than spots because faculae are more widely distributed on the stellar surface (Johnson et al. 2021), how to handle the contribution from faculae to the light curve is a future issue.
Nevertheless, our investigation of the spot characteristics of an M-dwarf, using the best current instruments and methods, is an important step forward.Employing this method, along with further improvements as mentioned above, may be crucial for more accurate observations of planetary atmospheres with space telescopes, and for further studies of the stellar activity of M-dwarfs.

DATA AVAILABILITY
The TESS data can be accessed through the MAST (Mikulski Archive for Space Telescopes) portal at https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html.The ZTF data can be accessed through https://irsa.ipac.caltech.edu/Missions/ztf.html.The data taken by Sinistro, MuSCAT2 and 3 are available upon request to the authors.The SVO Filter Profile Service and theoretical stellar spectra can be obtained from the SVO Theoretical Model Services http://svo2.cab.inta-csic.es/theory/main/.
This work has utilized these public software: LCOGT BANZAI pipeline (McCully et al. 2018)  A3 Analysis on the long-term change of the variability amplitude Kain et al. (2020) revealed that the modulation amplitude of K2-25 differs from year to year from the analyses of light curves observed with the MEarth Observatories (Irwin et al. 2015) over a three-year period.To see the recent long-term variation of the modulation amplitude, we divided the ZTF light curves into 8 different segments, each of which has a 200-day period and derived modulation amplitudes by the GP method.This time, the hyperparameters were fixed to the derived value described in Table A.1 for each ZTF band.The shape of the brightness modulation is changing in a timescale of ∼ 200 days, as shown in Figure A.3.The amplitude also varied at each time, with a relatively large amplitude during the period we observed with other instruments (MJD 59400 -59800), and almost no visible variations during the previous period (MJD 59000 -59400).In addition, the shape of the variation is consistent between   and   -band, with   -band always showing a greater amplitude.
This means the spot distribution has changed in a timescale of a few ×100 days, and it is difficult to constrain spot characteristics by the combination of different observations at the different observing epochs.

A4 Possibility of bright spots
In this section, we deal with the case of bright spots, which was not considered in the derivation of spot temperature and covering fraction from the brightness modulation amplitudes (Section 3.3).Bright spots can create similar modulations in light curves as described by Equation 3, although the bright spots are most visible when the star is brightest.We obtained a BT-Settl model with spot temperatures ranging from 3200 K to 5000 K in 100 K increments, which was interpolated for use.−11 K and  spot = 12 ± 2 % with 1 uncertainties.The TLSE caused by those bright spots results in smaller transit depths in shorter wavelength in the optical.We calculated the chi-squared values of the calculated transit depths to observed transit depths as done in Section 4.4, adopting reasonable values for the bright spots case ( spot = 3300 K and  spot = 15 %).Table A.3 presents similar chi-squared values for all cases, which do not strongly suggest the existence of bright spots, considering the large errors in observed transit depths.
Here, bright spots refer to regions that have a higher temperature than the photosphere, and their physical characteristics are not con- sidered.These could be faculae, but as explained in the Conclusion, faculae exhibit complex behaviour in light curves and are challenging to implement in spot mapping.This remains a future task, including how such regions, which slightly warmer than their surroundings, can be physically explained.

A5 Analysis on the transit depth variability
We analysed the transit depth variability to see if the TLSE causes the variability at the different stellar rotational phases.Transit analysis was performed following the method described in Section 3.6.However, in this analysis, the transit depth was set as a separate parameter for each observation night, while the other transit parameters were common.Fitting was done for the four MuS-CAT bands simultaneously.To reduce the computational cost, the GP hyperparameters were fixed with the best-fit values from the joint-fit results (Table 4).Therefore the total number of parameters was 36:  0 , , √  cos , √  sin , / ★ , , and (  / ★ )  for each light curve  = {1, 2, • • • 22}, and ( 1 ,  2 ) for each band  = {, , , }.Note that, since TESS transit data are not very high precision, we did not use them to examine transit depth variability in the first place.
Figure A.5 shows the derived transit depth variations in each band, along with the observing time and corresponding stellar rotational phase.The stellar rotational phase at each transit timing was calculated using the values of  rot = 1.87708 day, described in Section 3.1.
At the same time, we calculated the degree of TLSE (contamination factor)   at each band  per stellar rotation phase.In principle, the derivation of   (Equation 6) means the ratio of the 'flux from the spotless stellar surface' to the 'flux from the current stellar surface'.This means that we can calculate the   at each phase by dividing each phase-folded light curve by the maximum of the light curve.The dashed lines in

A7 Analysis on the flares from TESS light curve
From visual inspection, we found there are several flare-like signals in the TESS light curve.Frequent flares would not only lead to errors in transit depth measurements but would be another source of transit depth variations.Therefore, we derived the flare frequency (how often flares of each energy occur) of the system, following the five steps; (1) detrending the light curves, (2) searching for flare candidates, (3) fitting flare models to the candidate light curve, (4) deriving total energy for each flare, and (5) deriving flare frequency.
First, we detrended the light curve.Assuming that the flares occur almost independently of the stellar rotation phase, models fitted to light curves folded in the stellar rotation phase can be regarded as unaffected by flares.Therefore, we detrended the TESS light curve by the phase-folded GP model described in Section 3.2 and unfolded it with time as the horizontal axis.
Second, the flare candidates were searched using the criteria described in Chang et al. (2015): "There are three or more consecutive points where the flux is more than 3 off the overall median value".As a result, three flare candidates were detected.We also visually checked all the other data points that excess 3 from the median and did not confirm flare-like signals other than these three.We cut out a light curve at 0.02 days before and 0.08 days after the first rise of each flare-like signal.To model the flare light curves, we used the empirical flare template Llamaradas Estelares (Mendoza et al. 2022), which is an updated model from Davenport et al. (2014).The flare shape in the light curves  flare () is characterized by three parameters: the time of flare peak  peak , the time-scale of decay FWHM, and the peak amplitude of the flux  amp .Independently for each flare candidate, parameter estimation was performed by emcee to maximise the likelihood of the flare model to the observed flare light curve.The emcee chain was well-converged, indicating that the probability that these are typical flares is high.
Therefore, the equivalent duration (ED) of the flare and flare energy

Figure 1 .
Figure 1.Results of period analysis for Sinistro-, ZTF-, ZTF-, TESS, and Sinistro- bands.The vertical lines for each panel show the detected period (1.88 days).The strong peak at ∼ 2.14 days show the one-day alias of the signal.The horizontal lines indicate the false alarm probability at 0.1%.All light curves show clear periodicity.

Figure 2 .
Figure 2. Derived modulation amplitude for the five light curves by sinecurve method (Section 3.2).The bands are listed from left to right in order of effective wavelength, from shortest to longest.A clear trend is seen: the longer the wavelength, the smaller the modulation amplitude.The values are summarized in Table A.2.

Figure 3 .
Figure 3.The phase-folded light curves and fitted model curves.The solid line shows the best-fit curves and the shaded region shows a 1 deviation of model uncertainties, derived from the models drawn by the converged MCMC chains.

Figure 4 .
Figure 4. Relation of calculated spot contrast   to Δ for each observed band  = {,  , , , TESS}.Δ represents the difference of spot temperature from the stellar surface temperature fixed at 3200 K.The wavelength dependence of the value of   is greatest around Δ = 500 K.

Figure 5 .
Figure 5. Posterior distribution of the spot temperature difference Δ =  phot −  spot and the spot covering fraction  spot .The red solid lines show the median values and the dashed lines show the upper limit with the 2 level.

Figure 6 .
Figure 6.Posterior distribution of Δ and  spot , with two, three, four, and five-spot cases shown in red, orange, green and blue.The dashed lines indicate the 2 upper limit of Δ and  spot at each case.

Figure 7 .Figure 8 .
Figure7.Spot map with optimum spot parameters for the two, three, four, and five-spot cases shown in red, orange, green and blue, respectively.The spots are located at high altitudes in the two-spot case.In cases with a larger number of spots, the spots show a similar distribution in the mid-latitude zone.

Figure 9 .
Figure9shows the derived posterior distribution using corner (Foreman-Mackey 2016).The 2 upper limit of Δ and  spot values were 1427 K and 0.56, respectively.This means that little restriction could be placed on the spot characteristics from only the wavelength dependence of the transit depths.The posterior of Δ shows a weak bimodal distribution, which is attributed to the fact that the obtained transit depths are consistent with flat considering the uncertainties and does not show a strong wavelength dependence.This is because the wavelength dependence of the spot contrast values is small at both low and high Δ values, as shown in Figure4.However, for large values of Δ, either very small  spot values or  true values will

Figure 10 .Figure 11 .
Figure 10.Transit light curves of K2-25 b taken with MuSCAT2 & 3 in , , , -band (from top to bottom) on each observing night. T number is shown next to each light curve.In each subplot, observed data points are shown in grey points, and best-fit transit + baseline models derived by GP regression are indicated in coloured lines.The shaded regions indicate the 1 error of the GP model.

Figure 12 .
Figure 12.Relations between  phot and temperature difference Δ spot =  phot −  spot for the dwarfs between 3000 K to 5000 K.The empirical relation and its uncertainties from Equation 4 in Herbst et al. (2021) (derived from dwarfs and giants in Berdyugina (2005)) are indicated in the blue line and filled regions.Note that data points without error bars are those for which the errors are not estimated.For the stars with  phot ≲ 3500 K, Δ tend to show smaller values compared to the empirical relation.

Figure 13 .
Figure 13.Generated transmission spectrum of K2-25b overplotted with the observed transit depths in multi-band.The black dots with errorbars indicate the transit depths by this work in , , , TESS and -band from left to right, and grey dots with errorbars indicate the values from Thao et al. (2020).The observed transit depths are compared with the planetary atmosphere models.Green/blue lines are the flat atmosphere model and orange/pink lines are the forward model generated by PLATON, with 100× solar metalicity atmosphere, with (solid) and without (dashed) the effect of spots (TLSE).The y-axis offset of each model was optimized for the data obtained in this work.The top figure is the close-up version of the bottom figure.

Figure A. 2 .
Figure A.2.The phase-folded light curves and fitted model curves by the GP method.The solid line shows the best-fit curves and the shaded region shows the 1  uncertainties of the models.The white noise are taken into account by adding white noise kernels.
Figure A.4 shows the resulting posterior distribution of spot characteristics from the modulation amplitudes.The derived values are |Δ spot | = 78 +16

Figure A. 4 .
Figure A.4.Similar to Figure 5, but for the case of bright spots.The red lines represent the median values of the posterior distributions.
shows the normalized TESS light curves, with and without detrending, along with the times of three detected flare candidates.

Figure A. 5 .
Figure A.5. Derived transit depths according to the stellar rotation phase for , , , and -band from the top to the bottom.The grey lines and regions show the transit depths and their uncertainties derived by joint fit.The black dashed lines show the predicted transit depth variation from the TLSE.The derived uncertainties in the transit depths are larger than the predicted transit depth variation.

Figure A. 8 .
Figure A.7. Normalized PDCSAP light curve of K2-25 from TESS Sector 44 (top) and the detrended light curve by phase-folded GP model (bottom).Data points whose normalized flux values are more than 3 away from the median value of the entire detrended curve are marked with orange crosses.The flare candidates are indicated by the grey dotted lines.

Figure A. 9 .
Figure A.9. Raw light curves by MuSCAT2/3 before removing outliers, in , , , -band from top to bottom at each panel, at each night.The coloured dots indicate the data points used for transit analysis, and cross marks indicate the removed data points.The values at the upper-right corner of each panel indicate the  T values.There are several flare-like signals detected and removed as outliers.The details of outlier removal are described in Section 2.4.

Table 2 .
Summary of the observational data used in the analyses, where  T indicates the number of transits counted from the reference transit-centre time  0 .

Table 3 .
Optimum values for   ,   , and   in Equation 1 to fit the phasefolded TESS light curve.

Table 4 .
Priors and derived parameters of the transit analyses.U(A, B) indicates the uniform prior from A to B and N(X, Y 2 ) indicates the normal prior with the mean value X and standard deviation Y.
For parameters with normal prior distributions, the mean and variance are adopted fromStefansson et al. (2020) from [ 1 ,  2 ] values calculated by PyLDTk.

Table 5 .
Calculated chi-squared values and reduced chi-squared values for each atmosphere model with and without spot effect, assuming  spot = 3100 K and  spot = 15 %.The values are calculated for (1) only using MuSCATs and TESS data from this work and (2) combined data from this work and from Thao et al. (2020) (T20).