Placing the spotted T Tauri star LkCa 4 on an HR diagram

TL;DR

This study models the impact of starspots on the observed properties of the T Tauri star LkCa 4, revealing that starspots significantly affect age and mass estimates derived from stellar evolution models.

Contribution

Introduces a new spectral fitting methodology to constrain starspot properties and highlights the need for spotted stellar evolution models for accurate young star characterization.

Findings

Starspots cover about 80% of LkCa 4's surface.

Starspots cause the star to appear younger and less massive.

Spectral features indicate a hot and a cool component with distinct temperatures.

Abstract

Ages and masses of young stars are often estimated by comparing their luminosities and effective temperatures to pre-main sequence stellar evolution tracks, but magnetic fields and starspots complicate both the observations and evolution. To understand their influence, we study the heavily-spotted weak-lined T-Tauri star LkCa 4 by searching for spectral signatures of radiation originating from the starspot or starspot groups. We introduce a new methodology for constraining both the starspot filling factor and the spot temperature by fitting two-temperature stellar atmosphere models constructed from Phoenix synthetic spectra to a high-resolution near-IR IGRINS spectrum. Clearly discernable spectral features arise from both a hot photospheric component $T_{\mathrm{hot}} \sim4100$ K and to a cool component $T_{\mathrm{cool}} \sim2700-3000$ K, which covers $\sim80\%$ of the visible surface. This mix of hot and cool emission is supported by analyses of the spectral energy distribution, rotational modulation of colors and of TiO band strengths, and features in low-resolution optical/near-IR spectroscopy. Although the revised effective temperature and luminosity make LkCa 4 appear much younger and lower mass than previous estimates from unspotted stellar evolution models, appropriate estimates will require the production and adoption of spotted evolutionary models. Biases from starspots likely afflict most fully convective young stars and contribute to uncertainties in ages and age spreads of open clusters. In some spectral regions starspots act as a featureless veiling continuum owing to high rotational broadening and heavy line-blanketing in cool star spectra. Some evidence is also found for an anti-correlation between the velocities of the warm and cool components.

Tables

Table 1. Table 1 Estimated B ​ V ​ R 𝐵 𝑉 𝑅 BVR magnitudes

JD $-$ 2450000	Instrument	$B$	$V$	$R$
781.7106	2MASS	14.02	12.61	11.25
4830.5637	TripleSpec	$\cdots$	13.02	$\cdots$
4830.7719	DoubleSpec	$\cdots$	12.94	$\cdots$
6665.7204	ESPaDOnS	14.36	12.87	11.71
6666.8505	ESPaDOnS	14.15	12.69	11.65
6667.7727	ESPaDOnS	14.31	12.87	11.63
6668.8699	ESPaDOnS	14.44	12.94	11.75
6672.8995	ESPaDOnS	14.25	12.80	11.56
6673.8408	ESPaDOnS	14.12	12.67	11.62
6674.7746	ESPaDOnS	14.46	12.99	11.85
6675.7396	ESPaDOnS	14.40	12.90	11.73
6676.7954	ESPaDOnS	14.18	12.72	11.65
6677.8699	ESPaDOnS	14.29	12.86	11.61
6678.7419	ESPaDOnS	14.53	13.02	11.80
6678.8950	ESPaDOnS	14.47	12.97	11.75
7344.8610	IGRINS	14.24	12.84	11.38

Table 2. Table 2 Adopted values of LkCa 4 from IGRINS spectra

$T_{\text{hot}}$	$T_{\text{cool}}$	$f_{\Omega }$
K	K	%
$4100\pm 100$	${2750}_{-50}^{+250}$	$80\pm 5$

Table 3. Table 3 Results from one- and two- component IGRINS spectral fitting

			One comp.	Two components
Order	Instrument	${\lambda }_1-{\lambda }_2$	$T_{\mathrm{eff}}$	$T_{\mathrm{hot}}$	$T_{\mathrm{cool}}$	$f_{\mathrm{cool}}$
		$\mu$m	K	K	K
0	ESPaDoNs	5161$-$5200	${3734}_{-63}^{+183}$	$\cdots$	$\cdots$	$\cdots$
1	ESPaDoNs	5199$-$5263	${3991}_{-81}^{+70}$	$\cdots$	$\cdots$	$\cdots$
2	ESPaDoNs	5262$-$5328	${3990}_{-140}^{+196}$	$\cdots$	$\cdots$	$\cdots$
3	ESPaDoNs	5327$-$5395	${3923}_{-77}^{+72}$	$\cdots$	$\cdots$	$\cdots$
4	ESPaDoNs	5394$-$5465	${3812}_{-85}^{+54}$	$\cdots$	$\cdots$	$\cdots$
5	ESPaDoNs	5464$-$5537	${4165}_{-61}^{+59}$	$\cdots$	$\cdots$	$\cdots$
6	ESPaDoNs	5536$-$5612	${4129}_{-38}^{+37}$	$\cdots$	$\cdots$	$\cdots$
7	ESPaDoNs	5611$-$5689	${4057}_{-52}^{+64}$	$\cdots$	$\cdots$	$\cdots$
8	ESPaDoNs	5688$-$5770	${4038}_{-44}^{+50}$	$\cdots$	$\cdots$	$\cdots$
9	ESPaDoNs	5769$-$5854	${4212}_{-64}^{+71}$	$\cdots$	$\cdots$	$\cdots$
10	ESPaDoNs	5853$-$5941	${4119}_{-46}^{+52}$	$\cdots$	$\cdots$	$\cdots$
11	ESPaDoNs	5940$-$6032	${4091}_{-63}^{+83}$	$\cdots$	$\cdots$	$\cdots$
12	ESPaDoNs	6031$-$6127	${4108}_{-42}^{+46}$	$\cdots$	$\cdots$	$\cdots$
13	ESPaDoNs	6126$-$6225	${3941}_{-24}^{+28}$	$\cdots$	$\cdots$	$\cdots$
14	ESPaDoNs	6225$-$6329	${4012}_{-42}^{+33}$	$\cdots$	$\cdots$	$\cdots$
15	ESPaDoNs	6328$-$6436	${4238}_{-40}^{+49}$	$\cdots$	$\cdots$	$\cdots$
16	ESPaDoNs	6435$-$6549	${3880}_{-76}^{+61}$	$\cdots$	$\cdots$	$\cdots$
17	ESPaDoNs	6548$-$6667	${3874}_{-70}^{+72}$	$\cdots$	$\cdots$	$\cdots$
18	ESPaDoNs	6666$-$6791	${4054}_{-32}^{+42}$	$\cdots$	$\cdots$	$\cdots$
20	ESPaDoNs	6920$-$7057	${3907}_{-20}^{+32}$	$\cdots$	$\cdots$	$\cdots$
21	ESPaDoNs	7056$-$7201	${3910}_{-20}^{+25}$	$\cdots$	$\cdots$	$\cdots$
22	ESPaDoNs	7200$-$7352	${3878}_{-15}^{+28}$	$\cdots$	$\cdots$	$\cdots$
23	ESPaDoNs	7351$-$7512	${4085}_{-66}^{+51}$	$\cdots$	$\cdots$	$\cdots$
26	ESPaDoNs	7858$-$8051	${3852}_{-24}^{+24}$	$\cdots$	$\cdots$	$\cdots$
29	ESPaDoNs	8473$-$8707	${3504}_{-46}^{+49}$	$\cdots$	$\cdots$	$\cdots$
30	ESPaDoNs	8706$-$8954	${3821}_{-29}^{+33}$	$\cdots$	$\cdots$	$\cdots$
123	IGRINS	14519$-$14721	$\cdots$	${4182}_{-75}^{+107}$	${2773}_{-68}^{+254}$	${0.82}_{-0.72}^{+0.03}$
122	IGRINS	14634$-$14837	$\cdots$	${4185}_{-137}^{+187}$	${2819}_{-105}^{+209}$	${0.78}_{-0.22}^{+0.04}$
121	IGRINS	14750$-$14956	$\cdots$	${4084}_{-144}^{+153}$	${2762}_{-58}^{+83}$	${0.79}_{-0.05}^{+0.04}$
120	IGRINS	14868$-$15076	$\cdots$	${4027}_{-112}^{+76}$	${2739}_{-37}^{+94}$	${0.69}_{-0.08}^{+0.05}$
119	IGRINS	14988$-$15198	$\cdots$	${4093}_{-132}^{+123}$	${2719}_{-18}^{+67}$	${0.76}_{-0.05}^{+0.04}$
118	IGRINS	15111$-$15322	$\cdots$	${4124}_{-131}^{+142}$	${2730}_{-30}^{+83}$	${0.74}_{-0.08}^{+0.08}$
117	IGRINS	15235$-$15448	${4034}_{-122}^{+147}$	${4105}_{-73}^{+63}$	${2731}_{-29}^{+70}$	${0.73}_{-0.06}^{+0.05}$
116	IGRINS	15362$-$15577	${4239}_{-73}^{+55}$	${4121}_{-146}^{+147}$	${2795}_{-85}^{+169}$	${0.80}_{-0.09}^{+0.04}$
115	IGRINS	15490$-$15707	${3718}_{-130}^{+94}$	${3914}_{-151}^{+176}$	${2769}_{-64}^{+298}$	${0.55}_{-0.44}^{+0.14}$
114	IGRINS	15622$-$15841	$\cdots$	${4021}_{-122}^{+159}$	${2751}_{-43}^{+147}$	${0.64}_{-0.13}^{+0.10}$
113	IGRINS	15755$-$15976	$\cdots$	${3874}_{-170}^{+265}$	${2779}_{-75}^{+193}$	${0.77}_{-0.11}^{+0.08}$
112	IGRINS	15891$-$16114	$\cdots$	${3986}_{-206}^{+268}$	${2943}_{-198}^{+425}$	${0.53}_{-0.17}^{+0.16}$
111	IGRINS	16030$-$16255	${4101}_{-116}^{+100}$	$\cdots$	$\cdots$	$\cdots$
110	IGRINS	16171$-$16398	${4233}_{-104}^{+59}$	${4180}_{-193}^{+264}$	${2777}_{-68}^{+184}$	${0.83}_{-0.17}^{+0.06}$
109	IGRINS	16314$-$16543	${3933}_{-97}^{+126}$	$0_{-0}^{+0}$	$0_{-0}^{+0}$	${0.50}_{-0.00}^{+0.00}$
108	IGRINS	16461$-$16692	${4150}_{-176}^{+100}$	${4248}_{-193}^{+183}$	${2802}_{-84}^{+339}$	${0.71}_{-0.35}^{+0.12}$
107	IGRINS	16610$-$16843	${3750}_{-112}^{+186}$	${3813}_{-177}^{+521}$	${3237}_{-486}^{+437}$	${0.38}_{-0.36}^{+0.27}$
106	IGRINS	16762$-$16997	${4225}_{-92}^{+68}$	${4324}_{-168}^{+161}$	${2882}_{-175}^{+1227}$	${0.53}_{-0.53}^{+0.19}$
105	IGRINS	16917$-$17155		${3930}_{-152}^{+209}$	${2763}_{-59}^{+212}$	${0.66}_{-0.10}^{+0.06}$
104	IGRINS	17075$-$17315		${4300}_{-245}^{+177}$	${2870}_{-130}^{+115}$	${0.83}_{-0.05}^{+0.05}$
103	IGRINS	17236$-$17478		${3848}_{-334}^{+312}$	${2972}_{-234}^{+199}$	${0.65}_{-0.30}^{+0.18}$
102	IGRINS	17400$-$17645		${4336}_{-159}^{+136}$	${2751}_{-48}^{+113}$	${0.82}_{-0.05}^{+0.04}$
101	IGRINS	17568$-$17815		${4199}_{-574}^{+264}$	${2916}_{-158}^{+158}$	${0.79}_{-0.23}^{+0.05}$
100	IGRINS	17739$-$17988		${4312}_{-193}^{+161}$	${2765}_{-61}^{+135}$	${0.79}_{-0.06}^{+0.04}$
99	IGRINS	17914$-$18165		$\cdots$	$\cdots$	$\cdots$
94	IGRINS	18855$-$19117	$\cdots$	${4293}_{-1132}^{+185}$	${2899}_{-150}^{+304}$	${0.81}_{-0.39}^{+0.07}$
93	IGRINS	19053$-$19318	$\cdots$	${4044}_{-378}^{+353}$	${2740}_{-36}^{+119}$	${0.81}_{-0.19}^{+0.08}$
92	IGRINS	19256$-$19524	$\cdots$	${3770}_{-183}^{+155}$	${2859}_{-149}^{+510}$	${0.56}_{-0.21}^{+0.18}$
91	IGRINS	19462$-$19734	${3464}_{-67}^{+70}$	${3962}_{-426}^{+389}$	${3126}_{-325}^{+857}$	${0.68}_{-0.56}^{+0.18}$
90	IGRINS	19674$-$19948	${3594}_{-70}^{+75}$	${3907}_{-211}^{+251}$	${3463}_{-161}^{+90}$	${0.75}_{-0.28}^{+0.10}$
89	IGRINS	19890$-$20168	${3557}_{-72}^{+82}$	${3706}_{-151}^{+254}$	${3412}_{-617}^{+183}$	${0.45}_{-0.23}^{+0.20}$
88	IGRINS	20112$-$20392	${3500}_{-60}^{+61}$	${3917}_{-129}^{+208}$	${2812}_{-99}^{+250}$	${0.75}_{-0.06}^{+0.09}$
87	IGRINS	20338$-$20622	${3314}_{-90}^{+77}$	${3709}_{-207}^{+305}$	${2848}_{-119}^{+167}$	${0.75}_{-0.24}^{+0.09}$
86	IGRINS	20570$-$20857	${3544}_{-48}^{+51}$	${3915}_{-275}^{+475}$	${3135}_{-345}^{+281}$	${0.68}_{-0.18}^{+0.15}$
85	IGRINS	20808$-$21098	${3315}_{-107}^{+99}$	${4220}_{-322}^{+252}$	${2828}_{-112}^{+205}$	${0.85}_{-0.10}^{+0.05}$
84	IGRINS	21051$-$21345	${3475}_{-56}^{+62}$	${3681}_{-191}^{+444}$	${3414}_{-130}^{+125}$	${0.68}_{-0.12}^{+0.17}$
83	IGRINS	21300$-$21597	${3234}_{-108}^{+81}$	${4334}_{-235}^{+144}$	${2896}_{-137}^{+166}$	${0.90}_{-0.03}^{+0.03}$
82	IGRINS	21555$-$21856	${3589}_{-8}^{+2}$	${3535}_{-258}^{+582}$	${2964}_{-241}^{+297}$	${0.62}_{-0.38}^{+0.21}$
81	IGRINS	21817$-$22121	${3466}_{-87}^{+117}$	${3536}_{-124}^{+146}$	${3373}_{-150}^{+152}$	${0.34}_{-0.10}^{+0.16}$
80	IGRINS	22085$-$22393	${3157}_{-124}^{+132}$	${3543}_{-273}^{+404}$	${2797}_{-88}^{+336}$	${0.69}_{-0.21}^{+0.16}$
79	IGRINS	22360$-$22671	${3437}_{-60}^{+79}$	${3511}_{-92}^{+138}$	${3383}_{-120}^{+122}$	${0.30}_{-0.10}^{+0.13}$
78	IGRINS	22643$-$22957	${3197}_{-43}^{+43}$	${3808}_{-231}^{+321}$	${2727}_{-26}^{+159}$	${0.78}_{-0.06}^{+0.07}$
77	IGRINS	22932$-$23251	${3273}_{-49}^{+52}$	${4221}_{-392}^{+254}$	${2716}_{-15}^{+273}$	${0.85}_{-0.12}^{+0.03}$
76	IGRINS	23230$-$23552	${3321}_{-53}^{+54}$	${3428}_{-116}^{+317}$	${3123}_{-362}^{+193}$	${0.46}_{-0.21}^{+0.22}$
75	IGRINS	23535$-$23861	$\cdots$	${3649}_{-218}^{+528}$	${2956}_{-222}^{+394}$	${0.62}_{-0.17}^{+0.17}$
74	IGRINS	23849$-$24179	$\cdots$	${3523}_{-94}^{+133}$	${3053}_{-308}^{+388}$	${0.34}_{-0.19}^{+0.20}$
73	IGRINS	24172$-$24505	$\cdots$	${4119}_{-433}^{+345}$	${3029}_{-294}^{+244}$	${0.80}_{-0.13}^{+0.08}$
72	IGRINS	24503$-$24840	$\cdots$	${3897}_{-184}^{+398}$	${2881}_{-157}^{+425}$	${0.74}_{-0.26}^{+0.11}$

Table 4. Table 4 Previous studies of LkCa 4

Reference	Band(s)	R	SpT	$T_{\mathrm{eff}}$	$\log g$	$A_V$	$v\sin i$	$v_z$	$\log L/L_{\odot }$
		$\lambda /\delta \lambda$		K			km/s	km/s
Herbig et al. (1986)	$V$	100000	K7 V	-	-	0.68	26.1$\pm$2.4	13$\pm$4	0.02
Hartmann et al. (1987)	$U$	100000	-	-	-	-	26.1$\pm$2.4	16.9$\pm$2.6
Downes & Keyes (1988)	$V$	13$\AA$	Me	-	-	-	-	-
Strom et al. (1989a)	IRAS	-	K7:V	-	-	0.95	-	-	0.04
Strom et al. (1989b)	$V$	3500	K7:V	-	-	-	-	-
Stauffer et al. (1991)	$VRI$	500	M1	-	-	-	-	-
Strom & Strom (1994)	$V$	-	K7	4000	-	1.25	-	-	0.06
Kenyon & Hartmann (1995)	$U$-IRAS	-	K7	4060	-	0.69	-	-	-0.08
Hartigan et al. (1995)	V	25000	K7	4000	-	0.68	-	-	-0.07
White & Ghez (2001)	$VI$	-	K7	-	-	1.21	-	-	-0.02
Nguyen et al. (2009)	$BVRI$	60000	-	-	-	-	$30\pm 2$	-
Nguyen et al. (2012)	$BVRI$	60000	-	-	-	-	$30\pm 2$	$16.0\pm 4.0$
Grankin (2013)	$BVRI$	-	K7	4040	-	0.54	26.1	-	-0.13
Donati et al. (2014)	$UBVRI$	68000	-	4100$\pm 50$	3.8$\pm$0.1	0.68$\pm$0.15	28.0$\pm$0.5	16.8$\pm$0.1	-0.04$\pm$0.11
Herczeg & Hillenbrand (2014)	$UBVRI$	700	M1.3	3670	-	0.35	-	-	-0.29
This Work: Best	$HK$	45000	-	3180	$\sim$3.8	0.3	29$\pm$3	-	-0.26
This Work: Alternate	$HK$	45000	-	3330	$\sim$3.8	-	29$\pm$3	-	$\sim -0.26$

Equations

$$S_{\mathrm{mix}}={\Omega }_{\mathrm{hot}}B(T_{\text{hot}})+{\Omega }_{\mathrm{cool}}B(T_{\text{cool}})$$

$$f_{\mathrm{cool}}=\frac{{\Omega }_{\mathrm{cool}}}{{\Omega }_{\mathrm{hot}}+{\Omega }_{\mathrm{cool}}}$$

$$T_{\mathrm{eff}}=[T_{\mathrm{hot}}^4(1-f_{\mathrm{cool}})+T_{\mathrm{cool}}^4{f}_{\mathrm{cool}}{]}^{0.25}.$$

$${\bar{f}}_{\lambda }=\frac{f_{\lambda }}{{\int }_0^{\infty }{f}_{\lambda }d\lambda }=\frac{f_{\lambda }}{f}$$

$$f_{\lambda ,\mathrm{mix}}$$

$$q$$

$$f_{\lambda ,\mathrm{mix}}^{\prime }$$

$${\mathsf{M}}_{\mathrm{mix}}^{\prime }={\Omega }_{\mathrm{hot}}({\tilde{\xi }}_{\mu }+\mathbf{X}{\mu }_{\mathbf{w},\mathrm{hot}})+q{\Omega }_{\mathrm{cool}}({\tilde{\xi }}_{\mu }+\mathbf{X}{\mu }_{\mathbf{w},\mathrm{cool}})$$

$${\mathsf{C}}_{\mathrm{mix}}^{\prime }={\Omega }_{\mathrm{hot}}^2\mathbf{X}{\mathbf{\Sigma }}_{\mathbf{w},\mathrm{hot}}{\mathbf{X}}^T+q{\Omega }_{\mathrm{cool}}^2\mathbf{X}{\mathbf{\Sigma }}_{\mathbf{w},\mathrm{cool}}{\mathbf{X}}^T+\mathsf{C}$$

Full text

††slugcomment: Accepted to ApJ

1. Introduction

The causes of age uncertainties and of large luminosity spreads in young clusters are controversial (see reviews by Preibisch, 2012; Soderblom et al., 2014). Absolute ages are usually inferred from comparisons between observed stellar properties and pre-main sequence evolutionary models, which are sensitive to stellar birthlines and the input physics that controls the contraction and internal heating. Within clusters, any real differences in ages between stars may be masked by observational uncertainties in temperature and luminosity measurements (e.g. Hartmann, 2001; Slesnick et al., 2008).

Magnetic activity is a likely source of significant uncertainty in both the models and observations of young low-mass stars. Convection at these young ages generates strong magnetic fields, as measured by Zeeman broadening and polarimetry (e.g. Johns-Krull, 2007; Donati & Landstreet, 2009) and as seen in starspots (e.g. Stauffer et al., 2003; Grankin et al., 2008). Evolutionary models are just now starting to implement prescriptions for how spots and magnetic fields affect the interior structure (MacDonald & Mullan, 2013; Jackson & Jeffries, 2014; Somers & Pinsonneault, 2015a; Feiden, 2016). Stellar evolution models that include these effects can make a coeval 10 Myr population exhibit apparent age spreads of 3$-$10 Myr. Mass estimates may also be biased to lower masses, as measured for some eclipsing binaries (e.g. Stassun et al., 2014; Rizzuto et al., 2016). Observationally, starspots may also be responsible for biases in stellar effective temperatures derived by different methods. For example, the effective temperatures for 3493 young stars measured using the APOGEE spectrograph ($1.5-1.70\;\mu$m at $R=22,500$ Wilson et al., 2010) are offset from and usually cooler than prior measurements made at optical wavelengths (Cottaar et al., 2014). Similarly, near-IR spectral types of young stars have sometimes been evaluated as later than optical spectral types (Bouvier & Appenzeller, 1992; Vacca & Sandell, 2011).

Rotational modulation from starspots is understood as the cause of cyclic variations in weak-lined T Tauri stars (WTTSs) and some classical T Tauri stars (e.g. Vrba et al., 1986; Herbst et al., 1994). Spots have recently been observed with exceptional photometric precision from monitoring surveys targeting exoplanet transits (e.g Harrison et al., 2011; Davenport et al., 2015). The single band light curves from such surveys cannot separate the relative areas and temperatures of the emitting regions, and therefore cannot yield estimates for the filling factor of spots without making assumptions about the relative temperature contrast. Contemporaneous or near-contemporaneous panchromatic photometric monitoring (Herbst et al., 1994; Petrov et al., 1994; Bouvier et al., 1995; Grankin et al., 2007; Cody et al., 2014) can break the degeneracy between starspot filling factor and starspot temperature, but these measurements are still limited to detecting large longitudinally asymmetric starspots. Transiting planets or eclipsing binaries passing in front of starspots (Desert et al., 2011) can break geometric degeneracies, providing estimates for the characteristic sizes and lifetimes of starspots or starspot groups, but the number of such transiting systems and demands of high photometric precision limit wider application. Zeeman Doppler Imaging (Donati et al., 2014, ZDI) produces the latitudinal and longitudinal distribution of starspots in reconstructed brightness maps based on models for spectro-polarimetric line shapes but would be unable to identify uniform distributions of small spots. All of these techniques would tend to underestimate the fractional coverage of starspots distributed isotropically on the surface, stars viewed pole-on, or stars with longitudinally symmetric bands of starspots.

Spectroscopic strategies have the power to measure both the starspot areal coverage and temperature in any geometric configuration. For WTTSs, SED modeling indicated the presence of cool spots (Wolk & Walter, 1996), while spot temperatures have been estimated from rotational modulation of optical colors (e.g. Grankin, 1998; Venuti et al., 2015; Koen, 2016). The TiO bands in optical spectra have been used to measure spots on WTTSs (Petrov et al., 1994), stars with ages of 125 Myr in the Pleiades (Fang et al., 2016), and evolved stars (e.g. Neff et al., 1995; O’Neal et al., 2001, 2004). Recent spectroscopic detections of spots on young stars have focused on TW Hya, an active accretor, and DQ Tau, a close binary in which both components are accreting (Debes et al., 2013; Bary & Petersen, 2014).

The star LkCa 4, a WTTS member of the Taurus Molecular Cloud (Herbig et al., 1986; Strom et al., 1989a; Downes & Keyes, 1988; Strom et al., 1989b) is an ideal exemplar for a spotted pre-MS star because it does not have any mid-IR or mm excess (e.g. Andrews & Williams, 2005; Furlan et al., 2006; Buckle et al., 2015) and is not actively accreting (e.g. Edwards et al., 2006; Cauley et al., 2012). LkCa 4 has no evidence for a close companion from either direct imaging searches (Karr et al., 2010; Kraus et al., 2011; Daemgen et al., 2015)111The status for wide companions is less clear; see Stauffer et al. (1991); Itoh et al. (2008); Kraus & Hillenbrand (2009); Kraus et al. (2011); Herczeg & Hillenbrand (2014) or from RV searches (Nguyen et al., 2012) and high-precision Zeeman Doppler Imaging (Donati et al., 2014). This single star demonstrates a large amplitude of sinusoidal photometric variability indicative of large spots (Grankin et al., 2008; Xiao et al., 2012). The variability amplitude cannot arise from an eclipsing binary since the large RV variations would have been seen in spectroscopic monitoring. All of these observations indicate that the spectrum of LkCa 4 should be devoid of complicating factors like near-IR excess veiling, accretion excess, or close binarity.

LkCa 4 has recently been examined with ZDI (Donati et al., 2014), revealing a complex distribution of cool and warm spots in brightness map reconstructions. Dark polar spots extend to about $30^{\circ}$ from the pole, with about 5 appendages reaching down to about $60^{\circ}$ from the pole. There is also evidence for a hot spot in the reconstructed ZDI map. These spot maps demonstrate that multiple temperature components contribute to the broadband emission from the star, which may explain the large differences between the 4100 K effective temperature measurement obtained from high-resolution optical spectra (Donati et al., 2014) compared with 3600 K obtained from an analysis dominated by TiO bands (Herczeg & Hillenbrand, 2014).

In this work, we measure the starspot properties of LkCa 4 with four complementary techniques.222Some data, Python code, Starfish configuration files, and other documents related to this paper are available on the project’s GitHub repository at https://github.com/BrownDwarf/welter Section 2 describes the full spectroscopic and photometric dataset and quantifies the optical variability of LkCa 4 over the last 31 years based on all available photometric monitoring; a collection of spectral observations at various phases of variability is introduced. Section 3 and Appendix A describe a spectral inference technique aimed at generalizing line depth ratio analysis, which is applied to a high resolution $H$ and $K$ near-IR spectrum of LkCa 4 in Section 3.5. Section 4 shows results from SED fitting, polychromatic photometric monitoring, and optical TiO line-depth ratio fitting. All lines of evidence are ultimately combined to build a consistent picture of the spectral and temporal evolution of LkCa 4. Section 5 describes how the two-temperature fit affects the placement of LkCa 4 on the HR diagram and what can be understood about pre-main sequence stellar evolution from this exemplar. Finally, Section 6 describes how our results should be interpreted, the implications for analyses of spots in high-resolution spectra, and methodological limitations.

2. Observations

2.1. IGRINS Spectroscopy

We acquired observations of LkCa 4 with the Immersion Grating Infrared Spectrograph, IGRINS (Park et al., 2014) on the Harlan J. Smith Telescope at McDonald Observatory on 2015-11-18 $09^{h}$ UTC. IGRINS is a high resolution near-infrared echelle spectrograph providing simultaneous $R\simeq 45,000$ spectra over 1.48-2.48 $\mu$m. The spectrograph has two arms with 28 orders in $H-$band and 25 orders in $K-$band. The data were reduced with the Pipeline Package (Lee, 2015). Telluric correction was performed by dividing the spectrum of LkCa 4 by a spectrum of HR 1237, an A0V star observed immediately before LkCa 4 , both at airmass 1.1. The broad hydrogen lines in the A0V star produced excess residuals in the spectrum of LkCa 4. No effort was made to mitigate these flux excess residuals, which affect several spectral orders. These orders were ultimately excluded from further analyses.

2.2. ESPaDOnS Spectroscopy

We used ESPaDOnS on CFHT to obtain twelve high-resolution optical spectra of LkCa 4 from 8-21 Jan. 2014 as part of the MaTYSSE Large Program. These spectra cover 3900$-$10000 Å at $R\sim 68,000$ and were obtained in spectropolarimetry mode. The Zeeman Doppler Imaging obtained from these observations were analyzed by Donati et al. (2014). In this paper, we concentrate on the intensity spectrum. To calculate TiO indices, the relative flux calibration was obtained from an approximate instrument sensitivity calculated from observations of 72 Tau. The spectral shapes match the flux calibrated low-resolution spectra described in §2.3.

2.3. Low-resolution optical and near-IR spectra

A flux-calibrated low-resolution optical spectrum of LkCa 4 was obtained using Palomar/DBSP (DBSP, Oke & Gunn, 1982) on 30 Dec. 2008 at 06:30 UT. These spectra cover 3200–8700 Å at a spectral resolution $R\sim 500$. The data reduction and a spectral analysis are described by Herczeg & Hillenbrand (2014) within the context of a larger survey.

Near-infrared spectra of LkCa 4 were obtained at  01:30 UT on Dec. 30, 2008, a few hours earlier than the DBSP spectra, using the TripleSpec spectrograph on the Astrophysical Research Consortium (ARC) 3.5 meter telescope at Apache Point Observatory in Sunspot, New Mexico. LkCa 4 was observed for a total integration time of 90 s using an ABBA dither sequence along the slit to remove sky emission by differencing sequential exposures. With the 1.1$\arcsec$ slit used for these observations, TripleSpec provides nearly contiguous coverage from 0.95 to 2.46 $\mu$m at a resolution of R $\sim$3200 (Wilson et al., 2004). Spectra were reduced with a variant of the SpeXTool pipeline, originally developed by Cushing et al. (2004) and modified for use with TripleSpec data. Spectra were differenced, flattened, extracted, and wavelength calibrated prior to telluric correction and flux calibration; the latter two corrections were performed using the XTELLCOR IDL package (Vacca et al., 2003) and a spectrum of HD 25175, a nearby A0V star observed directly after LkCa 4 at a similar airmass.

2.4. Photometric monitoring

Photometry in $B$, $V$, and/or $R$-bands in 1216 distinct epochs was assembled from published and unpublished data sources spanning 31 years in 22 observing seasons. Published data targeting LkCa 4 is composed of 29 $BVR$ visits from 1985$-$1986 (Vrba et al., 1993), 26 $BVR$ epochs from 1990$-$1991 (Bouvier et al., 1993), 284 visits in $V$-band (278 with $B$ and 268 with $R$) from 1992 August - 2004 October (Grankin et al., 2008), and 10 $BVR$ visits from 2013 (Donati et al., 2014). ASAS3 (Pojmański, 2004) acquired 63 $V-$band measurements from 2002$-$2004. Integral-OMC obtained 138 $V-$band measurements from 2006$-$2008 (Alfonso-Garzón et al., 2012). The AAVSO archive (Kafka, 2016) includes 385 $V$, 23 $B$ and 10 $R$ measurements from 2013$-$2016. Unpublished data from the ASAS-SN survey (Shappee et al., 2014) were obtained from 2012 January – 2016 March in 186 visits. Finally, recent photometry from the ongoing Crimean Astrophysical Observatory (CrAO) ROTOR project (Grankin et al., 2008) is presented here for the first time. The CrAO data includes 43 $BVR$ visits from August 2015 - March 2016. Figure 1 displays all the $V-$band photometry over the interval 1985-2016. HJD, BJD, and JD are simply referred to as JD, which is accurate to the precision of available data.

Figure 2 shows all available $V-$band data grouped by the 22 observing seasons and phase-folded by the period $P=3.375$ days obtained from multiterm Lomb-Scargle periodograms (Ivezić et al., 2014). The general appearance of the phase-folded lightcurves does not change with perturbations to the period on the scale of 0.003 days. The estimated $BVR$ magnitudes at the time of the spectral observations are determined from regularized multiterm fits (VanderPlas & Ivezić, 2015, i.e. Fourier series truncated to the first $\sim 4$ components) shown as the solid blue lines in Figure 2. Table 1 lists the estimated $BVR$ photometry and the observing epoch for data from IGRINS, ESPaDOnS, 2MASS (Skrutskie et al., 2006), DBSP, and TripleSpec.

3. FITS TO HIGH RESOLUTION SPECTRA

Starspots possess a spectrum distinct from the ambient photosphere. Approaches like ZDI and line-bisector analysis (e.g. Prato et al., 2008; Donati et al., 2014) target the modulation of stellar photospheric lines as starspots enter and exit the observable stellar disk. In contrast, the technique described in this Section measures the spectrum that emerges from both the hot and the cool surfaces of the photosphere.

Starspot line depth ratio analysis has traditionally been limited to isolated portions of spectrum that possess absorption lines attributable only to cool photospheric components and spectral lines arising from the warmer ambient photosphere (e.g. Neff et al., 1995; O’Neal et al., 2001). The apparent veiling of the lines and their respective temperature dependences can be combined to solve for the cool photosphere and ambient photosphere temperatures and relative areal coverages. This line depth ratio method suffers from the need to identify portions of the spectrum that possess easily identifiable strong lines, and from the need to assemble large atlases of observed spotless spectral templates to which the spotted star spectrum can be compared. In this section, we introduce a generalization of the line depth ratio analysis which employs pixel-by-pixel modeling of all echelle spectral orders. By using all the spectral data, this strategy has the power to constrain starspot properties at relatively low filling factors and to identify weak lines that originate from the cool photosphere.

We took two approaches to characterizing the photospheric temperature from a spotted-star spectrum. First, using the spectral fitting approach introduced in Section 3.1, we separately fit the optical (ESPaDOnS) and near-IR (IGRINS) spectra with distinct, single photospheric component models; these fits are described in Sections 3.2 and 3.3 respectively, and then compared in 3.4. Second, we use a two temperature mixture model to perform a fit to the near-IR spectra in Section 3.5. The near-IR spectral range is preferred over the optical, since the cool photosphere emits most of its flux in the near-IR, enhancing the likelihood of direct detection of emission from the cool component.

3.1. Methodology

Czekala et al. (2015, hereafter C15) developed a modular framework333The open source codebase and its full revision history is available at https://github.com/iancze/Starfish. The experimental fork discussed in this paper is at https://github.com/gully/Starfish. to infer stellar properties from high-resolution spectra. The C15 technique forward models observed spectra with synthetic spectra from pre-computed model grids. The intra-grid-point spectra are “emulated” in a process similar but superior to interpolation since it seamlessly quantifies the uncertainty attributable to the coarsely sampled stellar intrinsic parameters (see Appendix of C15). The forward model includes calibration parameters, line spread functions, and a Gaussian process noise model to account for correlations in the residual spectrum. We employed the Phoenix grid of pre-computed synthetic stellar spectra, which are provided at high spectral resolution over a wide wavelength range at 100 K intervals in stellar photospheric temperature444Single component Phoenix photospheres are indexed by the label “effective temperature”, or $T_{\textrm{eff}}$, which presupposes that the emergent spectrum is comprised of only a single photospheric component. In this work we avoid labeling the photospheric components by $T_{\textrm{eff}}$, and instead, refer to them as simply “photospheric temperatures”. in our range of interest (Husser et al., 2013). The modular framework was altered to accommodate starspot measurements in two ways. First, the single photospheric component was updated to include two photospheric components. Second, the MCMC sampling strategy was altered to accommodate the additional free parameters added by the two component model.

The stellar photosphere is characterized as two photospheric components with a cool temperature $T_{\textrm{cool}}$ and hot temperature $T_{\textrm{hot}}$, with scalar solid angular coverages $\Omega_{\mathrm{cool}}$ and $\Omega_{\mathrm{hot}}$, respectively555For flux calibrated spectra, the total solid angle $\Omega$ can be constrained, but for typical echelle spectrographs only relative values of $\Omega_{\mathrm{cool}}$ and $\Omega_{\mathrm{hot}}$ can be inferred.. The choice of nomenclature depends on the physical origin of the emission region, which is not directly observable without ancillary information. For many applications, it could be reasonably assumed that the cool photospheric component corresponds to a sunspot or sunspot group, and the hot component corresponds to the ambient photosphere. However, a cool ambient photosphere possessing localized hot spots—like plages on the sun—could yield the same or similar two-temperature composite spectrum. In principle, the physical origin could be assessed with detailed comparison of average surface magnetic field strengths through Zeeman line broadening or other magnetic-field sensitive techniques. For now, we index the two photospheric components in the most general way possible, simply cool and hot.

The two photospheric components share the same average intrinsic and extrinsic stellar parameters $v\sin i,\log g,[{\rm Fe}/{\rm H}],v_{z}$. The composite mixture model for observed flux density is:

[TABLE]

where $B(T)$ is the spectral radiance from model spectra. The filling factor of the starspots is:

[TABLE]

The term $f_{\mathrm{cool}}$ represents an instantaneous, observational fill factor seen on one projected hemisphere, not $f_{\rm spot}$, the “ratio of the spotted surface to the total surface areas” (Somers & Pinsonneault, 2015a). In the limit of homogeneously distributed starspots and assuming the cool photosphere arises from starspots, $f_{\mathrm{cool}}$ tends to $f_{\rm spot}$, but in general $f_{\rm spot}$ is not a direct observable since some circumpolar regions of inclined stars will always face away from Earth; $f_{\mathrm{cool}}$ can vary cyclically through rotational modulation, whereas $f_{\rm spot}$ changes secularly through starspot evolution that is not yet understood.

The two component mixture model includes two more free parameters than the standard C15 model, namely $T_{\textrm{cool}}$ and $f_{\mathrm{cool}}$. The addition of these two parameters makes the MCMC sampling much more correlated than it was for single component models. This challenge motivated a second technical change to the spectral inference framework, switching from Gibbs sampling (c.f. C15) to ensemble sampling using emcee (Foreman-Mackey et al., 2013). The practical effect of this switch is that all 14 stellar and nuisance parameters are fit simultaneously in a single spectral order, making stellar parameter estimates ${\bm{\theta}}_{o}$ unique for each spectral order $o$, whereas the C15 strategy had the power to provide a single set of stellar parameters ${\bm{\theta}}$ that was based on all $N_{\rm ord}$ spectral orders.

Interpreting the two-temperature mixture model as arising from the same photosphere requires single stars, or at least double stars with extremely large (>100) luminosity ratios or wide separations. Unresolved single-lined or double-lined binaries could mimic a signal that would be misattributed to starspots. As pointed out earlier, LkCa 4 has no detected companion from direct imaging, and its RV variations are consistent with starspot modulation, so any binary companion would have to have an exceptionally large luminosity (or mass) ratio, rendering it undetectable with our spectral inference methodology with the finite signal to noise available in the IGRINS data. The absence of a detectable disk around LkCa 4 means that the spectral modeling requires no further parameters like veiling or accretion.

Further details of the spectral inference methodology are described in Appendix A.

3.2. Single temperature fitting to the ESPaDOnS spectrum

We performed spectral fitting on an ESPaDOnS spectrum acquired on 2014 January 11 (JD 2456668.9). The spectrum was split into subsets of $N=26$ chunks, corresponding approximately to spectral order boundaries. Fitting was performed separately on each spectral chunk. The spectral emulator was trained on stellar parameters in the range $\log g\in[3.5,4.0]$, $[{\rm Fe}/{\rm H}]\in[-0.5,0.5]$, and photospheric temperatures $T\in[3500,4200]$ K.

Standard MCMC convergence criteria were assessed. Typical spectral chunks show photospheric temperature point estimates in the range $4000\pm 130$ K. Figure 3 displays the point estimates with 5th and $95^{th}$ percentile error bars placed at the central wavelength of each spectral chunk. The spectral chunk from 8473$-$8707 Å is an outlier, showing an estimated photospheric temperature of $\sim 3500$ K.

3.3. Single temperature fitting to the IGRINS spectrum

We selected a subset of 32 of the 54 available IGRINS spectral orders666The prohibitive computational cost limited us from running all the spectral orders. with low telluric absorption artifacts. We performed unique spectral fitting on each order, with free stellar parameters ${\bm{\theta}}=(T_{\textrm{eff}},\log g,[{\rm Fe}/{\rm H}])$. The same analysis procedure described for the ESPaDOnS spectra is used here, though in the $K$-band the temperature search range is expanded to $T\in[3000,4200]$ K.

We find a larger dispersion in the point estimates for the stellar parameters derived from the IGRINS data than those derived in the optical spectrum. The most conspicuous trend is in the derived photospheric temperature as a function of wavelength shown in Figure 3. The photospheric temperature peaks at values of $\sim 4200$ K in the short wavelength end of $H-$band and saturates at $<3500$ K at the long wavelength end of $H-$band. The $K-$band shows even lower derived photospheric temperatures of $\sim 3300$ K, or 700 K cooler than estimated from optical. No single photospheric temperature can describe all the spectral lines present in the high resolution optical and IR spectra. Sources with such discrepancies have been seen previously, for example in Figures 4 and 5 in Bouvier & Appenzeller (1992). These discrepancies are circumstantial evidence for the detection of spectral features attributable to a second photospheric component.

3.4. Heightened sensitivity to starspot spectral lines in the infrared

The heightened sensitivity to starspot spectral lines as a function of wavelength can be understood in the following way. Starspots are cooler than their surrounding photosphere and will, therefore, have a longer wavelength of peak emission. The average ratio of flux density between starspot and bulk photosphere will increase with wavelength until asymptotically approaching a fixed value on the Rayleigh-Jeans tail (Wolk & Walter, 1996). The bottom panel of Figure 3 shows the flux density ratio for patches of 2800 K and 4100 K photospheres with equal areas (50% filling factor). The black-body ratios predict smooth flux ratio increases with wavelength, while the ratio of Phoenix models shows wavelength regions with heightened sensitivity to starspot fluxes, for example in $J-$band. The visible portion of the spectrum will be dominated by the warm patches. The infrared spectrum will have a lower overall spot contrast and will have more net starspot flux than optical emission.

Some care should be taken when directly comparing results between the ESPaDOnS and IGRINS spectra since they were not taken at the same time. The particular ESPaDOnS spectrum used in Section 3.2 was obtained when LkCa 4 had an estimated $V-$band brightness of 12.94 mag, compared with 12.84 mag at the time of the IGRINS spectrum. The fainter magnitude during the ESPaDOnS spectrum acquisition implies a greater coverage fraction of cool material than during the IGRINS spectrum acquisition. The photospheric temperature derived in the optical bands and short-wavelength end of $H$-band yield similar values of $\sim 4100$ K, suggesting the bulk777The bulk appearance and disappearance of spectral lines is mostly controlled by the temperature of the emitting region of the local photosphere. In other words, most of the variance in a spectrum is attributable to temperature, assuming near-solar metallicity. Starspots imbue dearths of flux in the line profiles of optical spectra, but these are secondary to the mere presence of temperature-sensitive lines. spectral features are broadly consistent with emission from a single temperature component. A single ESPaDOnS order surrounding the TiO bands shows an exceptionally low estimated $T_{\textrm{eff}}$. The long wavelength portion of $H-$band and all of $K-$band are more sensitive to starspot spectral signatures than the shorter wavelength portions.

3.5. Two-temperature fitting to IGRINS spectra

Full-spectrum fitting was performed for 48 of the 54 available IGRINS spectral orders, omitting only the orders with the most pathological telluric spectral artifacts. We applied the mixture model as described in Section 3.1 and Appendix A. Standard MCMC convergence tests were evaluated. The stellar parameter ranges were $\log g\in[3.5,4.0]$ and $[{\rm Fe}/{\rm H}]\in[-0.5,0.5]$; the Phoenix model spectrum temperature range was $T_{\textrm{hot}},T_{\textrm{cool}}\in[2700,4500]$. $H-$band fits had normal distribution priors in place: solar metallicity to $\pm 0.05$ dex, $logg=3.8\pm 0.1$, and $v\sin i=29\pm 5$ km/s.

The fit quality is first assessed by examining the consistency of the point estimates of $v_{z}$ and $v\sin i$ across the spectral orders. The distribution of $v_{z}$ and $v\sin i$ exposed extremely poor fits in two orders ($o=91$ and $94$), with all other orders demonstrating $v_{z}=12.4\pm 2.6$ km/s and $v\sin i=28.8\pm 2.0$ km/s.

Overplotting forward-modeled spectra with the observed IGRINS spectra yielded insights on why some spectral orders perform better than others in assessing stellar properties. Orders with extremely poor telluric correction residuals, large spectral line outliers, and uncorrected H line residuals from A0V standard division, are all excluded from our final stellar parameter compilation. Several orders in $K-$band were rejected due to poor spectral fits influenced by deep metal lines. Metal lines can be biased by Zeeman broadening, which could explain poorer fit quality with wavelength (Johns-Krull, 2007; Deen, 2013). Some orders were also discarded because their fits were mediocre or uninformative. The remaining orders are relatively devoid of spectral line outliers and include the most information rich portion of the spectrum.

The IGRINS spectrum demonstrates some features that are present only in the hot photosphere model, and some features that are present only in the cool photosphere model. Figure 4 shows a selection of six such spectral features for a range of plausible fill factors. In some cases (e.g. the top two panels), featureless cool photosphere spectra veil isolated spectral lines predicted in the hot photosphere models. In other cases (e.g. the middle row and lower left panel), the hot photosphere model veils sequences of shallow spectral features predicted in the cool photosphere model. The cool spectral models predict shallow features because line blanketing from multiple indistinct molecular bands overlaps from rotational broadening; any feature of interest will be biased to non-zero veiling.

The combination of line blanketing and high $v\sin i$ has probably hampered efforts to identify isolated spectral features suitable for line-depth-ratio analysis in the spectra of rotationally broadened young stars. One such feature, the OH 1.563 $\mu$m line analyzed by O’Neal et al. (2001), shows a clear pattern of 3 lines in our data, with the central line exceeding the depth of the adjacent two lines. The Phoenix models predict a non-negligible hot photosphere contribution to the middle line for our range of hot and cool temperatures. In comparison, the forward-modeling technique described in this work thrives in the presence of long sequences of indistinct yet predictably correlated spectral features. The level of veiling of hot photospheric lines is set by the cool photosphere filling factor and temperature, while the level of veiling of the cool lines is set by the filling factor and temperature of the hot photosphere lines. Similar forward-modelling strategies have successfully identified patterns of weak metal lines in the line-blanketed spectra of low metallicity stars (Kirby, 2011; Kirby et al., 2015). Figures 14 and 15 in the Appendix show 42 of the $H-$ and $K-$band spectra on a log scale with a single random composite model spectrum overplotted.

Figure 5 shows the distribution of $T_{\textrm{hot}}$ (blue shading), $T_{\textrm{cool}}$ (red shading), and $f_{\mathrm{cool}}$ (yellow shading) for all of the spectral orders, with the rejected orders grayed out, and the reliable order subset shown in bold. The point estimates for each spectral order are listed in Appendix Table 3.

Best fit values for the stellar parameters are listed in Table 2. Remarkably, the filling factor of cool photosphere exceeds 50%, with a best fit value of $f_{\mathrm{cool}}=80\pm 5\%$.

4. The Two Temperature Fit to the SED and Stellar Rotation

In the previous section, we established that the high resolution optical and near-IR spectra of LkCa 4 may be explained by a two-temperature photosphere. In this section, we test the two component fit by comparing the combined spectrum to observed spectral features and rotational modulation. We adopt the best-fit two temperature model to the IGRINS spectrum, with components of 2750 K covering 80% of the visible stellar surface and 4100 K covering 20% of the stellar surface, obtained when LkCa 4 had an estimated brightness of $V=12.84$ mag. The cool and hot temperatures and the filling factor are blindly adopted without any adjustments to attempt to improve fits to the SED or broadband spectra.

Figures 6–7 compare synthetic spectra from the two-component photosphere to the SED and flux-calibrated spectra of LkCa 4. The only free parameters in this comparison, the luminosity and extinction, are both scaled to match the spectrum (see §5.1-5.2). The optical-IR SED is obtained by estimating photometry from Grankin et al. (2008) during the 2MASS epoch, with $V\sim 12.61$ mag (see Table 1). The SED also includes mid-IR photometry from Spitzer/IRAC (Hartmann et al., 2005), without adjusting for epoch. In the spectroscopic comparison, some minor discrepancies between the optical and near-IR spectra may be introduced because the spot coverage changed in the $\sim 5$ hrs between observations ($\Delta V=0.08$ mag). The synthetic spectrum is obtained from the Phoenix models, as in §3, and is extended beyond the longest wavelength (5 $\mu$m) to calculate the bolometric luminosity.

The synthetic models match the full spectrum reasonably well. When scaled separately to the red-optical and near-IR spectra, the synthetic models match even better (bottom panels of Figure 7), while single temperature photospheres fail to reproduce large spectral features. The $3900$ K component fits the TiO bands at $<7400$ Å but fails to reproduce molecular bands at longer wavelengths. The $3500$ K component suffers from the opposite problem, yielding molecular bands at short wavelengths that are too deep. The two temperature photosphere with the parameters from the best-fit calculated in §3 reproduces the TiO bands, the hump in the $H$-band, and the jump in flux at the long-wavelength end of the $J$-band. This latter feature is caused by opacity in H2O bands and is only seen when the fraction of emission from the cool component is high (Figure 8).

These spectral features should be rotationally modulated as the filling factors of the visible hot and cool components change. In §4.1, we first use our fit to the IGRINS spectrum to convert the $V$-band brightness to the filling factor of the two components. We subsequently calculate the color changes expected based on this spot coverage and compare these results to historical data (§4.2), and confirm our basic approach by demonstrating that TiO band depths vary with rotation (§4.3).

4.1. The $V$-band lightcurve and spot coverage

The $V$-band brightness reflects the instantaneous filling factor of the cool and hot components. In the two-temperature photosphere, the $V$-band emission is dominated by the hotter component and is, therefore, a good proxy for the visible surface area of the hot component. Figure 9 shows the winter 2015–2016 ASAS-SN lightcurve, converted into a cool spot filling factor. During this period, the brightest and faintest phase corresponds to cool component fill factors of $74\%$ and $86\%$, respectively. Since 1992 the $V$ band photometry has been as bright as $12.3$ mag, which corresponds to a cool component fill factor of $65\%$, and as faint as $13.2$ mag, or a cool component fill factor of $87\%$. This drastic change in visible area of the hot spot (35% to 13%) is needed to explain the full $\sim 1$ mag range in brightness, assuming no temperature change in either component.

4.2. Rotational Modulation of Colors

The relative contributions from the hot and cool components change as LkCa 4 rotates, leading to color changes. We model this rotational modulation using the main sequence colors (compiled by Pecaut & Mamajek, 2013), which should have less contribution from spots than pre-main sequence colors. The Cousins $R$ photometry from Pecaut & Mamajek (2013) is converted to Johnson $R$ for comparison with the Grankin et al. (2008) photometry, following the color transformation prescribed by Landolt (1983) and applied by Grankin et al. Photometry from synthetic spectra produced by Phoenix and BT-Settl models cannot be used here because the predicted $B-V$ colors are much bluer than observed. Offsets between model and observed colors are at least in part attributed to extinction and are discussed in §5.1.

Limb darkening effects also become important if starspots are distributed in a finite number of large patches, not many small patches. The variations in the mean projection angle cause confounding color trends in the photometric modulations.

The optical emission is dominated by the hotter component while both components contribute equally to the infrared emission (see inset in the top panel of Figure 7). Figure 10 demonstrates that the $B-V$ color from Grankin et al. (2008) does not depend on $V$, which is consistent with expectations for a single hot component with no contributions from the cooler component. The standard deviation of 0.03 mag in $B-V$ color (which includes a $\sim 0.01$ mag uncertainty in both $V$ and $B$) indicates that the spot temperature is stable to $\lesssim 50$ K over decades. The stability of this temperature implies that the hotter component may be the ambient temperature of the photosphere.

The correlation between $V$ and $V-R$ indicates that the star becomes redder when the cool spot has a higher filling factor. Our simple model predicts that the correlation should be much weaker than is measured. Most likely, our two temperature model for the photosphere is overly simplistic, since the $V-R$ color could be reproduced with some contribution from intermediate temperatures.

Rotational modulation is expected to lead to much smaller brightness changes at near-IR wavelengths, with amplitudes of $0.15-0.2$ mag, than the $0.6-0.8$ mag amplitudes seen at optical wavelengths. The smaller amplitude of near-IR brightness is also seen in a few spotted stars in the optical/IR monitoring of NGC 2264 (Cody et al., 2014), although those stars have a much smaller $V$ band amplitude than LkCa 4.

4.3. Rotational Modulation of TiO Band Depths

The depths of TiO and other molecular bands are common diagnostics of spectral type and therefore temperature in optical spectra (e.g. Kirkpatrick et al., 1991). In a 2-temperature photosphere cool enough for molecules to form, the depth of observed TiO bands will depend on the fractional coverage of the components. The TiO band depths vary with spot coverage using 12 epochs of CFHT/ESPaDOnS spectra obtained over 14 days ($\sim 4$ rotation periods) by Donati et al. (2014).

The spectra obtained near the estimated maximum and minimum optical brightness are shown in the top panel of Figure 11. Small but significant changes are seen in the TiO bands. On the other hand, the blue spectrum does not change. The TiO band depths in these spectra are measured from spectral indices as described in Herczeg & Hillenbrand (2014). The TiO indices are also combined into a single TiO index by calculating the number of standard deviations each point is from the median value of the index and then averaging the standard deviation.

The TiO bands trace the cool component. When the hotter component has a larger filling fraction, the TiO bands are veiled by the hotter component and should, therefore, be shallower. As expected, the TiO indices all correlate strongly with the optical brightness. The TiO-7140 band depths are accurately reproduced by varying the two-temperature fit with rotation. Other bands have small offsets (see middle panel of Figure 11), likely caused by either contributions from intermediate temperatures or by small errors in the synthetic spectra. In the spectral type scheme derived by Herczeg & Hillenbrand (2014), the change in TiO-7140 index corresponds to a range from M1.5 (3640 K) to M2.2 (3530 K).

5. Placing LkCa 4 on the HR diagram

5.1. Extinction estimates to LkCa 4

An accurate extinction estimate from the full SED is challenging because of a strong dependence on the distribution of temperatures in the photosphere, which in our work is simplified to a two temperature fit. However, isolating specific colors constrains the possible range of extinctions. Based on the following analysis, an $A_{V}=0.3$ mag is adopted throughout this paper, although some uncertainty is introduced by confusion in the choice of template colors.

The $B-V$ color is stable and not affected by cooler components, and may therefore be compared with the expected colors for the hotter (4100 K) temperature. The average $B-V=1.41$ mag color from Grankin et al. (2008) has a color excess $E(B-V)=0.09$ mag relative to the empirical color expected for a main sequence star and $E(B-V)=0.20$ mag for colors of pre-main sequence stars, based on the temperature-color tables of (Pecaut & Mamajek, 2013). The main sequence colors are probably most appropriate in this analysis because they are freer of spots. The $B-V$ color excess leads to an extinction $A_{V}=0.28\pm 0.15$ mag relative to main sequence colors, or $0.62\pm 0.15$ mag relative to pre-main sequence colors. The listed uncertainties are calculated from extinction estimates relative to colors of 4000 and 4200 K photospheres. The $V-R_{J}$ color excess leads to $A_{V}=0.48$ mag, which should be less affected by gravity (based on $V-I_{C}$ colors in Pecaut et al.) but is likely to be overestimated because the cool component contributes to the $R_{J}$-band emission, as seen in Fig. 10. In the near-IR, the estimated $V=12.61$ mag at the time of the 2MASS observations leads to a spot-weighted template (main sequence) color $J-K_{s}=0.87$ mag, with a range of $\sim 0.05$ mag for acceptable solutions with other temperatures and covering fractions. The observed $J-K=0.93$ mag leads to $A_{V}=0.35$ mag, with a range of $0.29$ mag to account for the range of other acceptable photospheric solutions. The template near-IR color is more robust than optical colors to differences between the main sequence and pre-main sequence temperature scales because the pre-main sequence $J-K$ color is redder at 4100 K and bluer at 2750 K.

The near-IR extinction to LkCa 4, and perhaps to other TTSs, will be overestimated if emission is generated from a multi-temperature spectrum but template colors are based on only the hotter component. This mismatch likely explains why near-IR extinctions are often higher than optical extinctions for young stars (Herczeg & Hillenbrand, 2014), IR colors that are often redder than expected from models (Tottle & Mohanty, 2015), and general challenges in near-IR templates of T Tauri stars (e.g. Espaillat et al., 2010; Fischer et al., 2011). For LkCa 4, the excess $E(J-K)$ color leads to an extinction of $0.84$ mag for a 4100 K photosphere and Pecaut & Mamajek (2013) pre-main sequence colors. This extinction is consistent with the color and SED analysis of Furlan et al. (2006), where a K7 spectral type was adopted, and is higher than extinctions of $A_{V}=0.69$ mag measured from optical photometry (Kenyon & Hartmann, 1995) and $0.35$ mag measured by comparing optical spectra to young spectral templates (Herczeg & Hillenbrand, 2014).

Our results tentatively support the lower extinction obtained at optical wavelengths, although the template colors used to calculate extinction are uncertain and will need revision. Differences between main sequence and pre-main sequence temperature scales and colors in the Pecaut & Mamajek (2013) tables are introduced by a combination of spots, in addition to any color differences caused by gravity differences. Uncertainty in colors may be solved with self-consistent comparisons to young stars with similar spectral types if spot coverage is also similar. Discrepancies in colors, spectral type/effective temperatures, and extinctions for LkCa 4, and perhaps more broadly for TTSs, are dominated by ignoring spots and not by spectral type-effective temperature scales. Fundamental problems in molecular opacities and photospheric models are likely important secondary effects.

5.2. Spot-corrected HR diagram placement

The spectrum of LkCa 4 is well reproduced by a cool component of $2750$ K covering 80% of the surface and a $4100$ K component covering 20% of the surface. Scaling this composite spectrum of LkCa 4 to the observed 2MASS photometry ($\sim 74\%$ spot covering fraction) yields $R=2.3$ $R_{\odot}$ and $\log L/L_{\odot}=-0.20$, adopting $A_{V}=0.3$ mag and a distance of $131\pm 2$ pc888This distance is consistent with the Gaia distance of $130\pm 4$ pc measured to the nearby Taurus members V410 Tau, V773 Tau, and BP Tau (Gaia Collaboration et al., 2016). from the range of parallaxes of nearby Taurus members measured by Torres et al. (2012). At the median 80% spot covering fraction from 2015, the luminosity would be $\log L/L_{\odot}=-0.26$.

The effective temperature $T_{\rm eff}$ of these two components may be calculated as

[TABLE]

For the median fill factor of 80%, the observed effective temperature would be 3180 K, with a range of 3070–3280 for 86% and 74% spot fill factors. If the cool spot is 3000 K, at the limit of the acceptable temperature range, then these values would increase by 150 K. The apparent effective temperature in any single epoch depends on the relative fractions of the hot and cool spot. When accounting for the visible contributions of both hot and cool components, this wide range of possible observed effective temperatures in any single epoch underscores the difficulty in placing young stars on HR diagrams.

Figure 12 and Table 4 show that this effective temperature is much cooler than previous measurements from optical spectra, which are clustered around 4000 K. The luminosity does not change significantly because the broadband SED yields similar luminosities.

When the photosphere is assumed to be a single temperature, the effective temperature depends on the wavelength where the measurement is made. For LkCa 4, the $K$-band temperature is closest to the effective temperature of the entire photosphere. These discrepancies may also explain the large temperature, systematic offsets identified when comparing temperatures from $H$-band APOGEE spectra to those obtained from other (usually optical) methods (Cottaar et al., 2014). Two component photopheres would also explain deriving different spectral types in the optical and infrared for the same source. On the other hand, many GAIA-ESO spectra cover only a narrow wavelength around H$\alpha$ (e.g. Jeffries et al., 2014; Frasca et al., 2015), which would yield an effective temperature of $\sim 4000$ K. Spectral types and effective temperatures measured from TiO depths (e.g. Herczeg & Hillenbrand, 2014) are cooler than those estimated from blue spectra but still overestimate the effective temperature. Temperature uncertainties for young K and M dwarfs are larger than any formal uncertainties in individual measurements.

In principle, the improved characterization of the photospheric emission and radius of LkCa 4 should lead to more accurate estimates of mass and age from pre-main sequence tracks. However, most standard pre-main sequence stellar evolution models do not take into account the effects of starspots, and those that do (Somers & Pinsonneault, 2015a) limit cool spot areal coverage fractions to less than 50%. Comparison of LkCa 4 to unspotted models is inherently flawed, but for the purposes of gauging the effect size of spots on derived parameters, our effective temperature and luminosity yield a mass of $0.15-0.3$ $M_{\odot}$ and an age of $\sim 0.5$ Myr in the Baraffe et al. (2015) evolutionary models999The Baraffe et al. (2015) evolutionary tracks were calculated using BT-Settl photospheric spectra (Allard, 2014), which may introduce minor uncertainties when compared to results from the PHOENIX models used here.. This age is much younger and this mass is much lower than previous estimates made assuming the warm component solely accounted for the effective temperature.

Even though LkCa 4 is only a single data point, the cool spot coverage fraction may be evidence that strong magnetic fields are inhibiting convection in LkCa 4. As a result, the surface is cooler and releases less radiation, slowing the contraction rate. LkCa 4 then appears more luminous and cooler than expected for a star of its genuine age from evolutionary tracks that do not consider magnetic fields. This shift to lower temperature may also help to resolve some of the age discrepancies between intermediate mass and low mass young stars (e.g. Herczeg & Hillenbrand, 2015). The corrected placement of spotted stars on the HR diagram reflects the accurate luminosity, radius, and effective temperature, but evolutionary models are not yet equipped to interpret such heavily spotted stars (Somers & Pinsonneault, 2015a). In an HR diagram, spotted evolutionary model tracks at a given mass and age shift to higher luminosity and lower effective temperature. Such a spot-corrected HR diagram position would alter the age and mass estimates to $>0.2-0.3$ $M_{\odot}$ and $\>0.5$ Myr; in other words, back in the direction of previous naïve estimates. However, the full magnitude of the correction is an outstanding question.

The IGRINS and ESPaDOnS spectra indicate $v\sin i\sim$ 28 km s*-1* and a rotation period of 3.375 days. Combining these numbers gives $R\sin{i}\sim 1.87R_{\odot}$. Our HR diagram analysis gives $R\sim 2.3R_{\odot}$, which suggests $\sin{i}\sim 0.813$, or an inclination of about 35*∘* from edge on. These values show broad consistency between rotational properties, spectral fitting values, and our interpretation of a tilted star revealing a circumpolar region with large areal coverage of starspots.

As seen in Figures 1 and 2, the areal coverage of spots appears to wax and wane secularly through the observing seasons of the last 31 years. LkCa 4 is currently in a relatively faint phase, 0.2 mag fainter in $V$ than it was through much of the late 1980’s and early 1990’s. It is conceivable that LkCa 4 happens to be going through a short-lived episode of relatively high coverage of cool photosphere analogous to— albeit distinct from—solar maximum. The intensity and duration of such episodes would further confound the interpretation of any instantaneous position on the HR diagram.

6. Discussion

In this paper, we have demonstrated that the IGRINS spectra of LkCa 4 are well fit with two temperatures, a hot photosphere of 4100 K and large cool patches of 2700–3000 K that cover 80% of the stellar surface. This solution reproduces the high-resolution $HK$-band spectra obtained with IGRINS, the broadband SED, and rotational modulation of spectral features. Such a high covering fraction for a spot seems extraordinary, but may be common among young stars. In this section, we describe that the need for a large cool spot is robust, despite some simplifying limitations of our spectral inference technique. We also discuss how such a large spot could be interpreted, and suggest that even though LkCa 4 appears extreme, other young sources have similar characteristics.

6.1. An absolute lower limit of starspot filling factor

This simple thought experiment demonstrates that the high starspot filling factor must be approximately correct. A lower limit to the starspot filling factor can be set by the size of a black spot that can induce such large $\Delta V$. Consider the following extremely conservative assumptions:

The starspot emits zero flux. This statement is equivalent to the starspot being at absolute zero. 2. 2.

The stellar disk shows no starspots at maximum brightness. 3. 3.

The starspot transits the center of the stellar disk. 4. 4.

The star shows no feature brighter than the quiet photosphere at the visible surface

For LkCa 4 in 2015 the $\Delta V\sim 0.5$ mag amplitude corresponds to an absolute minimum fill factor $f>0.37$. At least 37% of the stellar disk is covered in starspots at minimum light. In the 2004 season, LkCa 4 exhibited $\Delta V\sim 0.8$, corresponding to a minimum fill factor of $f=0.52$. If any of assumptions 1-3 is relaxed, the minimum fill factor increases. If assumption 4 is relaxed, the minimum fill factor decreases. The large amplitude photometric variations in LkCa 4 are therefore evidence for either a hot spot in excess of the ambient photosphere or a large coverage fraction of cool starspots.

6.2. Interpretation of the large filling factor of cool starspots

We have shown that LkCa 4 exposes between 74% and 86% cool photosphere towards our viewing direction in winter 2015-2016, the most recent season for which there is data. In previous seasons, like fall 2004 (c.f. Figures 2 and 9), the filling factor of cool spots could have exceeded 90%. This exceptionally large fraction of surface area covered by spots defies the conventional wisdom that starspots are the minority constituent of the stellar surface. The evidence for this derived fill factor has been shown in the SED, high-resolution spectral inference, direct spectral detection, and TiO variability modeling. Perhaps the strongest evidence for the large fill factor of cool photosphere is the large amplitude of $V-$band modulation, which cannot be explained with small spots, even if a spot temperature of absolute zero is assigned to the cool regions.

One conceivable interpretation of the 74%–86% coverage fraction of cool photosphere is a stellar axis inclined towards our observing direction, with a large circum-polar spot covering much of the observable hemisphere, as described by Donati et al. (2014). The ambient hot photosphere would then be close to the equator, which would be limb-darkened due to projection effects. Such limb darkening would cause a non-linear transformation between our observed fill factor $f_{\mathrm{cool}}$, and the surface areal coverage fraction of spots. It is conceivable, and probably likely, that the value we report for $f_{\mathrm{cool}}$ is an overestimate of the surface area actually covered by spots. In this circum-polar spot scenario, local fingers of cool spots could reach down to near the equator, entering into and out of the field of view and giving rise to the large amplitude photometric variation seen on LkCa 4. This scenario is consistent with the inclination derived from relating $v\sin i$, rotation period, and estimates for the stellar radius, as previously noted. Such estimates are limited by estimates of the radius since uncertainties in effective temperature and luminosity cause a wide range of radius estimates.

The starspot geometry is also consistent with morphologies observed in interferometric maps of $\zeta$ Andromedae: large areal coverage of permanent polar spots, and smaller, transient networks of varying darkness which thread across the whole surface (Roettenbacher et al., 2016).

6.3. Is LkCa 4 an extreme source?

The large amplitude of cyclical optical variability of LkCa 4 exceeds all other counterparts from the long-term monitoring of WTTSs by Grankin et al. (2008). From a statistical perspective, it might not be a surprise that the source with the largest photometric variability also possesses among the largest coverage fraction of cool spots. When compared with spot covering fractions of 15–60% for evolved giants and subgiants (Chugainov, 1976; Berdyugina, 2005), the 80% spot coverage of LkCa 4 seems high. However, photometric and spectroscopic evidence demonstrates that many other young low-mass stars also have large spot covering fractions.

Several WTTS in the Taurus-Auriga star-forming region have large amplitudes in periodic light variations, achieving $0.4-0.8$ mag in the $V-$band: LkCa 4, V410 Tau (Herbst, 1989), LkCa 7, V827 Tau, V830 Tau, and V819 Tau (Grankin et al., 2008). In a sample that is not biased to variable stars, WTTSs observed in NGC 2264 show photometric variations similar to–but of lower average amplitude than– LkCa 4 (Cody et al., 2014; Venuti et al., 2015). Modeling of polychromatic photometric monitoring of the locus of WTTSs indicates a lower limit to the true coverage fraction $f_{\mathrm{eq}}\sim 10-30\%$ and a difference $T_{\textrm{hot}}-T_{\textrm{cool}}\sim 500$ K.

Large amplitudes of the light variation indicate the existence of very extended spotted regions on the stellar surface. To estimate the total area and mean temperature of the spots in the visible stellar hemisphere Grankin used a simple non-parametric model for analysis of V410 Tau, V819 Tau, V827 Tau, and V830 Tau (Grankin, 1998, 1999). For the stars with the highest light curve amplitudes, this model shows that the spotted regions cover from 17 to 73% of the visible stellar hemisphere and that the mean temperature of the cool regions is $500-1400$ K lower than the ambient photosphere. An analysis of the calculated spot parameters of V410 Tau revealed a correlation between the amplitude of periodicity and the spot distribution nonuniformity, which is defined as the difference of spot coverage in the visible hemispheres of the star at maximum and minimum light. As the amplitude increases from 0.39 to 0.63 mag, the degree of nonuniformity in the spot distribution increases from 21 to 37%. An alternative analysis of TiO bands indicate a spot coverage of $\sim 41$% for V410 Tau (Petrov et al., 1994). Besides the high amplitude of variability, these objects show the phenomenon of long-term stability of a brightness minimum in the interval from 5 to 19 years.

The Pleiades provides an older ($\sim 125$ Myr, Stauffer et al., 1998) sample of young stars, many of which have starspots that cover large fractions of the stellar surface. For example, anomalous colors of K dwarfs in the Pleiades can be explained by $\sim$50% fill factors of spots (Stauffer et al., 2003), while LAMOST spectra of the TiO band indicates that $\sim 30$% spot coverages may be common for 3500–5000 K stars (Fang et al., 2016). Recent photometric monitoring of the Pleiades indicates that the anomalous colors scale with rotation rate, strengthing the case for magnetic-field induced starspots as the cause of the anomaly (Covey et al., 2016). However, the Covey et al. (2016) Pleiades rotation rates do not scale with the $\Delta V$ photometric amplitude at the precision of the PTF data, suggesting that while the overall starspot coverage increases with rotation rate, the longitudinally asymmetric component does not. LkCa 4 represents an example of such an architecture, with more than 50% of the visible stellar surface covered in spots at all observational phases. The longitudinally symmetric component evades photometric modulation detection but contributes to otherwise anomalous color offsets.

Spots are also likely important in estimating the properties of generic low-mass young stars. As an example, the well-studied CTTS TW Hya has a $\sim 4100$ K photosphere measured from high-resolution optical spectra (e.g. Yang et al., 2005). However, the near-IR spectrum is consistent with $\sim$ M2, when compared with main sequence dwarf stars (Vacca & Sandell, 2011). Debes et al. (2013) reconciled this discrepancy by invoking spots, although the spectral contribution from accretion and disk continuum emission was uncertain. McClure et al. (2013) found that the TW Hya spectrum was well reproduced with near-IR spectra of stars with similar ($\sim$K7) optical spectral type, as long as the comparison spectra were young. Herczeg & Hillenbrand (2014) found an intermediate spectral type for TW Hya and other TTSs by focusing on TiO bands at redder wavelengths than previous spectra, and also by using young stars as photospheric templates. Similar discrepancies have been found between optical and near-IR spectral types (e.g. Bouvier & Appenzeller, 1992; Bary & Petersen, 2014; Cottaar et al., 2014) and in comparisons of spectral types measurement based on TiO bands and those based on high-resolution blue spectra (see discussion in Herczeg & Hillenbrand, 2014).

The mismatches between near-IR spectra of dwarf stars and young stars (along with some of the color discrepancies, see Section 5.1) is also likely explained by the prevalence and importance of cool spots in low mass stars. For classical T Tauri stars, the veiling at $\sim 1$ $\mu$m is usually larger than can be explained by either most accretion models101010Ingleby et al. (2013) invoked multiple accretion flows with a range of densities to explain this red excess. or warm disk emission (e.g. Basri & Batalha, 1990; Cieza et al., 2005; Fischer et al., 2011; McClure et al., 2013). Similarly, Cottaar et al. (2014) found that most cool stars in the Pleiades are veiled in the H-band, with no possibility for any disk or accretion emission. Differences in spot temperatures and coverages likely contribute to and may explain entirely this excess veiling.

6.4. The limitation of a two-temperature model

We assume only two temperature components are present in the photosphere. In reality, starspots are probably described by a range of spot temperatures analogous to umbra, penumbra, and plages on the solar surface. A model possessing three unique temperatures and two unique fill factors would better approximate the minor constituents of the emergent spectrum than our two temperature model can, but additional components will always appear to fit data better while leading to possible overfitting. The fit-quality of the existing two-temperature fits seems strong through much of $H-$band and many $K-$band orders, so the contribution of minor photospheric constituents is relatively low. Interestingly, the $K-$band has several spectral orders that converge towards hybrid temperatures in the vicinity of $\sim 3500$ K. Most of these orders are flagged as erroneous due to strong metal line spectral outliers or uncorrected calibration artifacts. A few orders just show mediocre fits where no permutation of hot and cool two-temperature models could describe all of the variance in the spectrum. These $K-$band orders could be examined for a penumbra-like temperature component. We caution that systematic errors in the models would ultimately limit the interpretability of models with more than two temperature components. Space-based high-precision planetary transit spectrophotometry could offer an avenue to measure the morphologies and temperature distributions of starspots for the subset of transiting exoplanet host stars or eclipsing binaries.

Analyses of optical emission also suggest that a two-temperature model may be an oversimplification. Figure 13 shows that the $V-R$ colors from Grankin et al. (2008) are better reproduced if the cooler component is $\sim 3050$ K, thereby contributing more flux at optical wavelengths. The ZDI imaging of Donati et al. (2014) also indicates a complex structure of spots that are hotter than what is interpreted here as the ambient photosphere, while also reproducing the optical photometry. Despite the need for more complicated emission structures to simultaneously reproduce all observations, the broadband SED, spectral shape, and rotational modulation of spectral indices are reasonably well fit with the two components presented here. The high and low-resolution near-IR spectra, including the H2O band absorption between the $J$ and $H$ bands (see Figure 8), require a cool (2700–3000 K) component while the global photometry requires a large covering fraction for this cool component.

6.5. Limitations of assumptions in the spectral inference methodology

We assumed that the Phoenix synthetic spectra are a good representation of the observed IGRINS spectrum of LkCa 4, a young star that possesses a strong magnetic field (Donati et al., 2014). Cool sunspots are associated with locally heightened magnetic fields on the sun, so synthetic spectral models employing negligible magnetic fields will fail to fit magnetic sensitive features. Magnetic fields will have a few effects on the emergent spectrum- line broadening attributable to the Zeeman effect, and pressure broadening attributable to seeing deeper into the photosphere in the line wings than in the line center. The former effect will preferentially impact lines with high magnetic sensitivity (e.g. Johns-Krull et al., 1999; Deen, 2013). The latter effect could masquerade as an apparent shift in surface gravity, as seen in mixture models of spotted stars with empirical templates composed of dwarfs and giants (O’Neal et al., 1996). Zeeman amplification can also be present (Basri et al., 1992).

Our tacit assumptions have been that these line-broadening effects are 1) relatively small compared to the gross appearance and disappearance of spectral lines at disparate temperatures, 2) many orders and lines are not affected at all, and 3) the Gaussian process noise model provides some resilience against unknown residual spectrum correlations induced by such under-fitting. We see tentative evidence for physics outside of the pre-computed model grids in the poor fits to the observed spectra in $K-$band. Strong metal lines demonstrate large departures in their best fit $[{\rm Fe}/{\rm H}]$, as the naïve fitting strategy futilely attempts to minimize the residuals from Zeeman-broadened lines. No pre-computed model grids including a range of surface magnetic fields are available at the wavelength range and spectral sampling needed to approach these questions within this spectral inference framework.

Uncertainties in the Phoenix synthetic spectra continuum opacity flow down to uncertainties in the absolute value of the fill factor. Erroneous bulk opacities could be responsible for some of the scatter of fill factor distributions seen across different spectral orders, as in Figures 5 and 4. The large number of spectral orders over which the fill factors are averaged offers resilience against these anomalies.

6.6. Evidence for intra-spectrum RV jitter

Starspots cause significant jitter in radial velocity (RV) surveys for young planets (e.g. Donati et al., 2014; Robertson et al., 2014). The two-temperature model assumed that the hot and cool components shared the same $v\sin i$ and $v_{z}$. This assumption is accurate for either homogeneously distributed spots, azimuthally symmetric bands of spots, or polar spots. LkCa 4 has clear optical line-profile variations (Nguyen et al., 2012; Donati et al., 2014). We should be able to detect an anticorrelation between the line center positions of spectral lines arising from the ambient photosphere and lines arising from cool starspots, since there is a zero-sum competition for solid-angle on the recessional and advancing stellar limbs. We found tentative, albeit weak, evidence for such anticorrelation in $v_{z}$ across spectral orders dominated by cool photosphere and those dominated by hot photosphere. For example, the cool-photosphere dominated spectral orders $100,102,104$ ($1.7$–$1.8$ $\mu$m) show the largest median $v_{z}\sim 14.5-17.0$ km/s, while the other orders that possess a mix of primarily hot- and some cool- photosphere features demonstrate $v_{z}\sim 9-13$ km/s. An independent analysis of the RV of the IGRINS spectrum yields a barycentric RV of $14.6\pm 0.2$ km s*-1*. The direct emission of starspots could constitute a heretofore unaccounted source of correlations in conventional RV signal processing.

Radial velocity (RV) planet searches in the infrared–where the modulation of RV signals is lessened (e.g. Prato et al., 2008; Crockett et al., 2012)– could potentially reduce RV jitter by characterizing and constraining the effect of direct emission from starspots, in addition to the flux deficits imbued in optical line profiles. Direct spectral emission from cool starspots is likely much weaker in sources with less starspot coverage than LkCa 4. As demands on RV precision increase, the characterization and modeling of direct emission from starspots should in principle decrease RV jitter attributable to starspots in the near-IR spectra of any spotted star.

7. Conclusions

This paper establishes that the optical and near-IR spectra of the weak-lined T Tauri star LkCa 4 are produced by multiple temperature components. When the temperature distribution is simplified to two temperatures, the photospheric model that best fits the spectrum includes a hot (4100 K) photospheric region that covers 20% of the stellar surface and a cool (2700–3000 K) photospheric that covers 80% of the stellar surface. The previously-established $V$-band brightness variability of $0.8$ mag is caused by the rotational modulation of the covering fraction of the hot photospheric component. This spectral modeling leads to an improved effective temperature measurement of LkCa 4. For LkCa 4 and other young low-mass stars, the standard approach of inferring masses and ages of young stars is inappropriate and will lead to large uncertainties, both because effective temperatures are systematically too hot and because most models of pre-main sequence evolution neglect spots.

New results from Fang et al. (2016) show that large starspot coverage fractions of $50\%$ are common even among stars in the Pleiades, which at $\sim 125$ Myr are much older than the 1–10 Myr old stars in the Taurus Molecular Cloud. LkCa 4 breaks the apparent ceiling at 50% coverage fraction of cool starspots, though it is not clear how to classify stars for which the minority of the surface is purportedly the ambient photosphere. A high fraction of the stellar surface may be covered in spots without inducing any photometric modulation, if the spots are located in circumpolar regions that always face towards our line-of-sight, in longitudinally symmetric bands, or in many isotropic small spots.

Coupled with the results of Fang et al. (2016) and Covey et al. (2016), the results shown here suggest that the contribution of starspots to pre-main sequence stellar evolution have been systematically underestimated. Some cool stars in the Pleiades have a large excess in Li abundance, some have inflated radii, some are fast rotators, and some have spots (e.g. Somers & Pinsonneault, 2015b). While the source overlap for these properties is unclear, these anomalies may all be explained if the pre-main sequence evolution of some stars is strongly affected by spots (Somers & Pinsonneault, 2015a). LkCa 4 and other young stars with similar properties may help to connect the properties of anomalous Pleiades stars with earlier stages of stellar evolution.

The spectral inference technique described in this paper lays the foundation for future studies that seek accurate constraints on multiple photospheric components. The near-IR has been shown to be the most sensitive region to search for cool emission from starspots. Starspots and disk or accretion veiling can be modeled in this framework, and combining the two will enable studies of the large number of classical T-Tauri stars. Applying this technique to large ensembles of T-Tauri stars across a wide range of mass and age will offer an observational picture of starspot evolution that can revise the current understanding of pre-main sequence stellar evolution.

Appendix A Further methodological details

When the synthetic stellar models used in Starfish are in absolute flux units, our mixture model is complete as stated in Equation 1. However C15 employ standardized $f_{\lambda}$ when constructing the forward model $\mathsf{M}$:

[TABLE]

The reason for standardizing the flux is a matter of practicality: the number of PCA eigenspectra components in the spectral emulator scales steeply with the pixel-to-pixel variance of $f_{\lambda}(\{{\bm{\theta}}_{\ast}\}^{\textrm{grid}})$, increasing the computational cost of spectral emulation. The choice to standardize fluxes makes no difference for modeling a single photospheric component. But for two-component photosphere models, the relative flux of the two model spectra needs to be accounted for to get an accurate estimate of the areal coverage fraction of the cool spots, $f_{\mathrm{cool}}\equiv\Omega_{\mathrm{cool}}/(\Omega_{\mathrm{hot}}+\Omega_{\mathrm{cool}})$. So we scale the mixture model in the following way:

[TABLE]

where the prime symbol in the final line indicates a re-standardized mixture model flux, where the relative fluxes of the model components are now correctly scaled.

We are then tasked with computing an estimator, $\hat{f}({\bm{\theta}}_{\ast})$, for estimating the scale factor $q$ in-between model gridpoints. One robust approach would be to follow C15 by training a Gaussian process on the $f(\{{\bm{\theta}}_{\ast}\}^{\textrm{grid}})$. Instead, we simply linearly interpolated between the model grid-points. Interpolation can cause pile-up near model grid-points, as noted in Cottaar et al. (2014), which motivated the spectral emulation procedure in C15. We assume the interpolation of $\hat{f}$ is smooth enough that we will not see such pileups and if even we did see pileups, they would mostly be discernable as kinks in the distribution of samples in the starspot coverage fraction $f_{\mathrm{cool}}$. One drawback of our estimator method compared to the Gaussian Process regression method is that we do not propagate the uncertainty associated with the absolute flux ratio interpolation into our estimate of $f_{\mathrm{cool}}$. We assume this uncertainty is relatively small and can be ignored.

The mixture model is a linear operation with all the same stellar extrinsic parameters, so we can re-use all the same post-processed eigenspectra $\widetilde{\mathbf{\Xi}}$, mean spectrum $\widetilde{\xi}_{\mu}$, and variance spectrum $\widetilde{\xi}_{\sigma}$, with the tildes representing all post processing except the $\Omega$ scaling. We calculate two sets of eigenspectra weights $\mathbf{w}_{\in(\mathrm{hot},\mathrm{cool})}$, and their associated mean and covariances following the Appendix of C15, and yielding:

[TABLE]

[TABLE]

Appendix B Further examples of full spectrum fitting

Figures 14 and 15 show spectra from 42 IGRINS orders.

Appendix C Full table of IGRINS best fits

Table 3 lists all the IGRINS spectral order results.

Appendix D Previous work

Table 4 lists measurements of LkCa 4 from previous studies, compared to this work.

We thank the anonymous referee, Chris Kochanek, and Daniel T. Jaffe for valuable comments on the draft. MG-S and GJH also thank Steve Saar, Lisa Prato, and Lynne Hillenbrand for valuable discussions about spots. Portions of the data, Python code, Jupyter notebooks, Starfish configuration files, MCMC samples, and revision history of this manuscript are available on the project’s GitHub repository at https://github.com/BrownDwarf/welter. MG-S and GJH are supported by general grant 11473005 awarded by the National Science Foundation of China. The ESPaDOnS observations are supported by the contribution to the MaTYSSE Large Project on CFHT obtained through the Telescope Access Program (TAP), which has been funded by the “the Strategic Priority Research Program—The Emergence of Cosmological Structures” of the Chinese Academy of Sciences (Grant No.11 XDB09000000) and the Special Fund for Astronomy from the Ministry of Finance. TW-SH is supported by the DOE Computational Science Graduate Fellowship, grant number DE-FG02-97ER25308. We thank LCOGT and its staff for their continued support of ASAS-SN. ASAS-SN is supported by NSF grant AST-1515927. Development of ASAS-SN has been supported by NSF grant AST-0908816, the Center for Cosmology and AstroParticle Physics at the Ohio State University, the Mt. Cuba Astronomical Foundation, and by George Skestos. BJS is supported by NASA through Hubble Fellowship grant HST-HF-51348.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. This work used the Immersion Grating Infrared Spectrograph (IGRINS) that was developed under a collaboration between the University of Texas at Austin and the Korea Astronomy and Space Science Institute (KASI) with the financial support of the US National Science Foundation under grant ASTR1229522, of the University of Texas at Austin, and of the Korean GMT Project of KASI. This research has made use of NASA’s Astrophysics Data System. The research is based on data from the OMC Archive at CAB (INTA-CSIC), pre-processed by ISDC. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. Facilities: Smith (IGRINS), AAVSO, CFHT (ESPaDOnS), INTEGRAL (OMC), ASAS, CrAO:1.25m, ARC (TripleSpec), Hale (DBSP), Gaia *Software: * pandas (McKinney, 2010), emcee (Foreman-Mackey et al., 2013), matplotlib (Hunter, 2007), numpy (Walt et al., 2011), scipy (Jones et al., 2001–), ipython (Pérez & Granger, 2007), gatspy (Vanderplas, 2015), starfish (Czekala et al., 2015), seaborn (Waskom et al., 2014)