Correlating Changes in Spot Filling Factors with Stellar Rotation: The Case of LkCa 4

We present a multi-epoch spectroscopic study of LkCa 4, a heavily spotted non-accreting T Tauri star. Using SpeX at NASA's Infrared Telescope Facility (IRTF), 12 spectra were collected over five consecutive nights, spanning $\approx$ 1.5 stellar rotations. Using the IRTF SpeX Spectral Library, we constructed empirical composite models of spotted stars by combining a warmer (photosphere) standard star spectrum with a cooler (spot) standard weighted by the spot filling factor, $f_{spot}$. The best-fit models spanned two photospheric component temperatures, $T_{phot}$ = 4100 K (K7V) and 4400 K (K5V), and one spot component temperature, $T_{spot}$ = 3060 K (M5V) with an $A_V$ of 0.3. We find values of $f_{spot}$ to vary between 0.77 and 0.94 with an average uncertainty of $\sim$0.04. The variability of $f_{spot}$ is periodic and correlates with its 3.374 day rotational period. Using a mean value for $f^{mean}_{spot}$ to represent the total spot coverage, we calculated spot corrected values for $T_{eff}$ and $L_\star$. Placing these values alongside evolutionary models developed for heavily spotted young stars, we infer mass and age ranges of 0.45-0.6 $M_\odot$ and 0.50-1.25 Myr, respectively. These inferred values represent a twofold increase in the mass and a twofold decrease in the age as compared to standard evolutionary models. Such a result highlights the need for constraining the contributions of cool and warm regions of young stellar atmospheres when estimating $T_{eff}$ and $L_\star$ to infer masses and ages as well as the necessity for models to account for the effects of these regions on the early evolution of low-mass stars.


INTRODUCTION
Accurate age determinations for pre-main-sequence (PMS) stars are essential to understanding the formation and evolution of stars and planetary systems. One popular method for constraining the ages of PMS stars relies on a direct comparison between the ob-served stellar luminosities and surface temperatures and those predicted by theoretical evolutionary models (e.g., D'Antona & Mazzitelli 1994;Soderblom et al. 2014; Baraffe et al. 2015). Unfortunately, young stars are complex systems often characterized by strong magnetic fields, rapid rotation rates, excess emission from circumstellar material, mass accretion onto the stellar surfaces, and mass outflow from disk and stellar winds (e.g., Hartmann et al. 2016). Chromospheric and coronal activity are heightened, leading to strong flares producing large fluxes of high-energy photons (e.g., Feigelson & Montmerle 1999;Petrov et al. 2011). As part of the heightened activity, the large-scale inhibition of convec-tion by strong magnetic fields in these systems is quite possible, resulting in the formation of starspots covering significant fractions of the stellar surfaces. Given the episodic or transient nature of these phenomena, young stellar systems exhibit variability on timescales as short as hours across all wavelengths. Such activity often limits our ability to constrain otherwise straightforward observable stellar parameters. Ages inferred from comparisons of effective temperatures and stellar luminosities to those predicted by standard evolutionary models (e.g., Baraffe et al. 2015) typically result in large spreads in the ages of stars residing in the same cluster. While some of this spread may be due to different star formation epochs that have occurred within the same region, ignoring the effects of spots on the observable quantities and on stellar evolution likely contributes to this spread, confusing our understanding of the star forming history. Therefore, the presence of large starspots on the surfaces of a sizable fraction of young stars and the lack of evolutionary models that account for spots are potentially responsible for some of the spread in the ages, masses, and evolutionary statuses inferred for stars in a given young cluster (Preibisch 2012;Soderblom et al. 2014).
Large complexes of cool spots rotating with the surfaces of low-mass PMS stars produce periodic variability, which can be used profitably to measure rotation periods for these objects (e.g., Bouvier et al. 1995;Herbst et al. 2007;Grankin et al. 2008). Such spots also provide a reasonable explanation for the systematic color anomalies and optical/infrared spectral type mismatches observed for T Tauri stars, both accreting and non-accreting (Gullbring et al. 1998;Vacca & Sandell 2011;Debes et al. 2013; Bary & Petersen 2014;Czekala et al. 2015;Kastner et al. 2015;Gully-Santiago et al. 2017). Debes et al. (2013) and Bary & Petersen (2014) demonstrate that the near-infrared (NIR) spectra of TW Hya and DQ Tau, respectively, can be modeled with empirical composite spectra made from a weighted average of two standard star spectra -a warmer standard representing the photosphere and a cooler one representing the spot. In both cases, the authors find that cool spots may cover over 50% of the surfaces of the stars. Donati et al. (2014) use multi-epoch spectropolarimetric observations (R∼65,000) of LkCa 4, a weak-line T Tauri star and the subject of the study presented herein, to construct tomographic maps of the stellar surface and to study the magnetic topology of the star. Their results indicate the presence of large cool spots covering more than 20% of the surface as well as the existence of large warm plages. Herczeg & Hillenbrand (2014) use TiO features in the near-IR to revise the spectral type of LkCa 4 from a K7 to a later M1.5, likely highlighting the effect of spots on single-band temperature measurements. Gully-Santiago et al. (2017) (from here on GS17) also observe LkCa 4 at high-spectral resolution (R∼45,000) in the near-IR with IGRINS. Applying a two-temperature atmospheric model to fit their data as well as the TiO bands in the spectra of Donati et al. (2014), GS17 establish the presence of a large spot or spot complex that covers nearly 80% of the stellar surface. The surprisingly large discrepancy between filling factors determined by Donati et al. (2014) and GS17 can be reconciled by the fact that the Zeeman Doppler Imaging (ZDI) technique employed by Donati et al. is insensitive to collections of smaller spots.
Using high-resolution iSHELL spectra (R∼47,000), Flores et al. (2019); Flores et al. (2022) measure magnetic field strengths on the surfaces of T Tauri stars and correlate magnetic field strengths to measurements of stellar temperatures. The results of Flores et al. indicate that a correlation likely exists between spots and effective temperatures of these sources.
Collectively, these studies highlight the uncertainty and complexity that starspots introduce when using evolutionary models to infer ages and masses of highly active PMS stars. They also illustrate the importance of developing new evolutionary models of spotted stars, incorporating the physical mechanisms that produce the spots as well as predicting their impact on the evolution of young stars (e.g., Feiden 2016;Somers et al. 2020). Such models will improve our efforts to confidently and accurately infer the ages of PMS stars and the clusters in which they form. Constraining these models will require a simple and direct method for determining spot filling factors and spot temperatures for large samples of PMS stars across the mass spectrum.
Toward this goal, we present a multi-epoch, mediumresolution, NIR spectroscopic study of LkCa 4 in which we constrain spot sizes and temperatures and correlate the changes in the spot filling factors with the rotational phase of the star. Using our best-fit model parameters for photospheric temperature (T phot ), spot temperature (T spot ), and spot filling factors (f spot ), we are able to reproduce the V -band variability observed during a time frame that overlaps with our spectral observations. We show how our results compare with the studies mentioned above, which were conducted at much higher spectral resolution (Donati et al. 2014;Gully-Santiago et al. 2017). The observations presented benefit from consistent temporal coverage over five consecutive nights or roughly 1.5 stellar rotations. Although the absolute value of the total spot coverage depends on the modeldependent photosphere and spot temperatures, the tem-poral coverage permits us to better constrain the total fraction of the stellar surface covered by spots than that of a single observation.
First, we outline the observations and calibration steps in Section 2. In Sections 3 & 3.1, we describe the empirical composite spectral models and determine the bestfit filling factors and temperature ranges for the photosphere and the spots. In Section 3.2, we present the strong correlation we find between the photometric variability and the changes observed in the spot filling factors suggesting our observations and spectral models are sensitive to the rotation of the star. In Section 4, we take a small digression to explore another spectral type indicator, the D CO spectral index based on the 2.29 µm CO bandhead (Mármol-Queraltó et al. 2008) to test its sensitivity to spots and its consistency with other temperature indicators. In Section 5, we use the best-fit model parameters to calculate spot-corrected T ef f and L . We then place these spot-corrected values on the HR diagram alongside evolutionary tracks and isochrones predicted by both standard and spotted star evolutionary models (Baraffe et al. 2015;Somers et al. 2020) and discuss the results.

OBSERVATIONS
We observed LkCa 4 using SpeX, a medium-resolution cross-dispersed NIR spectrograph at NASA's Infrared Telescope Facility (IRTF) atop Maunakea over five consecutive nights on UT January 6-10 2019. Using the short-wavelength cross-dispersed mode (SXD; Rayner et al. 2003) with the 0. 3 × 15 slit (R∼2000), we collected a total of 12 spectra of the target as part of a larger program to study starspots and accretion activity in PMS systems. The SXD setting provides continuous wavelength coverage from 0.7-2.55 µm. The goals of this project were achieved due to favorable weather conditions, which permitted consistent monitoring of more than one full rotation of LkCa 4.
The data were collected using an AB nod sequence typical of long-slit near-IR spectra acquisition. The subtraction of the 2D spectral image pairs allows for the efficient removal of terrestrial OH emission lines, background, and dark current. A0V telluric standards HD 27761 and HD 24000 were observed close in airmass (∆ sec z ≤ 0.1) to the target and were used to remove telluric absorption features and to calibrate the target spectra. Flat-field corrections, wavelength calibrations, spectral extraction, co-adding, telluric corrections, and merging of the spectral orders were performed using SpexTool v4.0.5, an IDL-based reduction package described by Cushing et al. (2004). Details of the observations can be found in Table 1. Sample spectra cov-  Figure 1, highlighting the nature of the variability observed in the shorter wavelength regions of the spectra.

EMPIRICAL MODELS OF SPOTTED STARS
We constructed two-temperature models of spotted stars as empirical composite spectra following the same procedure outlined in Debes et al. (2013) and Bary & Petersen (2014). We use the term empirical composite to clearly indicate that these model spectra are not generated from synthetic stellar atmospheric models. Instead, they are produced using the spectra of standards found in the SpeX IRTF Library (Rayner et al. 2009;Cushing et al. 2005). In the two-temperature models, the spectrum of the warmer standard represents the photosphere, F λ (T phot ), and the cooler standard represents the spots, F λ (T spot ). Therefore, changing model parameters T phot or T spot is achieved by selecting different spectral standards from the SpeX Library. The composite spectral models are generated using the following (1) where F λ (T phot ) and F λ (T spot ) are normalized template spectra, f spot is the instantaneous spot filling factor (i.e., the fraction of the observable stellar surface that is covered by spots following the convention adopted by GS17), and C bb is a scaling constant defined as the ratio of the Planck functions of the photosphere and the spot at 1.1 µm. C bb approximates the relative normalized flux units of the two spectra based on their effective temperatures. Before combining, the wavelength arrays of the standard star spectra are aligned through a onedimensional interpolation using the interp1d algorithm found in SciPy (Jones et al. 2001). The resulting composite model spectrum is renormalized before fitting to the target spectra.
We have chosen to use dwarf spectral standards when constructing the models similar to Debes et al. (2013) and Bary & Petersen (2014). The spectral standards used to construct the spotted star models as well as the spectral types, effective temperatures, variability status, and B − V color excesses for each source are listed in Table 2. We note that M dwarfs are well known for their aperiodic variability due to strong stochastic flares (Hartman et al. 2011). We find that six of the M dwarf standards selected to represent the spots are identified as eruptive variables. However, we proceed with using the M dwarfs as standards assuming that the shortlived nature of the eruptions is not likely to impact the single-epoch observations presented in the SpeX library. Rayner et al. (2009) measure color excesses for the library stars and do not present dereddened spectra for stars with E(B − V ) < 0.108. The color excesses quoted for Gl 406 (M6V) and Gl 466C (M7V) are significant. Therefore, these standard star spectra were dereddened with a standard interstellar extinction law (Fitzpatrick 1999) prior to using them to construct spectral models.
In general, optical and infrared TiO and FeH features are temperature sensitive and are considered good, yet complicated spot indicators (Herbst & Levreault 1990;Neff et al. 1995;O'Neal et al. 1996;Schiavon et al. 1997). Therefore, it is important to note different sensitivities between TiO and FeH that may impact the individual constraints they place on the best-fit composite spectra and values for f spot . For instance, the FeH Wing-Ford band at 0.99 µm is sensitive to changes in surface gravity with the feature appearing to be strongest in the spectra of the coolest dwarf stars (Schiavon et al. 1997). In fact, its presence and strength in NIR spectra of unresolved stellar populations have frequently been used to measure the contribution from cool, dwarf stars (e.g., Couture & Hardy 1993;Schiavon et al. 1997;Cenarro et al. 2003). However, with regards to its sensitivity to small changes in log g, Bary & Petersen (2014) found little difference in the strengths of FeH between synthetic spectra of a K5IV with log g = 3.5, representing a "puffy" T Tauri star, and a K5V dwarf star with log g = 4.5 (Coelho et al. 2005). Such a similarity likely indicates that larger differences in the surface gravity are required to produce an effect on the FeH band strengths to be detectable with moderate resolution spectroscopy.
In addition to surface gravity effects, the magnetic sensitivities of both TiO and FeH molecular states and transitions to Zeeman effects is also an important consideration. Absorption features associated with both bands have been shown to be quite sensitive to magnetic fields and have been used to measure magnetic field strengths on K-and M-type stars (e.g., Afram & Berdyugina 2015). FeH has gained considerable attention as a probe of magnetic field strengths on M dwarfs, which are too cool to possess strong, magnetically sensitive atomic features (Valenti & Johns-Krull 2001;Reiners & Basri 2006;Shulyak et al. 2014;Afram & Berdyugina 2019;Kochukhov 2021).
Given that the composite spectra are constructed with template spectra of M dwarfs representing the cooler spotted regions that likely possess magnetic fields that are stronger than the non-spot regions, it is important to acknowledge that the M dwarf templates possess TiO and FeH features that are affected by strong magnetic fields. For instance, Afram & Berdyugina (2019) find a range of 3-6 kG fields with an average of 5 kG for a a All temperatures are from Pecaut & Mamajek (2013).
b Color excesses reported in Rayner et al. (2009). Previous studies of spotted T Tauri stars have indicated that the spectral types representing T phot in twotemperature models will be similar to the optically derived spectral types for the stars (Debes et al. 2013;Herczeg & Hillenbrand 2014;Bary & Petersen 2014;Gully-Santiago et al. 2017). Therefore, the spectral templates chosen to represent F λ (T phot ) in our composite models bracket the K7V spectral type reported for LkCa 4 in the literature (Herbig et al. 1986;Strom & Strom 1994;Hartigan et al. 1995;White & Ghez 2001;Grankin 2013). We constrained the parameter F λ (T phot ) by selecting four photospheric templates: M0V p , K7V p , K5V p , and K4V p (3850 K ≤ T phot ≤ 4600 K). The spectral types of the templates representing F λ (T spot ) were confined to the range between M8V s and M1V s (2570 K ≤ T spot ≤ 3660 K). This range of spot temperatures encompasses the values suggested by previous studies (i.e., GS17) and fits the typical spot-tophotosphere temperature ratios (Strassmeier 2009;Fang et al. 2018). The contribution of the spot template to the composite spectrum is weighted by the spot filling factor, f spot , which we allow to vary from 0.0 to 1.0. It is important to note that cooler spots with smaller filling factors will mimic warmer spots with larger filling factors leading to an inherent degeneracy in these models.

Model Fitting
We searched for the best-fit model for each of the 12 epochs of LkCa 4 spectra by varying four model parameters: T phot , T spot , f spot , and A V . The values for T phot , T spot and f spot were constrained by the parameter space defined above in Section 3. We let A V vary between 0.0 and 1.0 in steps of 0.1 to encompass the two values of 0.35 (Gully-Santiago et al. 2017) and 0.69 (Kenyon & Hartmann 1995) reported in the literature. The LkCa 4 spectra were dereddened with the same standard interstellar extinction law used for the M dwarf standards.
We selected three spectral windows to constrain the best-fit models. The first is a broad spectral window stretching from 0.8 1 to 1.35 µm, which we will designate F 0.8−1.35 µm . The large wavelength coverage of this window will force the best-fit models to accurately reproduce the shape of the continuum in a region of the spectrum that is most significantly affected by interstellar reddening. The other two spectral windows center on two spot-sensitive molecular absorption features: TiO band (λ=0.845 0.870 µm) and the Wing-Ford FeH band (λ=0.9851.02 µm). We note that both of the molecular features are included within the larger spectral window. Given the comparatively large wavelength coverage of the F 0.8−1.35 µm removing either or both of the features from the window while performing the fitting routine (see below) affects the χ 2 red values 2 on the order of 0.01%. Therefore, we simply perform the fits to this window without excluding the TiO and FeH absorption features. We will refer to these three spectral regions as spot indicators or indicators as shorthand.
An initial round of fits between the composite models and the target data were performed using the Levenberg-Marquardt minimization algorithm found in lmfit (Newville et al. 2014). The χ 2 red values for these fits are systematically large due to the noise in the composite models and the consistently poor fits to the atomic absorption features. The χ 2 red values did not exactly follow a normal distribution. The two-temperature models with f spot that possessed the lowest χ 2 red values were then used as the starting point for the walkers in an emcee Markov Chain Monte Carlo sampler (Foreman-Mackey et al. 2013). The MCMC procedure probed the posterior probability density function of f spot . In practice, the posterior probability density functions resembled normal distributions with maxima that closely corresponded to the models with the lowest χ 2 red values. We have adopted the 1σ widths of these distributions as the uncertainties in f spot . This was performed for every observation and all parameter space, yielding the uncertainties in f spot .
For the three different spot indicators, we find good agreement between the best-fit model parameters, T phot , T spot , and f spot constrained by the F 0.8−1.35 µm spectral window and TiO feature.
In Figure 2, we graphically illustrate the goodness of fit by presenting comparisons of one LkCa 4 spectrum to three composite models with different f spot values in the 0.8-1.35 µm window. By visual inspection, the K5V p +M5V s model with f spot = 0.88 and χ 2 red ≈ 2.3 is an overall better fit than the K5V p +M4V s and K5V p +M6V s models. For a difference in χ 2 red of ∼1.1 and 0.44, significant differences in the residuals and the quality of the fits to the strengths of the TiO features and the general shapes of the 0.8-1.35 µm region of the spectra are evident. Similar to other two-temperature models presented in Bary & Petersen (2014), atomic absorption features are not well fit and contribute to the large χ 2 red as compared to the fits to the TiO or FeH features. Similar differences between the goodness of fit for the K7V p +M4V s , K7V p +M5V s , and K7V p +M6V s models are observed with the K7V p +M5V s producing the lowest χ 2 value.
It is apparent in Figure 2, that the models do not fit the FeH feature in a predictable manner or one that is consistent with the other spot indicators. On average over the 12 epochs of observations, the FeH fits point to three possible best-fit models with considerably different values for f spot . For the K5V p +M3V s model, the value of f spot is 1.0 essentially selecting an M3V spectrum as the best fit to all 12 epochs. By contrast, the K5V p +M5V s and K5V p +M6V s models were equally good fits with f spot values falling in the range of 0.72-0.85 over the 12 epochs. Similar behavior of the FeH fits was observed for the models in which the K5V was replaced with a K7V standard. The spot filling factors decreased as expected for a model with a cooler photosphere. To illustrate the goodness of fit to the FeH feature, we present a similar plot for the FeH comparing four model fits in Figure 3. The K5V p +M4V s fits are the poorest due to the mismatch between the model and the data between 0.985 and 0.990 µm. The other three models do a far better job of fitting these short wavelengths and differ mostly in the way they fits the small components of the feature. The similarities between the other three fits and the outlier nature of the K5V p +M4V s points to a potential problem with the M4V standard. If this is the case, then the FeH feature may fit all four of these models equally well over a range of filling factors rendering it a less useful spot indicator than F 0.8−1.35 µm and TiO. Therefore, we will treat the best-fit models based on these spot indicators as the most reliable, but will include model parameters associated with the FeH indicator in the following discussion where we believe it is useful and instructive.
In Table 3, we present the values for f spot values associated with the best-fit models constrained by the TiO and F 0.8−1.35 µm . The values listed for each of the 12 observations correspond to two photospheric templates, F λ (K5V p ) and F λ (K7V p ), combined with one spot template, F λ (M5V s ). The f spot values listed for FeH correspond to the same F (T phot ) and F (T spot ) and do not represent the model with the minimum χ 2 values obtained using this indicator. In the Appendix, Tables A1 and A2 we provide a more complete listing of the bestfit f spot values over a larger range of spot temperatures for all three spot indicators.
Why do the spot indicators potentially point toward significantly different f spot values for a given pair of T phot and T spot ? The most straightforward answer seems to be that the moderate resolution spectra of the FeH feature does not distinguish between the twotemperature models as well as the TiO feature and the large spectral window. We also suspect that the differences in the dissociation energies of TiO and FeH (D TiO = 7.26 eV; D FeH = 2.9 eV; Wahlbeck & Gilles 1967;Wang & Angelici 1996, respectively) may contribute to this discrepancy. For instance, FeH will not form within the photosphere of a K-type star and will only exist within the cooler spot regions. In addition, there may be some differences in magnetic sensitivity as well as the strength of the magnetic fields impacting the line strengths of these molecules. Again with the FeH feature being produced predominantly in the spotted regions and the TiO features formed both within and outside of the spot, we speculate that the relative line strengths of these features to be more complex than can  be described by our relatively simple two-temperature models that lack any specific magnetic field parameters. Finally, as described above, the unknown characteristics of the standard stars given their variability and surface magnetic activity may also affect the composite models.

Correlating f spot Variability with Stellar Rotation
In Figure 4, we present nearly seven years of AAVSO V -band photometry collected between UT 2013 December 24 and UT 2020 October 17 along with the bestfit values for f spot for TiO and F 0.8−1.35 µm as well as the corresponding values for FeH. All values have been phase-folded setting the time of our first observation taken at JD 2458488.71352 as φ = 0.0 and using P rot = 3.374 days (Grankin et al. 2008). The periodic variations in the spot filling factors are positively correlated with the periodic variability in the V -band light curve just as one would expect if the V -band variability were due to spots rotating with the surface of the star. Despite the differences in the absolute values of f spot derived from the TiO, F 0.8−1.35 µm , and FeH for models with similar values of T phot and T spot , the variations within those values correlate with rotational phase (Figure 4).
Given what is known about spot lifetimes on young, active stars, we assume that the spot complex(es) on LkCa 4 are stable during the five nights over which these data were collected. Therefore, we conclude that the observed variability in the instantaneous f spot values is due to the rotation of the star and not due to periodic changes in the total spot coverage or the spot temperatures. In addition, we also rule out contributions to the variability from circumstellar material including phenomena such as inner disk warps and/or accretion flares (Covey et al. 2021) based on the lack of evidence for a circumstellar disk and accretion activity in the LkCa 4 system (e.g., Andrews & Williams 2005).

Predicting ∆V from Best-fit Model Parameters
Next we test whether or not our best-fit spotted star models can reproduce the V -band variability observed in the AAVSO light curve from UT 2019 given that the time frame aligns with the SpeX observations. We calculated the minimum to maximum V -band variability, ∆V , using the following equation where F λ (T p ) and F λ (T s ) are the flux densities of the photosphere and spot components, respectively, f max and f min are the maximum and minimum spot filling factors, and S V (λ) is the transmissivity of the Johnson-Cousins V -band filter as defined in the General Catalog of Photometric Data and revised by Mann & von Braun (2015). We used BT-Settl(CIFIST) (Allard 2014) synthetic spectra over the relevant wavelengths to represent the F λ (T p ) and F λ (T s ) contributions to the V -band magnitude. We calculated four values of ∆V calc , one for each of the f spot values derived from the F 0.8−1.35µm and TiO indicators for models: K7V p +M5V s and K5V p +M5V s . Values for the mean spot filling factor, f mean spot , and its amplitude, f amp spot , were extracted from a sinusoidal fit to the phase-folded f spot curves in Figure 4. These values are listed in Table 4.
The observed value for the min-to-max variability in the 2019 AAVSO light curve, ∆V obs = 0.526 ± 0.032, agrees to within the uncertainties of ∆V calc = 0.56 +0.10 −0.07 and ∆V calc = 0.46 +0.09 −0.06 calculated for the f spot values associated with the TiO indicator and the best-fit models, K5V p +M5V s and K7V p +M5V s , respectively. The  The agreement between the amplitudes of variations in the spot sizes and the value of ∆V for the TiO indicator provides a reasonable consistency check for the twotemperature model parameters f spot , T spot , and T phot .

Year-to-year Variations of Spot Coverage
The 2013-2019 AAVSO light curve presented in Figure 4(d) possesses a few interesting aspects that suggest the modulation and possible evolution of the spot complexes during this time period. The phase-folded light curve possesses a vertical width of δV = 0.1-0.2 mag at all phases, though it appears to be largest near φ = 0.5. Such a modulation in δV indicates that the average spot coverage of the stellar surface has changed over this time frame. The slight asymmetry of the variations in the ver-tical width, particularly between φ =0.2 and 0.4 suggests an offset or shift in rotational phase, likely indicating a slight alteration in the timing of the minima and maxima over year-long timescales. We estimate the shift in rotational phase to vary between δφ ∼ 0.008 and δφ ∼ 0.084 over the 7yr period. Given the small uncertainties in the photometric measurements, these variations in the phased light curve are likely real features. We do not interpret these phase shifts necessarily as a change in the accepted 3.374 day rotation period of the star. Instead, we believe that they are most likely the result of the spots or spot complexes forming and dissipating at or migrating to different latitudes and/or longitudes. Similar and even more substantive secular changes to the spot coverage were highlighted in three decades of LkCa 4 observations (UT 1986(UT -2016 presented in GS17.
Working from the photometry and assuming that the T phot and T spot remain constant over time, we can use V mean and ∆V obs in the AAVSO data to estimate f mean spot and f amp spot in previous years. Assuming that the spot coverage is relatively constant on a year-long time scale, the AAVSO photometry from UT 2013-2019 was divided in 1yr periods. In Figure 5, we present each of the 7 yr long light curves with a best-fit sinusoidal function determined using the same minimization and MCMC algorithm adopted for the spectral fitting. For each light curve, V mean has been subtracted to permit a direct comparison of the year-to-year variations of the amplitude.
In addition, small shifts in rotational phase measured relative to the minima and maxima of the UT 2019 curve also have been removed. The values for V mean , ∆V obs , and ∆φ are listed in Table 5.  In order to determine f mean spot and f amp spot from previous years of V -band data, we associate the value for f mean spot derived from our models with V mean from 2019. Differences between the UT 2019 V mean and previous year V mean are used to determine the values of f mean spot . The ∆V value measured from the light curve for a previous year can then be used to determine f amp spot (see Table 5). Also presented are the ratios of the maximum to minimum photospheric filling factors, f max phot /f min phot . While directly related to the ratio of the values of f spot , we include this value to highlight the strong correlation between ∆V and the changes to f phot . One can understand the nature of this correlation better by considering two stars that have the same T phot and T spot , but possess the different mean filling factors of f mean spot = 0.7 and 0.8. For both stars, assume that the filling factors vary by the same amount, ∆f spot = 0.1 or 10%, such that the ranges of filling factors are 0.65 ≤ f spot ≤0.75 and 0.75 ≤ f spot ≤0.85, respectively. The star with the larger spot coverage would appear dimmer on average with a higher V mean and would exhibit a larger amplitude of variability despite the same change in f spot . The reason being that the photospheric emission dominates the brightness of the star in the V -band of even the most heavily spotted stars. Therefore, the amplitude of variability depends most directly on the ratio of the maximum and minimum filling factors of the photosphere, f max phot /f min phot and less on the absolute change in f spot . Hence, a 10% change in f spot for the star with a larger f mean spot leads to a larger fractional change in the portion of the star covered by the hotter photospheric region. Such a change leads to a larger overall dim- ming/brightening of the star and a larger amplitude of variability.
The secular changes present in the AAVSO data are depicted in Figure 6, which plots the values of (a) V mean , (b) ∆V obs , (c) the extrapolated minimum and maximum values for f spot , and (d) f max phot /f min phot by year. Within uncertainties, we find that the star is brightest when it displays the smallest ∆V . These values for the V mean , the ∆V values, and the filling factor ratios agree with those of GS17, where a roughly twofold increase in the photospheric filling factor is required to cause a ∆V = 0.6.

CO SPECTRAL INDEX TEMPERATURES
The SpeX SXD spectra provide excellent wavelength coverage of the K band, which includes several strong CO bandheads. The strength of the shortest wave- length CO band at 2.29 µm is used as an indicator for T ef f , [Fe/H], and log g and often applied to unresolved spectra of stellar populations in clusters and galaxies. Mármol-Queraltó et al. (2008) updated the definition of a CO spectral index for the 2.29 µm feature, which they designate as D CO . Given the impact of spots on spectral type determinations of heavily spotted PMS stars, we apply the D CO index to our spectra of LkCa 4 for comparison.  Table 6 lists all of the D CO and corresponding T ef f values.
The mean temperature is two subclasses later than a K7V (4100 K), the spectral type reported in Donati et al. (2014). The value more closely matches the temperature of 3670 K determined by Herczeg & Hillenbrand (2014). The D CO temperature, lying between the photospheric and spot temperatures and in relatively good agreement with Herczeg & Hillenbrand (2014) value, suggests that the D CO index is somewhat sensitive to the presence of the cool spots. However, we find in the following analysis of spot-corrected temperatures that the D CO temperature is significantly warmer.

SPOT-CORRECTED T EF F , L , AND SPOTS EVOLUTIONARY MODELS
Next we calculate spot-corrected values for T ef f and L from the values derived for f mean spot , T phot , and T spot . First, we adopt f mean spot values and uncertainties from the best-fit composite models, K7V p +M5V s and K5V p +M5V s , as estimates of the total spot coverage of LkCa 4 with the understanding that a portion of the star is never visible due to the inclination of the stellar rotation axis (i = 35°; GS17). Values for T ef f were calculated using the following for the model parameters associated with spot indicators: F 0.8−1.35 µm and TiO. Adopting a stellar radius of 2.3 R (GS17), we calculated the corrected values for L . The uncertainties on T ef f were estimated by assuming one-half of a subclass uncertainty on the spectral types and associated temperatures for both parameters T spot and T phot . The resulting extrema of T ef f values were then used to calculate the upper and lower bounds for the corrected L values. The values derived for T ef f and L are listed in Table 7. In Figure 8, the corrected T ef f and L values with the uncertainties obtained from the F 0.8−1.35 µm indicator are placed on an H-R diagram alongside literature values taken from Donati et al. (2014), Herczeg & Hillenbrand (2014), GS17, and two values derived using empirical relationships (detailed below) found in Flores et al. (2022). While our two data points with uncertainties reflect only the mean value of f spot with its uncertainties, the gray parallelogram encompasses the entire parameter space associated with f min spot and f max spot measured over the full rotation of the star. Therefore, any single-epoch observation used to correct T ef f and L should place the star within the gray parallelogram.
Overlaid are the isochrones and tracks from the Baraffe et al. (2015) standard evolutionary models (orange), as well as those from the Stellar Parameters of Tracks with Starspots (SPOTS) models for f spot = 0.85 (black, Somers et al. 2020). The SPOTS evolutionary models incorporate the structural effects of large cool starspots by accounting for phenomena such as the inhibition of convection by strong magnetic fields and the impact spots have on the pressure of the stellar photospheres. SPOTS evolutionary tracks are calculated for low-mass stars with spot filling factors as large as 0.85, making LkCa 4 a perfect candidate for comparison to these models.
First, we note the position of the non-corrected optical T ef f and L values from Donati et al. in Figure 8, which places LkCa 4 between the 0.7 and 0.8 M standard evolutionary tracks with an age in the range of ∼1-3 Myr. The cooler value for T ef f from Herczeg & Hillenbrand (2014) shifts the masses down to 0.4 M and the age to ∼ 0.5 Myr. The rest of the values plotted are in some sense corrected for the existence of spots or cooler regions of the star and shift the values of T ef f farther to the right and slightly downward on the H-R diagram.
Relative to the Baraffe+2015 standard evolutionary models, the corrected values yield lower mass ranges of 0.25-0.30 M with ages less than 0.5 Myr. These new values represent a shift in mass by a factor of two or slightly greater and a decrease in age by a factor between two and six.
Comparing the two sets of evolutionary models, the tracks for the SPOTS models are shifted to cooler temperatures and higher luminosities. Given that the corrected T ef f and L moves the star down and to the right, it appears that accounting for the evolutionary effects of spots mitigates some of the shift toward lower masses and younger ages. As such, the corrected placement of LkCa 4 in relation to the f spot = 0.85 SPOTS models increases the mass to a range of 0.45-0.60 M and the age to a range of 0.5-1.25 Myr, relative to the Baraffe+2015 models.
We now turn our attention to the two data points obtained using the empirical relationships derived in Flores et al. (2022). Based on the correlations these authors find between the optical and infrared temperatures as well as the optical temperatures and the strengths of the stellar magnetic fields, they develop two empirical equations for the shifts in temperature: (1) ∆T opt−ir = 0.36T opt −1170 K and (2) ∆T B = 206B − 135 K. Note that the uncertainties in the Flores+2022 relations are large, so we chose only to use the median values. For both equations, the values for ∆T opt−ir and ∆T B are to be subtracted from T opt to determine T ef f . Following Flores et al., we adopted T opt = 3670 K for LkCa 4 from Herczeg & Hillenbrand (2014). For the first relationship, we find a value of T ef f = 3520 K and correspnding to Log (L/L ) = −0.144. For the second relationship, we assumed the median B-field strength of 1.8 kG reported by Flores et al., likely to be a conservative value for LkCa 4, and find T ef f = 3435 K with Log (L/L ) = −0.187. As in all other cases, the associated values of L were calculated using R = 2.3 R . Given the large uncertainties on the empirical relations, we are somewhat surprised to see such reasonable agreement between the values these relations predict and the parameter space we find using our two-temperature models. Figure 9 gives a closer view of the corrected placement of LkCa 4, the GS17 values, and those derived from the Flores+22 empirical relations with respect to both SPOTS and standard evolutionary models.
For a subset of stars in the Flores et al. (2022) sample, dynamical masses were available from Atacama Large Millimeter/submillimeter Array observations and were compared to the masses inferred from both optical and infrared temperatures using the Feiden (2016) evolutionary models which incorporate magnetic effects. They find that the masses inferred from infrared temperatures are much closer to the dynamical masses of the stars than masses inferred from the optical temperatures. However, for the lowest mass stars in their sample, 0.25-0.4 M , the infrared masses are also overestimations when compared to the dynamical masses by as much as 50%. The interesting and potentially paradigmshifting implications of this work is that the "spot" component, F λ (T spot ), of the two-temperature models appears to be a better approximation of the star's spectral type than the photospheric component, F λ (T phot ). Per- Table 7. Spot-corrected T ef f and L Values
The results of these studies strongly suggest that previous determinations of temperatures and luminosities for many young stars are inherently flawed since the optical spectra appear to be dominated by small, warm, non-representative regions of the stellar surface. For such stars, the masses and ages inferred from the warmer temperatures will represent overestimates of masses and ages alike, whether compared to standard evolutionary models or those that take into account magnetic effects and spots. Whether or not these results hint at a new physical phenomenon associated with these stars is unclear, but for sources with large discrepancies between optical and infrared colors and temperatures, the spread in masses and ages for objects in a given cluster will be quite large if left uncorrected.

Impact of Short-and Long-term Variability on Mass and Age Estimates
The short-term, periodic variability of the filling factors due to stellar rotation and the long-term variability associated with changes in the total spot coverage both impact the placement of a star on the H-R diagram and the SPOTS models chosen for comparison. The parameter space enclosed within the solid parallelograms in both Figures 8 and 9 accounts for the full range of variability of the instantaneous filling factors observed over one full rotation of the star. The min-to-max amplitude of the short-term variability in the filling factors is ∼0.1. In the long-term V -band variability inferred from the AAVSO data ( Figure 5) spanning nearly a decade and several thousand stellar rotations, we find a comparable level of variability in the total spot filling factor. The parameter space bounded by the parallelogram also constrains the range of ages and masses caused by long- Figure 8. Plotted are the spot-corrected T ef f and L values associated with the K5Vp+M5Vs (orange circle) and K7Vp+M5Vs (blue circle). Error bars on the values come from assuming a half-subclass uncertainty in T phot and Tspot and from the uncertainty in the mean fspot. Overlaid are isochrones (dashed) and evolutionary tracks (solid) from Baraffe et al. (2015) (orange) and Somers et al. (2020) (black) with fspot = 0.85. The gray parallelogram encloses the full range of T ef f and L values from the instantaneous filling factors with uncertainties. The solitary black triangle represents the values from Donati et al. (2014). The black cross at T = 3670 K is from Herczeg & Hillenbrand (2014). The green (∆Topt−ir) and cyan (∆TB) diamonds represent the T ef f values derived from Flores+2022 empirical relations. The red and purple squares are the "best" (T phot = 4100 K, Tspot = 2750 K, fspot = 0.80) and "alternate" (T phot = 3100 K, Tspot = 3000 K, fspot = 0.80) estimates from GS17. term variations in spot coverage. Therefore, while the application of the SPOTS models is affected by changes to the spot coverage on the order of years, LkCa 4 suggests that over the period of the last decade these models would constrain the mass uncertainties to be ∆M = 0.10 M and age uncertainties to be ∆Age = 0.625 Myr.

SUMMARY & CONCLUSION
We illustrate the utility of multi-epoch, mediumresolution spectroscopy combined with two-temperature empirical composite models to accurately constrain (1) the photospheric and spot temperatures, (2) the visual extinction, and (3) the instantaneous and total spot filling factors for the heavily-spotted young star LkCa 4. Relying predominantly on the 0.8-1.35 µm region of the SpeX spectra that is sensitive to interstellar extinction and possesses absorption bands associated with TiO and FeH, we fit spectral models of spotted stars with four different photospheric temperatures, eight spot temperatures, and 11 visual extinctions to each of the 12 observations of LkCa 4 collected over five consecutive nights allowing f spot to vary freely between 0.0 and 1.0. Mini-mizing χ 2 over all possible models, we find two best-fit composite spectra, K5V p +M5V s and K7V p +M5V s , and A V = 0.3 consistently provide the best fits. Night-tonight variations in the filling factors positively correlate with the historic AAVSO light curves for the system with a rotational period of 3.374 days. Such a correlation demonstrates how multi-epoch spectroscopic observations with moderate resolution can detect the rotation of a spotted star and better constrain the total spot coverage than a single observation.
In addition, we have used the stellar parameters, T phot , T spot , and f spot returned by the best-fit composite models to predict the magnitudes of V -band variability observed in a similar time frame to when the spectroscopic data was collected and found good agreement. Whether we think of the star as possessing the warmer photosphere with a significant fraction of its surface mottled with cooler regions of suppressed convection activity or being a cooler star with warm spots, the two-temperature models simultaneously explain the observed spectroscopic and photometric variability of the source. Regardless, the observed correlation between magnetic field strengths, anomalous colors, and optical and infrared spectral type mismatches hint at a magnetic origin to these phenomena (Flores et al. 2022).
Assuming that F λ (T p ) and F λ (T s ) remain fairly constant over year to decade-long timescales, the 7yr time period of AAVSO data studied indicates that the total spot filling factor does not change by more than 5-10% over this time frame. Small shifts in the rotational phase over year-long time intervals indicate possible migration and evolution of the spots on this time frame.
In comparing the placement of the optical temperatures and luminosities and the corrected values to standard evolutionary models on the H-R diagram, we infer significantly lower masses and younger ages for LkCa 4; a result that agrees well with those from GS17. When compared to the SPOTS models, we find some of the shift to lower masses and younger ages is mitigated. However, the shifts are still significant and important in light of the spread in ages frequently observed for star forming regions. The range of the corrected values of T ef f and L associated with single-epoch observations alone shows considerable spread in mass, ∆M = 0.2 M and ∆ Age ≈ 1 Myr and makes the case for more multiepoch studies of young spotted stars.
Given the apparent ubiquity of optical vs. infrared color and spectral type discrepancies for young, lowmass stars, characterizing the impact on previous measurements of T ef f and L will be useful for revising and, possibly, refining age and mass inferences for most young stellar clusters. In general, these results should also problematize the notion that a single spectral type can be assigned to most young stars. For instance, the D CO index points to a T ef f that more closely aligns with that determined from optical and NIR TiO features. This value is a few hundred Kelvin warmer than the corrected values from the two-temperature empirical composite models, but several hundred Kelvin cooler than the a spectral type based on both optical spectra and color indexes.
The work of Flores et al. (2022) shows that masses inferred from infrared temperatures much more closely align with the dynamical masses of the star, yet they too are likely overestimates for the lowest mass stars in their sample. Their work highlights the need for more comparisons between mass inferences involving evolutionary models and dynamically determined stellar masses. Nonetheless, these results intriguingly point toward a fundamental shift in the phenomenological understanding of the color and spectral type discrepancies.
A collective effort by the star forming community combining multiwavelength and multi-epoch observations of sources in several nearby star forming regions would help to characterize total spot coverage, spot and photospheric temperatures, and determine spot-corrected values for T ef f and L across the optical and infrared. Such studies would also permit the application of the newer SPOTS evolutionary models to an extensive sample of PMS stars. As illustrated by Flores et al. (2022), important to these efforts will be the availability of dynamical masses, which will provide useful and necessary constraints on the models as well as additional confirmation or refutation of the phenomenon of large cool starspots. The availability of medium-resolution NIR spectrographs on 3 and 4 meter class telescopes makes this methodology an attractive approach to determining spot characteristics for a large number of nearby PMS.
The authors wish to thank Kevin Covey and Christian Flores for early and insightful discussions related to this work. We also thank spot expert Steve Saar for his insights related to the nature of spots and the magnetic effects on the spot indicators. We thank the staff at NASA's IRTF for generously supporting the observations presented herein. We thank the anonymous referee for helpful feedback that has certainly improved this work. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. H.M. acknowledges support from the National Science Foundation through the Keck Northeast Astronomy Consortium's REU program through grant AST-1950797. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. The authors also acknowledge a publication grant from Colgate Univerisity's Research Council. RED VALUES In order to determine the best-fit models and values for F λ (T phot ), F λ (T spot ), f spot , and A V , we probed a large parameter space spanning four photospheric templates, eight spot templates, 11 visual extinctions, and filling factors freely varying between 0.0 and 1.0. In Tables A1 & A2, we report the best-fit f spot and associated χ 2 red values determined from the three spot indicators for each of the 12 observations. Values are listed for the two best photospheric templates K5V p and K7V p combined with six of the eight spot templates. All have been calculated using A V = 0.3.
For both the F λ (K5V p ) and F λ (K7V p ) model fits, we find that the warmest (M1V and M2V) and coolest (M7V and M8V) spot templates produce the largest χ 2 red values and the poorest fits for both TiO and F 0.8−1.35µm spot indicators. Models constructed with M4V s , M5V s , and M6V s spot templates produce substantially better fits. For each of the six possible composite models and each of the three spot indicators, we calculated the mean χ 2 red values over the 12 observations. The K5V p +M5V s and K7V p +M5V s models possess both the lowest mean χ 2 red values (with the value for the K5V p models being marginally better than the K7V p ) and the smallest standard deviations from the mean suggesting these models are consistently the best fits to the 12 observations spread over five nights.
As discussed in 3.1, the FeH indicator does not easily distinguish between the best-fit models, selecting models with warm spots and extremely large filling factors (i.e., essentially single-temperature fits) and those with cooler spots and smaller filling factors equally over nearly all the epochs. In addition, the poorer fits made to the models using an M4V standard to model the spots seem to be an outlier and may reveal something about the nature of the standard (see Figure 3).