A Relationship Between Stellar Age and Spot Coverage

In this section we define what we call the “smoothed amplitude” of each light curve. This quantity was first published by Douglas et al. 2017 for Praesepe members. The smoothed amplitude is the difference between the maximum and minimum flux after the light curve has been phase-folded and smoothed with a Gaussian kernel.

USCO is a $10\pm 2$ Myr old part of the nearby Sco-Cen star-forming region (Pecaut & Mamajek, 2016). We queried for K2 photometry from FGK members of the young association listed by Gagné et al. (2018), and found 19 sources. We measure the stellar rotation period for each star by estimating the peak power in the Lomb-Scargle periodogram (Lomb, 1976; Press & Rybicki, 1989). We then phase-fold each light curve on the best period, smooth the light curve with a Gaussian kernel of width 50-cadences, and report smoothed amplitudes. We visually inspected each light curve for hints of binarity in the periodogram (multiple, non-aliased periods), and discarded any possible binaries and stars with ambiguous rotation periods.

UCL is a $16\pm 2$ Myr old part of the nearby Sco-Cen star-forming region (Pecaut & Mamajek, 2016). We queried for TESS photometry of sources listed as UCL members by Gagné et al. (2018), and found 34 sources with TESS Input Catalog (TIC) masses $M>0.6M_{\odot}$ in the full-frame images (FFIs). We query the FFI database for a square region 10 pixels per side, centered on the coordinates of each UCL member. We subtracted the median flux in a 3 pixel radius circular aperture from each FFI, and remove a quadratic trend from each FFI light curve. As in the previous section, for each light curve, we measure the stellar rotation period for each star by estimating the peak power in the Lomb-Scargle periodogram, limiting the maximum period to half of the TESS sector duration. We phase-fold each light curve on the best period, smooth the light curve with a Gaussian kernel of width 50-cadences, and report the smoothed amplitude. Again, we visually inspected each light curve for signs of binarity in the periodogram, and discarded any possible binaries. See Figure 1 for a visual representation of this process.

Psc-Eri is a 120 Myr old stellar stream extending 120^∘ across the sky (Meingast et al., 2019). We followed a similar procedure to the previous subsection to produce light curves for each of 100 FGK targets which were identified as members of the Psc-Eri stream by Curtis et al. (2019b). In addition to using their membership list, we also used the rotation periods reported by Curtis et al. (2019b), and simply measured the smoothed amplitudes of each light curve after phase folding the light curve and smoothing with a Gaussian kernel with width 100 cadences.

Praesepe is a well-studied, nearby, 650 Myr old cluster. Douglas et al. (2017) measured the amplitudes of rotational modulation of many apparently single stars in Praesepe with K2. Here we will focus on stars that are not classified as binaries or blends. The authors calculated “smoothed amplitudes” (which are reported as semi-amplitudes) of the rotational modulation for each star, in which the maximum and minimum flux are measured in smoothed, phase-folded light curves. We adopt the authors’ smoothed amplitudes and rotation periods without modification for 220 FGK stars in Praesepe.

NGC 6811 is a 1 Gyr old cluster in the Kepler field (Meibom et al., 2011). Curtis et al. (2019a) found what they called a temporary epoch of stalled spin-down for low-mass stars in NGC 6811. We use the Curtis et al. (2019a) membership list and rotation periods to build a sample of 167 FGK stars in NGC 6811, and measure smoothed amplitudes from the PDCSAP fluxes for each star in Kepler Quarter 2 using a Gaussian kernel with width 100 cadences. We select only the second quarter of Kepler observations to more closely approximate the variability on timescales similar to the K2 and TESS observations; including the full Kepler light curve tends to average over many spot evolutions and decrease smoothed amplitudes.

M67 is a $4.2\pm 0.2$ Gyr old cluster in K2 Campaign 5 (Gonzalez, 2016a, b). We use the cluster membership and rotation periods for 96 stars listed in Gonzalez (2016b) (which are in good agreement with the periods of Barnes et al., 2016), and measure smoothed amplitudes from the PDCSAP fluxes for each star using a Gaussian kernel with width 250 cadences.

Figure 2 shows a trend in the observable properties of the light curves: there is an anti-correlation between the typical rotation periods of stars in each cluster and the smoothed amplitude of the light curves. One useful perspective encoded in this plot has to do with the axisymmetry of the starspot distributions. If a star has several starspots which are distributed uniformly in longitude, the rotational modulation amplitude will be relatively small; whereas if spots are concentrated into a small region on one stellar hemisphere, the rotational modulation amplitude will be relatively large. Young stars have short rotation periods and large smoothed amplitudes, corresponding to significant concentrations of dark spots, or significant deviations from uniformly distributed spots. As stars age they drift towards the upper left of the plot (following the direction of the silver arrow); their rotation periods increase and their smoothed amplitudes decrease, or their spots become distributed more uniformly.

Figure 3 shows another view of the stellar samples in Figure 2, demonstrating the decay of the smoothed amplitude as a function of cluster age. The linear regression trend line (gray) indicates that smoothed amplitudes of light curves generally decline with increasing stellar age, as

where and $\alpha={0.56}^{+1.00}_{-0.31}$ and $m=-0.50\pm 0.17$.

This power-law index $m$ is remarkably close to the $t^{-1/2}$ decline in chromospheric emission with age discovered by Skumanich (1972). Perhaps this result suggests there may be a simple relationship between the area in chromospherically active regions and the area in starspots, causing this simple metric for the spot coverage, the smoothed amplitudes, to have the same age dependence as the chromospheric emission index (such relations already exist for magnetic field strength and Ca emission, for example: Schrijver et al., 1989).

We simulate ensembles of light curves of stars through a full rotation and measure their smoothed amplitudes using fleck ¹¹1https://github.com/bmorris3/fleck. fleck is a pure Python software package which simulates starspots as circular dark regions on the surfaces of rotating stars, accounting for foreshortening towards the limb, and limb darkening, which is an efficient, vectorized iteration of earlier codes used in Morris et al. (2018, 2019). The fleck algorithm is outlined as follows: suppose we have $N$ stars, each with $M$ starspots, distributed randomly above maximum latitudes ${\rm\ell_{max}}$, observed at $N$ inclinations $\vec{i}_{\star}$ (one unique inclination per star), observed at $O$ phases throughout a full rotation $\vec{\phi}\sim\mathcal{U}(0,90^{\circ})$.

We initialize each star such that its rotation axis is aligned with the $\hat{z}$ axis, and set the observer at $x\rightarrow\infty$, viewing down the $\hat{x}$ axis towards the origin.

We define the rotation matrices about the $\hat{y}$ and $\hat{z}$ axes for a rotation by angle $\theta$:

We begin with the matrix of starspot positions in Cartesian coordinates ${\bf C_{i}}$,

for $i=1$ to $N$ with shape $(3,M)$, which we collect into the array ${\bf S}$,

of shape $(3,M,N)$. We rotate the starspot positions through each angle in $\phi_{j}$ for $j=1$ to $O$ by multiplying $\bf S$ by the rotation array

with shape $(O,1,1,3,3)$. Using Einstein notation, we transform the Cartesian coordinates array $\bf C$ with:

to produce a array with shape $(3,O,M,N)$, where $lm…$ indicates an optional additional set of dimensions. Then after each star has been rotated about its rotation axis in $\hat{z}$, we rotate each star about the $\hat{y}$ axis to represent the stellar inclinations $i_{\star,k}$ for $k=1$ to $N$, using the rotation array

with shape $(N,3,3)$, by doing

which produces another array of shape $(3,O,M,N)$. Now we extract the second and third axes of the first dimension, which correspond to the $y$ and $z$ (sky-plane) coordinates, and compute the radial position of the starspot ${\bf\rho}=\sqrt{y^{2}+z^{2}}$, where $\bf\rho$ has shape $(O,M,N)$. We now mask the array so that any spots with $x<0$ are masked from further computations, as these spots will not be visible to the observer. We use $\rho$ to compute the quadratic limb darkening

for $\mu=\sqrt{1-\rho^{2}}$. We compute the flux missing due to starspots of radii $\bf R_{\rm spot}$, which has shape $(M,N)$:

The unspotted flux of the star is

so the spotted flux is

The model presented above works best for spots that are small. The array masking step for $x<0$ does not account for the change in stellar flux due to large starspots which straddle the limb of the star. Large starspots also have differential limb-darkening across their extent, which is not computed by fleck.

Comparison with STSP²²2https://github.com/lesliehebb/STSP is shown in Figure 4. The models reproduce consistent rotational light curves at the 50 ppm level – similar to the Kepler noise floor on 1 hour timescales for bright stars (Borucki et al., 2011; Christiansen et al., 2012). The maximum divergence between models occurs when the spots are near the limb, where STSP accounts for the spot which straddles the limb and fleck does not. The differences between the models are small compared to the uncertainties in flux of the K2 photometry, for instance.

The relationship between starspot coverage and age can be deduced from the slope and intercept in the right panel of Figure 7. We find the hemispheric starspot covering fraction $f_{S}$ is related to the stellar age $t$ in Gyr by the simple relation

where $a={0.014}^{+0.022}_{-0.008}$ and $n=-0.37\pm 0.16$. The power law index $n$ is statistically consistent with the approximate inverse square root relation between chromospheric activity and stellar age by Skumanich (1972).

One must take care not to infer upper or lower limits to the spotted coverage of individual stars of a given age from Equation 10. The light curve ensemble modeling technique constrains a typical spot coverage for stars in each sample, and outliers are likely to exist which will not fall neatly within the confidence intervals of Equation 10.

LkCa 4, for example, is a weak-lined T Tauri star in the Taurus Molecular Cloud ($<10$ Myr; Kenyon & Hartmann, 1995). High-resolution near-infrared IGRINS spectra from Gully-Santiago et al. (2017) constrained the star’s spot coverage to $f_{S}\sim 0.8$. The authors argue that such a large coverage by dark regions challenges our notion of “spot coverage,” since the majority of the stellar photosphere emits at a cooler temperature than its spectral type (Pettersen et al., 1992). The ensemble light curve modeling technique assumed three starspots which cover a minority of the stellar surface, so the results (Equation 10) should not be applied to stars covered in mostly cool regions, like LkCa 4.

In this subsection we discuss several factors which may have small but significant affects on the conclusions drawn from the ABC analysis.

In addition to ensemble light curve modeling, we can also use fleck to directly model the photometry of individual stars. In this section, we fit the fleck model to the light curves of spotted stars with Markov Chain Monte Carlo (MCMC) to validate the spot coverage relation in Equation 10 for stars with precise photometry and ages, listed in Table 5.3.

As discussed earlier, the positions, radii, and contrasts of starspots are degenerate with one another, making it difficult to extract precise spot properties from the rotational modulation of planet hosting stars. Thus we make several simplifying assumptions that allow us to fit for the spot coverage on these stars. First, we assume that the stellar inclination is $i_{\star}=90^{\circ}$, which may be a good approximation since each of these systems host (often multiple) transiting exoplanets. Next, we fix the spot contrast to $c=0.7$; this is compatible with the area-weighted spot coverage of sunspots, the starspot contrasts of HAT-P-11 (Morris et al., 2017), and a valid approximation to the observed spot contrasts of several stars extrapolated into the Kepler and TESS bandpasses (Morris et al., 2018). We also place a uniform bounded prior on the spot latitudes $|\ell|<60^{\circ}$. It is necessary to impose this latitude prior because without it the Markov chains occasionally prefer spots with implausibly large radii, skewing the fits towards large spot coverage; this prior is consistent with by solar observations since sunspots are not observed above $\sim 45^{\circ}$ (Morris et al., 2019). Finally, we also enforce the prior that $\rm 0<R_{spot}/R_{star}<1$, ensuring that the largest spots are still small compared to the stellar radius, to forbid spots with radii several orders of magnitude larger than the largest sunspot groups.

We choose to fix the number of starspots to three in each fit. We find that fitting less than three spots does not accurately reproduce the observed rotational modulation, and fitting additional spots does not significantly improve the fits.

The photometry, maximum-likelihood models, spot coverage posterior distributions, and spot maps are shown in Figure 8. The three-spot model reasonably approximates the rotational modulation in each light curve. Then we compare the posterior distributions for the expected spot coverage inferred from the fleck models with the spot coverage estimate from Equation 10 in Figure 9. Most stars fall within the $1\sigma$ confidence interval for the predicted spot coverage.

We include an upper-limit for the spot coverage of the Sun in Figure 9 using the Mount Wilson Observatory sunspot catalog (Howard et al., 1984), analyzed in the stellar context by Morris et al. (2017).

Equation 10 tends to over-predict the amplitudes of rotational modulation in the oldest stars: Kepler-21, Kepler-50 and the Sun (Figure 9). Some possible explanations for overestimate from include: (1) cluster stars like those of M67, which are the anchors for the spot coverage relation near 4 Gyr, are more active than field stars like Kepler-21 and Kepler-50; (2) spots are more uniformly distributed on old field stars, diminishing the rotational modulation amplitude; (3) the field stars could have been observed near activity minimum; or (4) we may be viewing Kepler-21 and Kepler-50 at low inclinations, producing small light curve amplitudes. A larger sample of $>4$ Gyr stars with precise ages and clear rotational modulation is required to determine whether or not a break in the power-law near $1$ Gyr is justified.

Exoplanet radii are measured from the depths of their transit light curves. Unocculted dark starspots slightly increase the depths of transit light curves, giving the appearance of slightly larger planets. In this section we quantify the extent of exoplanet radius amplification by dark starspots.

The Mandel & Agol (2002) transit model for a uniform source is computed

where $\lambda^{e}$(t) is the fraction of the host star eclipsed. However, this is assuming the eclipsed body is not changing in brightness. If the host star is changing in brightness, Equation 13 becomes

where $F_{\star}(t)$ is the brightness of the star as a function of time relative to the unspotted stellar flux. In practice, observations of individual transits are often normalized by the out-of-transit flux immediately preceding and following the transit event, so what is observed is

For a star with quadratic limb-darkening, for example, the maximum flux obscured during the transit event is

where $b$ is the impact parameter, so the minimum observed flux during a transit event $F_{\rm obs,min}$ is

Since a spotted star has $F_{\star}(t)\leq 1$, it is clear that the transit depth $\delta_{\rm obs}=\lambda^{e}_{\rm max}/F_{\star}$ is amplified by a dimmer, more spotted surface. It is straightforward to compute $F_{\star}(t)$ for spotted stars with fleck, so we can easily compute $\delta_{\rm obs}$ given parameters for the planet and starspots.

We simulate a star with three large spots of radius $\rm R_{spot}/R_{\star}=0.4$, with stellar inclination $i_{\star}=90^{\circ}$, and a planet with $\rm R_{p}/R_{\star}=0.1$ at impact parameter $b=0$ with solar quadratic limb-darkening parameters. This is equivalent to the very large spot coverage $f_{S}=0.12$. The results are shown in Figure 10. The transit depth can be amplified by as much as $10\%$ by the modulation due to starspots. This is large enough to be detected given typical observational uncertainties on the radii of small planets orbiting FGK stars, so the variation in transit depth as a function of time is likely an important factor for accurately measuring exoplanet radii orbiting stars with ages $\lesssim 20$ Myr.