The effect of spots on the luminosity spread of the Pleiades

Cool spots on the surface of magnetically-active stars modulate their observed brightnesses and temperatures, thereby affecting the stellar locus on the H-R diagram. Recent high precision space-based photometric surveys reveal the rotational modulation from spots on stars in young clusters, including K2 monitoring of the 125-Myr-old Pleiades cluster. However, lightcurves reveal only the asymmetries in the visible spot distributions rather than the total sizes of spots on stellar surfaces, which causes a discrepancy between the spot coverage measured by photometric and spectroscopic observations. In this paper, we simulate photometric variability introduced by randomly-distributed starspots on a 125-Myr-old coeval cluster. Our simulation results show that randomly distributed small spots on the stellar surface would explain this discrepancy that the photometric observations only reveal 10% to 40% of the spot coverage measured by spectra. The colors and luminosities of photospheres are modeled for a range of photospheric temperature, spot coverage, and spot temperature. The colors and luminosities of a simulated population are then compared to the luminosity spread of Pleiades members, excluding the 25% of candidates that are identified as non-members with Gaia DR2 astrometry. The observed luminosities of Pleiades members have a standard deviation of 0.05 dex, which could be entirely explained by spots with a star-to-star standard deviation of spot coverage of 10%, but with an average coverage area that is not well constrained.


Introduction
Starspots are the visible manifestation of internal magnetic activity (see reviews by Schrijver & Zwaan 2000;Strassmeier 2009). The magnetic dynamos on the Sun and other stars are produced by the rotation of convective plasma (Browning et al. 2010;Charbonneau 2014). Single stars rotate fastest when they are young, generating strong magnetic activity (see review by Bouvier et al. 2014). While solar magnetic activity has a negligible effect on the total irradiance (Balmaceda et al. 2009) and flux transport of the Sun, at young ages the internal magnetic fields can change the structure of the star and the efficiency of convection (Somers & Pinsonneault 2015;Feiden 2016;MacDonald & Mullan 2017). Starspots, one of the consequences of this magnetic activity, also complicates the measurement of stellar properties by introducing additional temperature components.
Since magnetically active stars are rarely resolved (see, e.g., Roettenbacher et al. 2016), starspots are usually detected with photometric monitoring (e.g., Hall 1972;Strassmeier et al. 1997;Meibom et al. 2009). The recent space-based CoRoT and Kepler missions provide sensitive light curves to hunt for exoplanets (Baglin 2003;Borucki et al. 2010), and can also be used to measure stellar rotation periods for analyzing angular momentum evolution (e.g., Matt et al. 2015). The rotation amplitudes of the light curves are measured as consequences of starspots (Savanov & Dmitrienko 2017). In addition, the stellar surfaces are mapped by light curve inversion techniques that infer the spot geometry and evolution through the morphologies of light curves (e.g., Savanov & Dmitrienko 2012;Roettenbacher et al. 2013). The K2 light curves of stars in the 125 Myr old Pleiades cluster typically have amplitudes of 1%-10% (Rebull et al. 2016a), and less than the 30%-50% amplitudes are often seen in at least some 1-10 Myr stars (e.g., Grankin et al. 2008;Alencar et al. 2010;Lanza et al. 2016). However, the total spot coverage only from light curve amplitudes may be underestimated because symmetric morphologies, such as polar spots, will not cause variability (see discussions in Rebull et al. 2016b;Rackham et al. 2018).
In contrast to light curves, spectroscopic techniques can provide estimates of the total spot coverage by comparing (usually molecular) features of a star with those from an unspotted template (e.g., Petrov et al. 1994;Neff et al. 1995;O'Neal et al. 1996;Fang et al. 2016;Gully-Santiago et al. 2017). The magnetic fields themselves can be measured either through Zeeman broadening, which yields an averaged magnetic field strength over the visible stellar surface (e.g., Johns-Krull et al. 2004;Lavail et al. 2017), or using polarized light, which when combined with high-resolution spectroscopy yields maps of the largest magnetic structures (Zeeman Doppler Imaging; see the review by Donati & Landstreet 2009).
These detection methods all leverage the effect that spots have on the emission from the star. While many studies have investigated the subsequent limitations on our ability to measure the presence of planets (e.g., Desort et al. 2007;Reiners & Christensen 2010) and exoplanet atmospheres (Rackham et al. 2018), the spots also interfere with our ability to measure the properties of spotted young stars. Spots alter the evolution of pre-main-sequence stars (Somers & Pinsonneault 2015), require radius inflation relative to predictions from standard models to emit the same amount of energy (e.g., Jackson & Jeffries 2013;Somers et al. 2017), and introduce cooler components into the emission from the star (Gully-Santiago et al. 2017). Differences in spot properties among stars may introduce a spread in the observed properties of coeval cluster members, which could help to explain the luminosity spreads seen in every star-forming region (see reviews by Preibisch 2012;Soderblom et al. 2014).
In this work, we investigate the spot properties of stars in the Pleiades cluster (Subaru in Japanese and Mǎo in Chinese), the closest young open cluster with a low foreground extinction (Mermilliod 1981). The Pleiades has an age of 125 Myr, as estimated by the lithium depletion boundary (e.g., Stauffer et al. 1998), which places this cluster near the peak of the angular velocity evolution of low-mass stars (see review by Bouvier et al. 2014). In Section 2, we describe the differences between photometric and spectroscopic measurements of spots on the Pleiades. In Section 3, we simulate light curves for different spot properties and demonstrate that the spectroscopic and photometric differences may be explained if these stars have many small spots. In Section 4, we then investigate how spots affect the loci of low-mass young stars in colormagnitude diagrams and discuss observational biases introduced when starspots are not considered. In Section 5, we explain the observed luminosity spread in the Pleiades cluster by the appearance of cool spots. We then summarize our results in Section 6.

Observations of Spots on Stars in the Pleiades
As the nearest young open cluster, the Pleiades is a cornerstone for studies of the early evolution of young stars and their rotation. The cluster includes ∼2100 known members (Stauffer et al. 2007;Lodieu et al. 2012;Bouy et al. 2015) with a mass function described by a log-normal distribution and a mean characteristic mass of 0.2 M e (Hambly et al. 1999;Moraux et al. 2003Moraux et al. , 2004Deacon & Hambly 2004;Lodieu et al. 2012). The distances of the Pleiades members have an average of 137 pc from Gaia DR2 parallaxes with a 4 pc standard deviation that includes measurement uncertainties and a real depth (Gaia Collaboration 2018). Spots on Pleiades stars have been extensively detected through photometric monitoring (e.g., Stauffer et al. 1987;van Leeuwen et al. 1987;Hartman et al. 2010;Covey et al. 2016). Redder V−K s colors 5 are observed on fast-rotating low-mass members, suggesting cool spots or stronger magnetic activities on fast rotators. In this section, we first revisit the photometric and spectroscopic observations for the Pleiades members, and then compare these results to obtain a geometric view of the starspots on Pleiades members.
The photometric light curves for the Pleiades members analyzed here were obtained in the K2 extension (Howell et al. 2014) of the Kepler Space Telescope. The 72 day long K2 monitoring campaign (C04, 2015 February 7-April 23) included 826 candidate Pleiades members, of which 92% have accurate periods measured from spot modulation (Rebull et al. 2016a). 6 For low-mass stars, the photometric amplitude, ΔF phot , defined as the 10th to 90th percentile range of the K2 light curve, is typically 1%-10% of the average emission. Estimating the spot coverage from photometric light curves alone captures solely asymmetric structures, resulting in significant underestimates of the starspot area (Rackham et al. 2018). In this work, we adopt the photometric amplitudes of Pleiades members measured by Rebull et al. (2016a) from K2 light curves obtained in Cycle 4.
The spot contributions for 304 Pleiades members were measured by Fang et al. (2016) through low-resolution (R∼1000) spectra spanning 3700-9000 Å that were obtained with LAMOST (Large sky Area Multi-Object fiber Spectroscopic Telescope; Zhao et al. 2012). The spot temperature and covering area were measured based on empirical relationships between TiO absorption band ratios. Fang et al. (2016) fit the observed TiO features by a warmer photospheric temperature, derived from V−I colors, and a cooler spot temperature left as a free parameter. This spot measurement includes a degeneracy between spot coverage and spot temperature, so Fang et al. (2016) adopt the minimum spot size that can explain the TiO feature depths, hence a lower limit of spot temperature. The LAMOST-based estimates of f spot are more uncertain when the stellar effective temperature (T eff ) is smaller; for example Fang et al. (2016) reports typical measurements of 0%-50%, with comparable uncertainties of 10% for T eff >3800 K and >15% for T eff <3800 K. The measured spot parameters from Fang et al. (2016) is provided through private communication. The photometric observations of spots in the Pleiades are combined in this paper. The spectroscopic amplitude (ΔF spec ) is defined here, based on the spot and stellar parameters from the LAMOST spectra (Fang et al. 2016), as where f spot represents the spot coverage. F spot and F phot are the fluxes from the spot and photospheric regions in a unit area through the Kepler band, and both are generated from BT-Settl models with solar metallicity (see Section 3.1 for more information). The spectroscopic amplitude ΔF spec corresponds to the change in the absolute brightness introduced by the spots at the time of the observation. The definition of ΔF spec is the slightly mismatched ΔF phot . In the extreme scenario where the photometric amplitude is described by a single spot on the visible hemisphere, then ΔF spec would equal 1.25×ΔF phot for a given spot filling factor and temperature. The ratio between ΔF phot and ΔF spec offers a simplified geometrical assessment of the level of symmetry of the spots. When the star is covered by spots with a high degree of longitudinal symmetry, ΔF phot would be much smaller than ΔF spec . For instance, polar spots, circumpolar spots on inclined stars, Jupiter-like bands, and a uniform distribution of small spots would all induce negligible temporal photometric modulation despite the coverage of the stellar surface. On the other hand, in some cases, the spectral amplitude is equal to the photometric amplitude, which indicates that the star may have a single large spot visible on one hemisphere. Figure 1 compares the spot modulation amplitudes of 113 Pleiades sources that have both a photometric amplitude measured from K2 (Rebull et al. 2016a(Rebull et al. , 2016b and spot parameters estimated from LAMOST spectra (Fang et al. 2016). Most stars sit closer to the symmetric or "circumpolar spot" regime than to the "single equatorial spot" regime. The large uncertainty in LAMOSTbased f spot estimation results in an±8% error bar on ΔF spec in Figure 1, especially since Fang et al. (2016) used a spot temperature that leads to the minimum spot coverage. However, when taken in aggregate, the spectroscopic amplitudes are much larger than the photometric amplitudes.

Simulating Spotted Photospheres
In the previous section, we compared the results from Fang et al. (2016) and Rebull et al. (2016a) to establish that spots cover a large fraction of the surface of young stars yet produce light curves with small amplitudes. In this section, we simulate broadband light curves of spotted young stars for a variety of spot distributions and parameters to help interpret this difference between the photometric and spectroscopic measurements of spots.
The simulations of inhomogeneous stellar photospheres in this work follow similar methodologies of stellar spots by Desort et al. (2007) and Rackham et al. (2018) to understand their effects on exoplanet detection and characterization. Our Monte Carlo simulations first generate coeval 125 Myr old stellar clusters, each containing 2000 stars with masses evenly spanning from 0.08 to 1.32 solar mass, with stellar effective temperatures given by stellar evolutionary models . Then, cool spots with certain temperatures are added on each stellar surface, occupying a range of regions and geometries (more details of spot configurations are introduced in Section 3.2). Finally, light curves in photometric bands (B, V, R, I, K s ) are integrated from the synthetic BT-Settl models (CIFIST2011 2015 Allard 2014) and stellar rotations.
We use these simulations to investigate the impact of a spread in spot coverage on the luminosity and color spreads on coeval star clusters. Then, the spreads in spot areal surface coverage are ultimately constrained based on the observed spread in the pre-main-sequence H-R diagram. Hence, the observed luminosity and color spreads may be introduced by existing spots rather than by any bona fide age spread. In this context, our starspot simulations inform the analysis of H-R diagrams by helping to reveal the different possible configurations that could simultaneously explain the photometric and spectroscopic detection of spots.

Stellar Samples
The Pleiades cluster is simulated with 2000 stars at an age of 125 Myr, with temperatures and luminosities adopted from the evolutionary tracks . The stellar masses are distributed between 0.08 and 1.32 M e 7 to span the range from the approximate detection limit of K2 to the approximate temperature above which spot modulation mixes with stellar pulsations (Balona et al. 2015). The inclination angles (i, j) of the stellar axes are randomly assigned between 0°<i<180°a nd −90°<j<90°, where i is the inclination of the rotation axis toward the observer, and j is the rotational angle perpendicular to i and on the horizontal plane that contains the direction of line of sight. When applying the stellar inclination, we rotate the star under the rotational coordinate transformation by i first, then by j.
Starspots are then introduced to the stellar surface by displacing warm photosphere with cooler (albeit nonzero) emitting regions, making the star appear redder and fainter. The radius is kept fixed to the Baraffe et al. (2015) models, so the total luminosity of the star is less than a spot-free stellar surface. Diminishing the emergent luminosity creates an artificial mismatch between the internal energy production from contraction, fusion, and the total luminosity emitted through the stellar surface. This mismatch would ordinarily lead to radius inflation, which has recently been observed in the Pleiades cluster (Somers et al. 2017;Jackson et al. 2018) and in rapidly rotating M dwarfs (Kesseli et al. 2018), with stars outsizing predictions by about 10%-18%. We assume for now that radius inflation would translate to isochrones but would not introduce additional scatter on the H-R diagram. Faculae, accurate limb-darkening effect, and warm chromospheric activity are not considered in these simulations.
The emission from both cool spots and the photosphere are calculated from the BT-Settl models (version CIFIST2011_2015, Allard 2014 with solar abundances from Caffau et al. 2011). The synthetic spectra are generated based on the T phot , T spot , and gravity from the Baraffe et al. (2015) pre-main-sequence evolution grid. The stellar emission is hence defined by combining two temperature components modulated through the spot filling factor. The flux in the Johnson's B, V, I, and 2MASS-K s photometric bands in any single epoch is calculated by integrating the synthetic spectra over generic filter transmission curves (from SVO service; Rodrigo et al. 2012) in flux space (W/(m 2 μm)).

Spot Configurations
In this section, we introduce how cool spots are simulated on stars. First, the stellar surfaces are constructed by twodimensional maps with latitude θ ä(−90, 90) and longitude fä(0, 180) with a resolution of one degree on each angular axis. Then, spots are introduced with ranging size and temperature. The stellar rotation is simulated by moving pixels step by step along the longitude axis on the (θ, f) map. The (θ, f) map is then converted to an observed surface map (θ obs , f obs ) via rotational coordinate transformation according to inclination angles. For an inclination of 90°, each point on the stellar surface is visible for half of the revolution, while for an inclination of 0°, only one hemisphere is ever visible and the light curve would be constant.
The simulated spots on the (θ, f) plane are controlled by five variables: the number of spots on the stellar surface (N spot ), and the size (S i ), the temperature (T spot,i ), as well as location (L spot (θ i , f i ) ) of each individual spot. The total spot filling factor f spot on Cartesian coordinates constrains the N spot and S i , reducing the number of independent variables by one, where S å represents the surface area of the star. The range of spot coverages in the simulations is f spot ä(30-50)% for latetype stars (T phot <3800 K) and (1-40)% for warmer stars, notionally consistent with the measurements from the TiO band absorption seen with LAMOST (Fang et al. 2016). In our initial simulations, T spot is a free parameter between 2500 and 50 K below the T phot for each test star, according to previous spot temperature measurements summarized in Berdyugina (2005) and Strassmeier (2009). We assume three different spot morphologies: (a) a single giant spot, (b) multiple spots, with 5-10 medium size spots randomly distributed on the stellar surface, and (c) small spots, with 100-200 spots that cover the stellar surface. The last case is comparable to solar spots, which are dominated by small-size spots ( S 10 6 1 2 - ) and follow a log-normal distribution (Bogdan et al. 1988;Solanki et al. 2006) as is the mean value of the spot sizes, and σ S is the standard deviation of spot sizes, with σ S =10 ppm from estimates of magnetic active stars (Solanki 1999;Solanki & Unruh 2004;Barnes et al. 2011).
In our simulations, the shapes of each individual spot are circular for the "single giant spot" and "small spots" configurations before placing them onto the stellar surface. For the "multiple spots" configuration, each spot is an ellipse with eccentricity e<0.5. The axes of the ellipse are along the longitude and latitude directions. The upper limit of eccentricity is set as 0.5 to prevent highly elongated spots. The extremely elongated spot along a large range of longitude is rarely detected on the Sun and is not considered in our simulations. The key parameters that modulate the light curves are the spot coverage and their general distribution on the stellar surface. The shapes of individual spots do not affect the large structures of the light curves. Since some spots overlap on the stellar surface, the total spot coverage is slightly smaller than the given f spot when counting these overlapped regions. Figure 2 shows four maps, one of a pure photosphere and three examples of spot morphologies. Morphologies that are azimuthally symmetric, including polar spots and rings, are not shown and would not contribute to any variability. Such longitudinally symmetric structures would produce color and brightness offsets from a spot-free photosphere detectable in spectral decomposition approaches and are discussed in Section 3.4.
Finally, the observed stellar emission (F obs ) is calculated following the "double-cosine" rule (Rackham et al. 2018 ) is the ratio of spot-blended stellar emission (F spot ) and the photospheric emission (F phot ) of each pixel. Light curves are generated by recording F obs step by step with stellar rotations. The darkening process in this paper is simulated by this "double-cosine" rule. The rule itself is a geometric effect when converting the flux from the (θ, f) map. This simplified limb-darkening does not consider the radiative transfer within the stellar atmosphere and hence does not change between the photometric temperatures. Since the core argument of this paper focuses on overall color and magnitude changes introduced by cool spots, rather than on the fine Figure 2. Four two-dimensional (2D) maps of the Monte Carlo simulations, with parameters T phot =4500 K, T spot =4200 K, and f spot =10%. Upper left: a photospheric map with an arrow showing the rotation axis. Upper right: A map for the "single spot" case. Bottom left: A map for "multiple spots," with an example of eight spots on stellar surface, of which three are located at the front side of the star. Bottom right: An example of the "small spots" morphology. The shapes of spots in the "single spot" and "small spots" morphologies are set as circular. Spots in the "multiple spots" case are ellipses with e<0.5. structure of the light curves, we choose to keep this "double-cosine" rule as an approximation of the limbdarkening effect.
In the following analysis, two values are defined to describe the variations of light curves in linear space. ΔF LC is the 10th to 90th percentile amplitude of the simulated light curve, while ΔF spot is the amplitude of the minimum observed stellar flux compared to the flux emerging from an entirely spot-free stellar photosphere with otherwise identical properties. The term ΔF LC is defined similarly as the photometric amplitude in Figure 1, while ΔF spot resembles ΔF spec defined in Equation (1).

Light curve Analysis
The simulated light curves have amplitudes and power spectra that depend on how spots are distributed on the star. The light curves simulated with the simplest spot morphology, the "single spot" case, are shown in three photometric bands in Figure 3. At a fixed spot temperature, the spot filling factor ( f spot ) determines the photometric amplitudes. Meanwhile, with similar f spot and T spot /T phot , spots create larger luminosity variation on cooler stars than warmer stars. Among the light curves, the V-band always shows larger variation amplitude as an indicator of color dependency (see Section 4.2 for more information).
Although the "multiple spots" morphology is geometrically more complicated than the "single spot" case, the methodology of generating light curves is the same. As described in Section 3.2, multiple spots are put onto each simulated stellar surface with a range of spot parameters, and the light curves are generated through synthetic spectra and stellar rotation. Figure 4 presents a few examples of spot configurations with their output light curves. In general, the distribution of spots controls the shapes and the varying amplitudes of the observed light curves. More azimuthal clustering of the spots generates light curves with larger amplitudes, as the top panel of Figure 4. When spots are distributed symmetrically in the longitudinal direction, as in the third panel, the photometric variability is low because one spot becomes visible as another spot rotates to the unseen side of the star. The bottom panel of Figure 4 presents an example of the "small spots" case with 191 randomly distributed spots on the stellar surface. The small asymmetry in the spot distribution leads to a periodic synthetic light curve with a small amplitude, here 2% in the V-band, which is the smallest among all examples.
The frequency analysis of synthetic light curves (the right column of Figure 4) is obtained using the generalized Lomb-Scargle periodogram (Zechmeister & Kürster 2009), which calculates the power over a range of frequencies, assuming that the light curve is well explained by sinusoidal functions. In all cases, the input frequency is the smallest frequency rather than the highest peak. For most simulated light curves, the frequency analysis yields peaks in the power at several different periods because of coincidental resonances. The . Light curves modulated by a single spot with spot filling factors ( f spot =4%, 20%, and 60%) and a fixed spot temperature T spot =3500 K on a K-type star (T phot =4000 K). Three optical bands are adopted in the light curve, shown by different colors (V: black, I: orange, K s : red). As the definition of spot coverage is S spot /S å , the 60% coverage of a single spot means more than half of the star is covered by a spot where F min ∼(T spot /T phot ). 4 identification of multiple periods in light curves reproduces the general behavior of the K2 light curves of the Pleiades analyzed by Rebull et al. (2016b). For K2 light curves with multiple peaks in the power spectra, the light curve amplitudes are generally smaller than stars with single periods, which is well explained with our simulations. Multiple periods may also be detected in unresolved multiple systems.

Comparing Detection Ratio between Morphologies
To quantify the relationship between the spot distribution and the amplitude of the light curve, we define a parameter ζ=ΔF LC /ΔF spot as the ratio between the 10th to 90th percentile amplitude of the light curve divided by the amplitude of the minimum stellar flux against a pure photosphere. The value of ζ quantifies the symmetry in the distribution of spots on the stellar surface by the shape and depth of the light curves, ranging from 0 for symmetric distributions, such as a ring or polar spot, to 1 for highly asymmetric distributions, such as a single spot. This ratio also depends on inclination, since viewing angles that are closer to pole-on reduce variability and lead to a lower value of ζ. Figure 5 compares the amplitude of the light curve to the spot coverage for our simulations. For the "multiple spots" case, the median amplitude ratio ζ is 0.53. Simulations with ΔF spot >50% require clustered spots, leading to large ζ in all such cases. The median ζ-value represents that the variation detected by the light curve is only two-thirds of the realistic brightness changing on the star due to the symmetric distribution of spots. The gray dots in Figure 5 present the "small spots" case with smaller amplitude as well as ζ-values. In other words, if explained by small spots, a 1% amplitude in a K2 light curve corresponds on average to a 10% coverage fraction of starspots.
We add polar spots to these simulations, referred to here as "multiple+polar" simulations, to approximate the large high-latitude spots that are often detected on magnetically active stars (Donati et al. 2008; Barnes et al. 2015). The sizes of the polar spots shown in the "multiple+polar" configuration are fixed to 50% of the total spot coverage on each star. Figure 6 summarizes the ζ-values in log space for three simulated morphologies (multiple, multiple+polar, and small spots) and for the observational results of K2 and LAMOST presented in Figure 1. The distribution of observed log 10 (ζ obs ) is located between the "small spots" and "multiple+polar" case, while the "multiple spot" case results in the largest detection ratio significantly offset from the observed sample. By comparing the "multiple+polar" and "multiple spot" cases, the median detection ratio (ζ) is decreased by 0.22 dex when 50% of the spots are located at the polar region. By enlarging the sizes of polar spots, the median detection ratio becomes smaller, and reaching the observed value (K2/LAMOST) when 70% of the spots are polar spots. Meanwhile, the detection ratio of the "small spots" morphology is an average of 0.1 dex lower than the observed sample. The observed distribution of detection ratios can be explained by a highly symmetric spot morphology, like the "small spots" configuration. The detection ratios for different morphologies are listed in Table 1.
The simulation results confirm our intuition that the bulk of starspots do not induce temporal modulation because of longitudinal resonances. One spot exits the observable stellar hemisphere as another spot enters, balancing the loss of diminished flux. Realistic spot distributions are likely more complicated than our assumptions and likely include large Figure 5. Comparison between the spot covering amplitude, ΔF spot , and the observed variation amplitude in the V-band. The "multiple spots" case is shown by the black dots, while the "small spots" are shown by gray dots. The black and red dashed lines represent the detection rate ζ=ΔF LC /ΔF spot = 1.0 and 0.53, and the blue dashed line is ζ=0.08.

The Observational Effects of Starspots
Spots affect our ability to accurately measure the radius and effective temperature of the star. In this section, we simulate color-magnitude and H-R diagrams for different spot parameters to evaluate how spots will change the measured properties of stars in a cluster. A specific "multiple spots" configuration is adopted in this section since the following discussions are independent of the spot distribution and stellar rotation. We first evaluate the changes in luminosities and colors caused by spots and then create populations to determine the effects on color-magnitude and H-R diagrams.

Luminosity Variations Introduced by Spots
The first consequence of adding cool spots on the stellar surface is the decay of observed stellar luminosity. At a fixed radius, spots decrease the stellar luminosity by where L phot and L spot depend on the respective temperatures. Figure 7 presents the luminosity variation ( versus spot temperature (T spot ) on three example stars, with photospheric temperature (T phot ) of 3500, 4500, and 5000 K. The luminosity variation is sensitive to spot coverage and asymptotically approaches the maximum when T spot /T phot is lower than 0.6. The maximum luminosity variation occurs when no flux is emitted from the spot region, or the so-called non-emitting spot case, resulting in ). In our simulations, the stellar luminosities are calculated by integrating the flux from BT-Settl models over the (θ, f) space, as shown in Equation (4). When considering the fainter photospheric region at the edges (see Figure 4) and limbdarkening effects, the realistic luminosity variation would be slightly larger than the calculation above. However, the basic relationships between L L log D  ( ) and spot parameters are the same.

Color Variation Introduced by Spot Parameters
In this section, we reveal the observational outcomes of various spot parameters, fspot and T spot , on obtained stellar color. To minimize the contribution of the geometric distribution of spots, we only consider the "single spot" configuration here. As an example, the color variation between V and K s bands are displayed in this section, as Δ(V−K s ), which represents the largest color change introduced by certain stellar and spot parameters throughout the stellar rotation. The Δ(V−K s ) is calculated from the synthetic light curves of two colors, and the results are displayed in Figure 8. The synthetic color-color diagram is then compared to the observed colors of stars in the Pleiades for understanding the spot behaviors on cluster members.
The spot coverage dominates the variations of observed colors within certain ranges of T spot . The curves of Δ(V−K s ) versus f spot (the right panels, Figure 8) are similar for different sets of photospheric temperatures. For instance, the maximum Δ(V−K s ) reaches 0.06, 0.15, and 0.28 mag for spot filling factor f spot =0.2, 0.4, and 0.6 when T phot =5500 K.
On the other hand, with a fixed spot filling factor, Δ(V−K s ) always reaches a maximum at T spot that is around 80%-90% of T phot . At higher spot temperatures, the spot and photosphere colors are similar. At lower spot temperatures, the spot flux becomes increasingly negligible and would cause the star to have a fainter luminosity without affecting the colors. The heuristic for T spot scaling matches previous simulations of starspots, in which the choice of T spot is either treated as a free parameter in a certain range (Desort et al. 2007), or as a fixed temperature difference against the photosphere in linear (Barnes et al. 2011) or logarithmic space (Rackham et al. 2018). For realistic measurements of spot temperatures (see Figure 7 and Table5 from Berdyugina 2005), cool stars (T phot <4500 K) always have a spot temperature as T spot ∼0.85T phot that would result in maximum color variations. On the contrary, for warmer stars, only a few detections have such high T spot , while others are around T spot ∼0.7T phot .
Our simulations of the stellar members in a coeval cluster compare the color and luminosity changes in different bands introduced by starspots (Figure 9). Since the peak of the blackbody of the cold spot is located at the longer wavelength, the amplitude changes in bluer bands are larger than redder bands. The relative color changes between B, V, I, and  The effects of existing cool spots on changing the stellar loci on color-color diagrams are presented in Figure 10. The color-changing effects are shown by adding color variations on each sample in the simulation, shown in Figure 9, onto their photospheric colors calculated from the BHAC stellar models . Observed Pleiades members (see Section 5.2 for a description of Pleiades membership) are also shown for comparisons including B-and V-band photometry from Kamai et al. (2014) and K s -band from 2MASS (Skrutskie et al. 2006). Large observed color spreads are seen on cooler stars, indicating strong spots and plages events, which are generated from active stellar magnetic fields.
The empirical colors from Pecaut & Mamajek ( 2013) are also shown in Figure 10 for comparison. A zoomed-in view of low-mass stars (e.g., (V−K s )>4) in the middle panel of Figure 10 shows a notable discrepancy between the empirical main-sequence and photospheric colors. Here, a spot-covered star cluster is generated by adding the simulated color changes from cool spots onto the modeled photospheric colors. The binned simulated samples fit the empirical color of main-sequence stars, suggesting that the discrepancy between empirical colors of low-mass main-sequence stars and photospheric models might come from starspots.
The empirical color difference between young and MS stars found by Pecaut & Mamajek ( 2013) might be related in part to the presence of cool spots. However, to recover this color difference, a 90% coverage of the cool spot is required for the warm photosphere, or (V−K s )<4. For cooler stars, a bluer excess in (B−V ) is needed to fit the lower end of the isochrone. This might relate to faculae or chromospheric activity (Kowalski et al. 2016) that most luminous in B-band. However, the 90% spot coverage is abnormal even on the most magnetic active stars, though an extreme spot coverage is indeed detected on a specific WTTS, LkCa 4 (Gully-Santiago et al. 2017). In addition to stellar activities, gravity might be another cause of the color degeneracy but is beyond the scope of this paper.

Footprints on Color-Magnitude and H-R Diagrams
Starspots cause observational consequences that bias the age determination of pre-main-sequence stars via stellar luminosity. In this section, we discuss the footprints of starspots on the color-magnitude and H-R diagrams, which are sensitive to the spot parameters ( f spot and T spot ). Since the geometric distributions of spots are not important when studying the observed "snapshots" on the color-magnitude and H-R diagrams, only the "small spots" configuration is displayed as other morphologies result in similar luminosity and color spread. All colors applied in the following analysis are based on the empirical colors of main-sequence stars (Pecaut & Mamajek 2013) plus simulated color changes. However, the luminosities are directly calculated from BT-Settl models following Equation (5).
Along with the simulations, six groups of spot parameters are defined to understand how spots with different parameters affect the stellar position for a star cluster on the H-R diagram (listed in Table 2). The parameters of Group 1 and 2 are selected to show the maximum luminosity and temperature spreads introduced by spots, where f spot =60% is close to the maximum detection value in the Pleiades (Fang et al. 2016), while the T spot is chosen as 0.6 T phot and 0.85 T phot according to our previous calculation in Sections 4.1 and 4.2. The median and relatively low spot coverages ( f spot =40% and 20%) are applied in Groups 3-6.
Color-magnitude ((B−V ) versus V ) diagrams are shown in the left and middle panels of Figure 11 with simulated samples and five out of six models (Group 1-5, shown by solid color and dashed lines. 8 ). For warmer stars, the spot-covered stellar positions are all fainter than the isochrones, since starspots affect stellar luminosities more than colors. The modeled curves demonstrate that the variation on V-band magnitude is more sensitive to T spot /T phot (0.6-0.85) than f spot (20%-60%) at the warmest end of the isochrone. However, toward the cooler end, the brightness variation is controlled by f spot .
A color-luminosity diagram is shown in the right panel of Figure 11, in which the stellar luminosities are calculated based on V-band magnitudes and bolometric corrections from Pecaut & Mamajek (2013) based on (B−V ) colors. The absolute ages measured from H-R diagrams are affected by spots. First, evolutionary models that include spots will inflate the radii of stars and raise the observed luminosity, relative to unspotted stars (Jackson & Jeffries 2014;Somers & Pinsonneault 2015). However, counter to this effect, when optical photometry or spectra of spotted stars are interpreted as emission from photospheres with a single temperature, then the measured temperature is hotter than the effective temperature and the luminosity is fainter than the real luminosity. 9 In this way, the measured position of a spotted star should appear older than their age and more massive if the star is fit with multiple components.

The Contribution of Spots to Luminosity Spreads of Young Clusters
The duration of star formation in a single cluster should lead to an age and therefore a luminosity spread on H-R diagrams. While luminosity spreads are measured for all young clusters (e.g., Da Rio et al. 2010;Beccari et al. 2017;Jose et al. 2017), the interpretation of the luminosity spread as an age spread is degenerate with any differences in evolution, stellar variability, and observational errors. The differences in evolution may be caused by a distribution of accretion histories (e.g., Hartmann et al. 1997;Hosokawa et al. 2011;Baraffe et al. 2017), which should be minimal by the age of Pleiades, or by differences in interior structure caused by spots and by magnetic fields inhibiting convections (Somers & Pinsonneault 2015;Feiden 2016). In active star-forming regions, variability in accretion and extinction can also result as luminosity variation or spread (e.g., Venuti et al. 2015;Rodriguez et al. 2016;Contreras Peña et al. 2017;Guo et al. 2017). However, the observational errors of the Pleiades should be minimal, since the Pleiades has no differential extinction, an accurate accounting of binarity (Stauffer et al. 2007), robust membership (Lodieu et al. 2012), and accurate distances from Gaia DR2 astrometry (Gaia Collaboration 2018). The range of spot properties is likely the most significant uncertainty in both the observational measurements of stellar properties in the Pleiades and in the evolutionary models for stars of Pleiades age.
In this section, we first measure the empirical luminosity spread of the Pleiades. We then compare our spot simulations to the empirical luminosity spread to constrain the range in spot properties of low-mass stars in the Pleiades. Finally, we discuss challenges in applying these results to younger regions.

The Empirical Luminosity Spread of the Pleiades
The Pleiades cluster includes over 2000 members (Stauffer et al. 2007;Lodieu et al. 2012;Bouy et al. 2015). In this work, we use the B, V, and I c -band photometry obtained by Kamai et al. (2014). Of the 383 observed members, 61 are identified as binary or multiple systems from their location on the colormagnitude diagram and are ignored here.
We revise the membership classification of the stellar sample in Kamai et al. (2014) using the recent Gaia DR2 parallaxes and proper motions (Gaia Collaboration et al. 2016Luri et al. 2018). By cross-matching the observed samples in Kamai et al. (2014) and Gaia DR2 database by on-sky stellar position, 311 single stars are detected with Gaia parallaxes. The crossmatched targets have parallaxes consistent with a Gaussian distribution at 7.3±0.2 mas, or 137±4 pc in distance, in which the errors are 1σ of the real depths including measurements uncertainties.
To avoid background and foreground contamination on the luminosity spread, we exclude the stars located outside 12 pc (3σ) from the median distance of the cluster of 137 pc. An additional proper motion selection excludes two stars with proper motion 10 mas yr −1 (5σ) away from the median values of 20 and −45 mas yr −1 . In the end, 234 out of 311 single stars from Kamai et al. (2014) are identified as possible cluster  Baraffe et al. 2015). The left panel shows a comparison between simulated and observed colors, while the middle panel provides a zoom-in value of low-mass stars in the simulations. The orange dashed line in the middle panel is the binned simulated stellar color with error bars as the bin size (horizontal) and standard deviation within each bin (vertical). The right panel includes a calculated isochrone with 40% spot coverage on hot stars ((V−K s )<1.7) and 90% spot coverage on others, represented by the red solid line.

Table 2
Spot Parameters members. The contamination rate from the non-members of the Kamai et al. (2014) members is at least ∼25%.
The selected 234 single Pleiades members are located in a color range of 0.4<(V−I)<3.3, or earlier than M5V type. 10 The absolute magnitude is calculated from the individual Gaia DR2 parallax distances for each star. The luminosity spread introduced by treating all members to a uniform 137 pc distance is 0.06 dex, according to the intrinsic spread of σ=4 pc, while the luminosity spread based on the parallax uncertainties for each star is 0.02 dex. The stellar luminosity is then calculated from the V-band absolute magnitude with a bolometric correction from the (V−I) color. Then, an empirical isochrone is fit on the L L log - ) space. The luminosity spread is calculated as the distance on the H-R diagram between each observed sample to the empirical isochrone on the luminosity axis. The distribution of observed luminosities, shown in the middle panel of Figure 12, is well fitted with a Gaussian profile with σ=0.05 in the L L log  ( ) space.

Constraining the Range of Spot Properties of Pleiades Stars and Younger Clusters
To find out a possible origination of the observed luminosity spread, we compare two simulated luminosity spreads to observed samples in Figure 12. The first simulated cluster contains 2000 samples with spot coverages on each star evenly distributed between 1% and 40%, as described in Section 3.1. The distribution of simulated luminosities (the left panel of Figure 12) has a standard deviation, 0.1 dex, two times larger than the observed spread. The other simulated cluster marked as reconstruction in Figure 12 contains 234 samples with spots covering an average of 40% of the stellar surface in a Gaussian distribution with a standard deviation of 10%. The spot   (Kamai et al. 2014) calculated based on V-band magnitude, individual distances from Gaia parallax, and bolometric correction based on (V−I) color. Right: A special simulation to reconstruct the observed luminosity spread. The corresponding spot coverage is around 40% under a normal distribution and the spot temperature is fixed to 60% of the photospheric temperature. 10 A large spread of stellar positions on the V versus B−V diagrams is seen at the lower mass end (Kamai et al. 2014), which might be generated by chromospheric activities and the weak correlation between the color and stellar brightness at the very low-mass end. A smaller spread is also found on (V−I). To avoid this bias, we set a lower mass limit at (V−I)=2.4 (M3V), when calculating the luminosity spread.
temperatures are fixed to 0.6T phot . The "reconstruction" cluster yields a distribution of luminosities that is similar to the Pleiades. From this test, we learned that the luminosity distribution is controlled by the standard deviation in spot coverage, while the average spot coverage is not constrained. The low scatter in the stellar luminosity suggests that for the Pleiades, the spot coverage is similar for all members of a given stellar mass.
The standard deviation of 0.05 dex in luminosity at any single age is likely the minimum luminosity spread for young clusters. The contribution of spots at younger ages may be larger. The variation amplitudes of the light curve of some weak-lined T Tauri stars can reach 0.3-0.5 mag, while others have little variation (e.g., Herbst et al. 1994;Grankin et al. 2008). These large variations may be explained if younger stars have strong dipole magnetic fields (Gregory et al. 2012) and therefore larger dispersions in estimated luminosities. Such large spots could help to explain the 0.2-0.3 dex luminosity spread seen in all young clusters (e.g., Da Rio et al. 2010;Preibisch 2012;Herczeg & Hillenbrand 2015). By the age of the Pleiades, the magnetic structures for all stars of a given mass may have converged to similar spot coverages, leading to a smaller spread in luminosities.

Summary
A large observational discrepancy is seen on spot coverage between different methodologies. The spot coverage obtained from light curve amplitude is strongly underestimated. In this paper, we presented a series of Monte Carlo simulations for spot distributions on 125 Myr old young cluster members with masses ranging between 0.08 and 1.32 M e , including morphologies such as single, multiple, circumpolar, and small spots.
In a simple one-star simulation, we learned that color and luminosity changes are dominated by spot filling factor. For multi-band light curves, the B and V bands are more sensitive to spots than the T and K s bands and the spot filling-factor is crucial to determine the variation amplitude. The color variations reach the maximum when T spot /T phot ∼0.85.
For "multiple spots" and "small spots" configurations, a detection ratio (ζ=ΔF LC /ΔF spot ) is defined to quantify how much spot coverage is seen from the light curve. In the "multiple spots" morphology, half of the spot variations are measured on light curves. When the spot is distributed symmetrically along the longitude, the "small spots" case, 90% of the spot behavior is hidden. The detection ratio between K2 light curve is 18% and LAMOST spectra, close to "small spots," suggesting that most of the spot coverages are not detected.
The spot-modified stellar loci are shown on color-color, color-magnitude, and color-luminosity (H-R) diagrams. In general, the spot-covered stars are redder and fainter than the photosphere. The luminosity spread of the observed samples is well reconstructed by cool spots on star surfaces with 10% standard deviation in spot coverage. Large spots introduce luminosity spreads up to 0.1 dex on young clusters, while the measurements on spot coverage through K2 light curves are biased by their symmetrical distribution.
Several effects, including faculae associated with dark spots, the inflation of stellar radii as a result of the energy conservation, the latitudinal uneven distribution of starspots, and the lifetime of spots, as well as other complicated configurations beyond and between our default settings, are not included in this work. However, our Monte Carlo simulations provide a view of spot distribution, temperature, coverage, and the corresponding luminosity spreads that well match the observational results.