Possible Bright Starspots on TRAPPIST-1

The M8V star TRAPPIST-1 hosts seven roughly Earth-sized planets and is a promising target for exoplanet characterization. Kepler/K2 Campaign 12 observations of TRAPPIST-1 in the optical show an apparent rotational modulation with a 3.3 day period, though that rotational signal is not readily detected in the Spitzer light curve at 4.5 $\mu$m. If the rotational modulation is due to starspots, persistent dark spots can be excluded from the lack of photometric variability in the Spitzer light curve. We construct a photometric model for rotational modulation due to photospheric bright spots on TRAPPIST-1 which is consistent with both the Kepler and Spitzer light curves. The maximum-likelihood model with three spots has typical spot sizes of $R_\mathrm{spot}/R_\star \approx 0.004$ at temperature $T_\mathrm{spot} \gtrsim 5300 \pm 200$ K. We also find that large flares are observed more often when the brightest spot is facing the observer, suggesting a correlation between the position of the bright spots and flare events. In addition, these flares may occur preferentially when the spots are increasing in brightness, which suggests that the 3.3 d periodicity may not be a rotational signal, but rather a characteristic timescale of active regions.


INTRODUCTION
TRAPPIST-1 is an M8V star and host to at least seven Earth-sized planets (Gillon et al. 2016(Gillon et al. , 2017. It was observed with Spitzer at the IRAC-2 4.5 µm band for 20 days (Gillon et al. 2017;Delrez et al. 2018), and 76 days later it was observed by NASA's Kepler spacecraft during K2 Campaign 12 for another 79 days Luger bmmorris@uw.edu * Guggenheim Fellow † DIRAC Fellow NSF Astronomy and Astrophysics Postdoctoral Fellow (2017). The K2 observations of TRAPPIST-1 show quasi-periodic flux modulations with a period of 3.3 days Delrez et al. 2018), which has been interpreted as a signal of stellar rotation as starspots rotate into and out of view (Vida et al. 2017). Modeling these variations faces the vexing computational problem of fitting stellar rotational flux modulation with a forward model (Aigrain et al. 2015;Davenport et al. 2015). One challenge in modeling rotational modulation is that often a light curve can be fit equally well with a few bright spots or with a few dark spots. This degeneracy has been referred to as the "zebra effect" (Pettersen et al. 1992): is the star bright with a few dark spots, or dark with a few bright spots? Although star spots are usually dark due to magnetic field pressure balancing thermal pressure, one star with a mostly dark surface and a few bright regions is T Tauri star LkCa 4 (Gully-Santiago et al. 2017). Near-infrared spectroscopy of LkCa 4 shows a heterogeneous atmosphere -80% of the star has T eff ∼ 2800 K, and the other 20% has T eff ∼ 4000 K. Modern Zeeman Doppler Imaging is most sensitive to large scale features in stellar magnetic fields, so it is generally unclear what small-scale magnetic structures may exist on fully-convective stars (see e.g. Morin et al. 2013).
Another difficulty in interpreting the TRAPPIST-1 brightness variations is whether the periodicity is due to rotation at all. Reiners & Basri (2010) measure a rotational velocity of 6±2 km sec −1 , while the 3.3-day period would predict a maximum rotational velocity of 1.8 km sec −1 . The spots on TRAPPIST-1 evolve with time, and so it will require long-term monitoring of the 3.3-day periodicity to determine whether it is time-steady, which would argue for a rotational origin, and would call into question the measurement of the rotational velocity.
Another means of studying spots is via occultation during planetary transits, which has been studied in detail on HAT-P-11 (Morris et al. 2017). To date, no significant spot occultations have been observed in either the Kepler or Spitzer transit light curves of TRAPPIST-1 (Delrez et al. 2018) -this implies that if spots are responsible for the rotational modulation, they must either reside outside of the transit chords of the planets, or be smaller than the size of the planets, which would reduce the probability and amplitude of a spot occultation.
In this paper we present a model for the K2 and Spitzer variations based upon a spot variability model, and examine its relation to the stellar flares that are frequently seen in both datasets. In section 2 we discuss the spot model for stellar variability. In section 3 we compare this model with extant K2 and Spitzer data. In section 4 we discuss a possible relation between bright flares and the bright spots. We end with discussion and conclusions.

Monochromatic variability
We compute spot variability at a given wavelength with a simplified model assuming small spots with radii smaller than 10% of the stellar radius. We first integrate the total stellar flux of the unspotted, limb-darkened star, where r is the radial coordinate from the center of the star, D is the distance to the star, I(r) is the specific intensity at radius r. We model the specific intensity with a quadratic limb-darkening law, using the limb-darkening parameters of Luger et al. (2017): (u1, u2) = (1.0, −0.04). We describe each spot with an ellipse with centroid r i = (x i , y i ) at radial coordinate r i = |r i |. We then compute the flux contribution from each spot by computing the approximate spot area and spot contrast based on the ratio of Phoenix model atmospheres integrated over a given bandpass, which we assume is independent of inclination angle of observation. From the observer's perspective, each spot has semimajor axis R spot along the azimuthal direction, and semiminor axis R spot 1 − (r i /R ) 2 in the radial direction, due to foreshortening. Since these spots are small compared to the stellar radius, R spot /R << 1, we adopt a constant limb-darkened contrast for the entire spot, c ld = (c−1)I(r i )/I(0), where c is the monochromatic contrast in the spot relative to the adjacent photosphere flux; c = 1 for an unspotted region. The integrated spot flux of each spot is then given by where the spotted flux of the star is As the star rotates, the flux of the spots varies due to foreshortening, limb-darkening, and disappearance behind the edge of the limb. Spots straddling the limb are ignored until they rotate onto the observer-facing hemisphere. This approximation is valid for spots that are small compared to the stellar radius, and small compared to the scale of limb-darkening variation across the stellar disk. We then compute the flux observed as the star rotates, by rotating the positions of the spots and recalculating the fluxes.
We find the posterior uncertainties on spot positions (latitude and longitude), radii and contrast in the Kepler bandpass with Markov Chain Monte Carlo (Foreman-Mackey et al. 2013).

Spot contrasts in Kepler and Spitzer
One might expect that the periodicity observed with Kepler at 3.3 days would be present in Spitzer observations at 4.5 µm, but after removing flares and transits from the Spitzer light curve, little signal of rotation is present Delrez et al. (2018), see Figure 1. We compute the autocorrelation function and Lomb-Scargle periodogram for the Spitzer light curve and find one periodic component near 0.5 days -similar to the timescale of super-granulation on the Sun (Aigrain et al. 2003) -in addition to periodicity near 4 days.
To measure the significance of the peaks in the Lomb-Scargle periodogram, we generate 100 simulated light curves -each light curve was a Gaussian process sample drawn from the Upper: Spitzer observations of TRAPPIST-1 after transits and flares have been removed (black circles) (Delrez et al. 2018), and a fit with a Gaussian process assuming a simpleharmonic oscillator kernel (blue curve). Lower: The autocorrelation function and Lomb-Scargle periodogram for the Spitzer observations show a characteristic peak at periods near 0.5 days.
maximum-likelihood kernel fit to the Spitzer observations -and found that 33% of the random light curves produced LS peaks with at least as much power as the 4 d peak, suggesting that all signals in the LS periodogram are insignificant.
We choose a small section of the Kepler/K2 EVEREST light curve to fit with the spot model, see Figure 2 (Luger 2017). We begin by normalizing the light curve by a quadratic, and median-filter the fluxes over a five-cadence kernel to remove most flares. We select a portion of the light curve (black points) that has a repeated pattern across more than one stellar rotation. There are very few repeated flux patterns in the Kepler light curve, indicating that the spot evolution timescale is likely shorter than the rotation period. We construct a spot model with three spots and allow their positions, radii and contrast to vary, with the contrast bounded on 0 < c < 1. The blue curves in Figure 2 are models drawn from the posterior samples, indicating that the three spots have radii near R spot /R ∼ 0.12 and contrast c = 0.6.
The prior applied to ensure that starspots are dark (c < 1) is informed by our studies of the Sun and sunlike stars, which have dark starspots (see e.g. Solanki 2003; Morris et al. 2017). However, bright spots have been observed on nearby brown dwarf the Luhman 16A via Doppler imaging (Crossfield et al. 2014), and on the T Tauri star LkCa 4 (Gully-Santiago et al. 2017).
To investigate the theoretical effects of bright spots or dark spots, we compute the contrast of a spot with spectrum F spot on a star with spectrum F photosphere , using spectra drawn from the PHOENIX+BT Settl model (Husser et al. 2013). We find that where λ is wavelength, and T is the filter transmission curve for Kepler or Spitzer. We compute c for spots on a star with T eff = 2500 K in Figure 4. For Kepler contrasts 0 < c < 1, the Spitzer contrast is 0.5 < c < 1 -in other words, the spot contrast of dark spots in the Kepler band is similar to the spot contrast of dark spots in the Spitzer band. This result is in contradiction with the observed spot modulation in the Kepler bandpass and the undetected modulation in the Spitzer light curve.
Alternatively, the contrast of bright spots in the Spitzer 4.5 µm band increases very slowly as the spot contrast in the Kepler bandpass exceeds c > 1 (Figure 4). A spot 100x brighter than the photosphere in the Kepler band would only be 3x brighter in the Spitzer band. This slow increase in the contrast of bright spots at   long wavelengths could produce large flux modulations in the Kepler bandpass and very small modulations in the Spitzer bandpass for very small spots.

SIMULTANEOUS FITS TO THE Kepler
AND Spitzer LIGHT CURVES We fit the spot model to the Kepler light curve and a segment of the Spitzer observations of the same duration, using the empirical relation from Figure 4 to convert spot contrasts in the Kepler bandpass to contrasts at 4.5 µm. We fit both light curves simultaneously, even though the observations were not simultaneous, so that the Spitzer light curve will constrain the spot radii and contrasts to be consistent with the typical variability in the infrared. The Kepler observations ended 76 days before the Spitzer observations started. We assume that the variability of the star did not change significantly over those 76 days because: (1) the variability does not seem to change significantly throughout the duration of either the Kepler or Spitzer light curves, and (2)  We must first choose the minimum number of spots required to reproduce the K2 light curve. We fit the K2 light curve for n = 1, 2, ..., 6 spots, allowing the spot positions, radii, and contrast to vary, and evaluate the reduced χ 2 and total spotted area of each fit. We find that the reduced χ 2 plateaus near a minimum after 3 or more spots have been added. If more than three spots are modelled, the radii of the spots are decreased in order to keep the total spotted area approximately constant. We therefore fix our spot number to 3 spots since that minimizes the χ 2 and the number of free parameters.
If there are many small spots distributed isotropically on the stellar surface, we would not detect them through rotational modulation. As a result, the rotational modulation fits are only sensitive to the longitudinal asymmetries in the spot distribution. Consequently, the spotted areas inferred from this model should be considered lower limits on the spotted area on the star.
We impose a prior to ensure that the addition of bright, hot spots does not change the color of TRAPPIST-1 significantly. We measure the color of the star from the optical spectrum V − I c = 4.7 (Burgasser et al. 2015). At each step in our Markov chains, we add to the prior a penalty for significant color deviations from the observed color, The maximum-likelihood bright spot model fit to the Kepler and Spitzer light curves 1 is shown in Figure 5. The relative flux is normalized to its minimum, which assumes that the minimum flux within this segment of the light curve represents the unspotted flux. The variability in the Kepler band is reproduced by the spot model, and the corresponding variability in the Spitzer band is comparable to the observational uncertainties.
One spot has radius R spot /R = 0.02 ± 0.002 and the other two have R spot /R = 0.013 ± 0.002. We measure the Kepler contrast c k = 230 ± 40, which corresponds to Spitzer 4.5 µm contrast of c s = 3.7 ± 0.1, and a spot temperature of T spot 5300 ± 200. The uncertainty in the minimum spot temperature is likely underestimated, since the Markov chains prefer a narrow range of spot contrasts to fit the Kepler light curve, but the Spitzer light curve weakly constrains the lower limit on the spot temperature. The total bright spot area coverage is 16 microhemispheres (one hemisphere is half the surface area of the star) -which is small compared to the typical dark spot area coverage on the Sun (Morris et al. 2017).
We repeat this analysis using a different segment of the K2 light curve to ensure that the results are reproducible at different times throughout the rotation of TRAPPIST-1. We choose the fluxes over two rotations spanning 2457762 < BJD < 2457769 and find spot sizes and contrasts consistent with the results from the other segment of the K2 light curve (both regions are labeled on Figure 7).
The autocorrelation function of the Spitzer observations plus the maximum-likelihood spot model serves as a sanity-check on the amplitude of the signal introduced by the inferred bright spots. If injecting the maximum-likelihood spot model into the Spitzer light curve introduces significant periodicity, we should see an uptick in the autocorrelation function at P rot = 3.3 d. Shown in Figure 6, the autocorrelation function is relatively unchanged by the injection of the spot modulation from the maximum-likelihood spot model, and which suggests that this spot model is plausibly consistent with the Spitzer observations.

CORRELATION BETWEEN BRIGHT SPOTS AND FLARES
The flares of active M4 dwarf GJ 1243 have been studied in detail Davenport et al. 2014). The authors searched for a correlation between flare occurrence and starspot phase by comparing the quiescent flux of the star just before a flare to the mean flux. If flares occur near starspots, the quiescent flux of the star just before the flare should be less than the mean. No such correlation emerged, which the authors suggest indicates that the positions of flares and spots are uncorrelated. They also searched for correlation with rotational phase, which could be connected to long-lived polar spots, and found no correlation.
We carried out a similar analysis to investigate if the flares are correlated with starspot phase for TRAPPIST-1, as was briefly noted in Vida et al. (2017). We manually identified the nine largest flares in the K2 EVEREST light curve, and masked them out -see Figure 7. We removed a fifth order polynomial to remove systematic trends without removing stellar vari- ability. To measure the flux of the star at the time of the flares, we fit the rotational modulation with a Gaussian process using a simple harmonic oscillator kernel, and extrapolate the model to the times of flares.
We find that the stellar flux at the time of flares is greater than the typical flux, see In other words, the star is typically brighter just before a flare event than the mean flux. This might suggest that the bright active regions are spatially correlated with the flares on TRAPPIST-1.
The flares also tend to occur when the change in brightness is most positive (see Figure 7). If we interpret the 3.3 day periodicity as rotational modulation, then the preference for flares to occur on the leading edge of the brightening events would indicate that flares are most likely to occur at a particular stellar longitude, in contradiction with the flare analysis of Vida et al. (2017). If we instead interpret the 3.3 day periodicity as the characteristic lifetimes of bright active regions on the star, then the 1% brightening events may be localized brightening associated in time and space with the flaring activity. Until a robust measurement of the stellar rotation period has been made, we cannot rule out the possibility that the 3.3 d period is not the stellar rotation period.  (Foreman-Mackey et al. 2017). We find that flares occur preferentially at times of high flux (see Figure 8). The region labeled "a" is the segment of the light curve in Figures 2 and 5, the region labeled "b" is the validation segment where we repeated the bright spot analysis. Lower: It also appears that the flares occur on the leading edge of the brightening events. We show the numerical derivative of the Gaussian process model in the middle panel, and find that the flares tend to occur when the flux is increasing most rapidly (see also Figure 9).
Bright spots on the Sun referred to as faculae arise due to the magnetic fields and viewing geometry of convective granules (Spruit 1976,  1977). Solar faculae have sizes comparable to convective cells (Keller et al. 2004), and intensity contrasts relative to the mean photosphere that are small (c ∼ 1.5) compared to the spot contrast we observe on TRAPPIST-1 in the Kepler band. The characteristic granule size on late-M dwarfs is likely to be ∼ 80 km across (Ludwig et al. 2002), while the spots on TRAPPIST-1 are ∼ 600 km in diameter. We suggest that the spots on TRAPPIST-1 should not be considered faculae. Few observations have revealed small scale active regions on fully-convective stars. The global-scale magnetic fields of late-M dwarfs have been studied in detail using Zeeman-Doppler Imaging (Donati et al. 2003;Morin et al. 2010Morin et al. , 2011Morin et al. , 2013. One exception is the fully convective star V374 Peg which shows 2% dark spot coverage (Morin et al. 2008).
We can compare our results with Rackham et al. (2017) and Zhang et al. (2018), who computed the spot and faculae covering fractions for TRAPPIST-1. Both studies find that the stellar spectrum is best described by a heterogeneous surface, with persistent spectral components at  Figure 9. K2 light curve (black circles) and the maximum-likelihood Gaussian process models (gray curves), folded at the times of flares. The red curve is the mean of the Gaussian process models. It appears that the flares occur preferentially during the rise in brightness (see also Figure 7). different temperatures. Our constraints, based upon two-band photometry, also imply that the star is best described as a heterogeneous mixture, described by bright spots and a dimmer photospheric component. We are proposing bright spots with higher temperatures, smaller covering factor, and shorter lifetimes than the steady spotted and facular components of Rackham et al. (2017) or Zhang et al. (2018). It is possible that the star has some combination of all of these features (bright spots, dark spots, and faculae).
In the previous section, we proposed that the apparent time-correlation between the occurrence of luminous flares and the brightness of the star is due to a physical association between the positions of bright spots which are associated with flares. An alternative hypothesis is that the ∼ 1% flux variations are not rotational modulation, but rather transient, bright active regions which accompany flares. If that is the case, then the 3.3 d periodicity in the K2 light curve should be interpreted as a characteristic active region timescale rather than a rotation period. In fact, if the bright spots were due to rotational variability, this might imply that the luminous flares preferentially occur at the same stellar longitude. This seems implausible given that magnetic activity on the surface of the star shouldn't be connected to the inertial frame. Thus, we feel that the possible correlation between the bright flares and spots may argue against the spot variation being due to rotation.
Hydrogen Balmer line emission (Hα) is common in late M dwarfs and variable Hα emission could potentially explain the 3.3-day modulations observed by Kepler; as Hα is not within the Spitzer bandpass, this would naturally explain the lack of a 3.3-day feature in the Spitzer lightcurve. The optical spectrum of TRAPPIST-1 of Burgasser et al. (2015) shows that the flux in Hα is only 0.3% of the flux integrated over the Kepler bandpass. However, Kruse et al. (2010) found that Hα equivalent widths vary by up to a factor of five for M8V stars, indicating that it is plausible that extreme Hα variability (by a factor of ∼ 5) could explain the Kepler variability. Further observations of Hα emission from TRAPPIST-1 are required to investigate this possibility.

CONCLUSIONS
We simultaneously model photometry of TRAPPIST-1 from Kepler and Spitzer to measure the properties of its putative starspots. We find that if the 3.3 day periodicity is due to stellar rotation, TRAPPIST-1 likely has a few bright spots rather than dark spots, and the spots have characteristic temperatures T ef f 5300 K and radii R spot /R ∼ 0.004.
The bright spots add a source of flux dilution to the transit light curves of each planet. We provide a correction factor for the transit depths of each planet, and propagate those depths to revise planet radii and densities; however, we note that other sources of variability and/or stellar inhomogeneity likely dominate this correction in the infrared (Delrez et al. 2018;Zhang et al. 2018).
We note that flares occur preferentially when the star is bright, and when the brightness is increasing most rapidly. This may suggest that that the flares are associated with the hot spots. Alternatively, the brightness variations could be the growth and decay of bright active regions on the stellar surface with a characteristic timescale of 3.3 days.
Though the nature of the proposed stellar activity is still unclear, the observations suggest that TRAPPIST-1 has bright spots rather than dark ones. Even if the continual spot variability observed by K2 were due to transient photospheric spots, rather than stellar rotation, the rapid appearance and disappearance of dark photospheric spots would produce a signal which was not observed by Spitzer.

A. THE EFFECT OF BRIGHT SPOTS ON TRANSIT LIGHT CURVES
Typically the transit depth is assumed to be the ratio of the cross-sectional areas of the projected planet on the star, where we use the equality sign because we are ignoring limb-darkening in this example. If there are bright unocculted spots on the star, the measured depth will be where c is the spot contrast relative to the photosphere. Rearranging, we find This correction for the unocculted bright spots in the Kepler and Spitzer bands are δ unspotted /δ spotted = 1.004 ± 0.001 and 1.00006 ± 0.00001, respectively. These transit depth dilution corrections allow us to update the observed planet radii reported by Gillon et al. (2017) with K2, which we list in Table 1.
If the are of bright spots is confined to the three spots we modeled, then the dominant systematic affecting planet radii from the Spitzer observations is the microvariability observed in Figure 1, rather than the bright spot variability. However, we note that this may be a lower limit on the effect of bright spots as our three-spot model only measures the variable component of the bright spots. We have carried out a test in which we added numerous small spots distributed in longitude, which produces an equally good fit to the variable light curve, but with much larger areal coverage of bright spots. The flux dilution due to the bright spots on TRAPPIST-1 will cause the transit depths to appear a shallower than they truly are. In Table 1 we list the revised planet properties -all revisions are within the uncertainties of the measurements. The effect as a function of wavelength is plotted in Figure 10.

B. POSTERIOR DISTRIBUTIONS
The full posterior distributions on all fit parameters in the bright spot model are shown in Figure 11.

C. PROBABILITY OF SPOT OCCULTATIONS
Occultations of the bright spots by the planets would allow us to infer the spot properties independently, as we showed for HAT-P-11 in Morris et al. (2017). We calculate an upper-limit on the probability that the proposed bright spots will be occulted by planets in the TRAPPIST-1 system. If we assume that the spots are infinitely long-lived, each spot is visible for half of the stellar rotation. We convert the posterior samples for spot latitudes to posterior samples of impact parameters p b , and integrate the posterior samples for b over the range of impact parameters occulted by each planet (b 0 , b 1 ), normalized by the integral of the posterior samples over all impact parameters, (this work) (G17) (this work) (G17) (this work) b 0.727 ± 0.009 0.729 ± 0.009 1.086 ± 0.035 1.090 ± 0.034 0.660 ± 0.560 0.656 ± 0.556 c 0.687 ± 0.010 0.690 ± 0.010 1.056 ± 0.035 1.060 ± 0.034 1.170 ± 0.530 1.159 ± 0.512 d 0.367 ± 0.017 0.368 ± 0.017 0.772 ± 0.030 0.775 ± 0.030 0.890 ± 0.600 0.882 ± 0.581 e 0.519 ± 0.026 0.521 ± 0.026 0.918 ± 0.039 0.921 ± 0.037 0.800 ± 0.760 0.793 ± 0.742 f 0.673 ± 0.023 0.676 ± 0.023 1.045 ± 0.038 1.049 ± 0.037 0.600 ± 0.170 0.589 ± 0.156 g 0.782 ± 0.027 0.785 ± 0.027 1.127 ± 0.041 1.131 ± 0.040 0.940 ± 0.630 0.927 ± 0.609 h 0.352 ± 0.033 0.353 ± 0.033 0.755 ± 0.034 0.759 ± 0.042 -- Table 1. Revised planet properties in the Kepler bandpass accounting for the flux dilution due to bright spots on TRAPPIST-1. Here we take planet masses and the definition of Depth = (R p /R ) 2 as in Gillon et al. (2017). The spot latitude solutions are degenerate: the spots could be on either side of the stellar equator and produce the same rotational modulation, and the planet impact parameter is similarly degenerate in that the planets could be occulting the northern or southern stellar hemisphere. As a result, we compute each probability of occultation twice: once assuming the planets occult the northern stellar hemisphere and once assuming they occult the southern hemisphere, and report the maximum spot occultation probability in Table 2.
The K2 light curve shows evolution of the flux modulation on timescales of days, indicating that spots are not infinitely long-lived. Thus these conservative upper limits are likely much greater than the probability of a spot occultation for transient spots on TRAPPIST-1.   Table 2. Maximum probability of a spot occultation for each planet in the TRAPPIST-1 system, assuming infinitely long-lived spots.