https://orcid.org/0000-0003-2528-3409
ABSTRACT
We use archival daily spot coverage measurements from Howard et al. to study the rotational modulation of the Sun as though it were a distant star. A quasi-periodic Gaussian process measures the solar rotation period Prot = 26.3 ± 0.1 d, and activity cycle period Pcyc = 10.7 ± 0.3 yr. We attempt to search for evidence of differential rotation in variations of the apparent rotation period throughout the activity cycle and do not detect a clear signal of differential rotation, consistent with the null results of the hare-and-hounds exercise of Aigrain et al. The full reconstructed solar light curve is available online.
1 INTRODUCTION
For decades astronomers have endeavoured to study the ‘Sun as a star’, measuring properties of the Sun that we typically measure on distant stars, with the goal of putting the Sun into context (e.g. Livingston 1991; Tayler 1996; Chaplin et al. 2004; Livingston et al. 2007; Hall et al. 2009; Bertello, Pevtsov & Pietarila 2012; Hall 2015; Egeland et al. 2017). These efforts are valuable, for example, for understanding the Sun’s activity through time, by observing Sun-like stars of different ages or at different phases in their activity cycles.
We are entering a new era for the study of rotational modulation of stars. Kepler has measured rotational modulation of tens of thousands of stars for four consecutive years, and K2 has measured rotation periods for many more stars, albeit over a shorter baseline. TESS will measure precision light curves for bright nearby stars, for a maximum duration of 355 consecutive days near the ecliptic poles in the primary mission (Ricker et al. 2014; Sullivan et al. 2015). Gaia will measure rotation periods for >105 stars (see e.g. Lanzafame et al. 2018). ESA’s PLAnetarty Transits and Oscillations (PLATO) mission may observe targets for up to 8 yr (Rauer et al. 2014), potentially allowing us to probe variations in the stellar rotational modulation of stars as a function of phase in their activity cycles. Having a solar benchmark light curve to compare these future, long-term light curves will be an important data product for the community.
Morris et al. (2018) developed tools for measuring the apparent stellar centroid offsets due to starspots that affect Gaia astrometry. In particular, a framework was developed for reconstructing archival spot maps of the Sun using the Mount Wilson Observatory (MWO) spot coverage catalogue published in Howard, Gilman & Gilman (1984). The MWO spot catalogue is a digitized representation of ‘white light’ photographic plate images of the solar disc taken from 1917 to 1985, denoting the apparent positions (latitude and longitude) and areas of penumbrae in sunspot groups. In this work, we use the same software and spot coverage archive as Morris et al. (2018), to reconstruct artificial time-series photometry of the Sun with one day cadence.
In Section 2, we introduce our approximation of the solar rotational light curve, and measure its properties as though it were a distant star. We will then recover several properties of the Sun using the reconstructed light curve. First and foremost we seek to recover the solar rotation period and activity cycle period, which are 25–34 d and 10.9 yr, respectively (Howe et al. 2000; Hathaway 2015).
We also follow the technique of Giles, Collier Cameron & Haywood (2017) to estimate the sunspot lifetimes. High-resolution observations of sunspots show that they have lifetimes ranging from hours to months (Solanki 2003). There is a roughly linear relationship between active region areas and their lifetimes, as described by Gnevyshev (1938) and Waldmeier (1955) (see also, e.g. Petrovay & van Driel-Gesztelyi 1997).
The broad range of possible rotation periods for the Sun is the result of differential rotation – the Sun rotates faster at the equator than at the poles (Miesch 2005). The pursuit to detect differential rotation from photometric rotational modulation of Sun-like stars in Kepler light curves has proven very difficult (Aigrain et al. 2015). Setting the perils aside, we will naively attempt to search for differential rotation by its effect on the solar light curve in Section 3.
2 THE SOLAR LIGHT CURVE
2.1 Constructing the light curve
We define Cartesian sky-plane coordinates (x, y), with the origin placed at the centre of the star, |$\hat{x}$| aligned with the stellar equator, and |$\hat{y}$| aligned with the stellar rotation axis. We describe each starspot with an ellipse with centroid ri = (xi, yi), and ri = |ri|. We can compute the negative flux contribution from each spot by computing the approximate spot area and contrast. A circular spot will be foreshortened near the stellar limb. The foreshortened circular spot can be approximated with an ellipse with semimajor axis Rspot and semiminor axis |$R_{\mathrm{spot}} \sqrt{1 - (r_i/\mathrm{R}_\odot)^2}$|.
Since these spots are small compared to the solar radius (Rspot/R⊙ < 0.1), we adopt one limb-darkened contrast for the entire spot, cld = (1 − c)I(r), where c is the flux contrast in the spot relative to the photosphere flux.
The spot group coverage catalogue of Howard et al. (1984) describes the daily areas and positions of sunspot groups from 1917 to 1985, see Fig. 1. We approximate each spot group as a single circular spot with the area of the entire spot group. We fix the spot contrast in the Kepler band at c = 1 − Ispot/Iphot = 0.7, which is the mean sunspot intensity averaged over the penumbra and umbra, assuming their typical penumbra covers a factor of 5 more area than the umbra (Solanki 2003).
The reconstructed solar light curve is shown in Figs 2 and 3. This very long-term view of the solar light curve shows periods of high variance separated by relatively quiet times, corresponding to the phase in the activity cycle. During solar maximum, there can be as many as 14 spot groups on the visible hemisphere of the Sun at once, leading to typical dips in flux of order ∼500 ppm. Near solar minimum, the spotless solar surface had no spot groups, and we have filled in those dates with no spot group entries with flux equal to unity.
Reconstructed solar light curve zoomed into cycle 19 to show fine structure. This light curve is unlike Kepler light curves for several reasons: our reconstruction has no photon noise, no p-mode oscillations, no granulation ‘flicker’, and no instrumental artefacts. In addition, unlike Kepler targets, we know the true unspotted flux of the Sun in these reconstructed light curves (i.e. we know when the Sun was truly spotless), so the light curve has a maximum of unity, rather than a median of unity.
The full reconstructed solar light curve is available online (Morris 2018).1
2.2 Constraining the effects of faculae
The Mount Wilson Observatory sunspot catalogue only tracked the positions and areas of dark sunspots, but did not measure the positions or sizes of faculae, which are small bright regions of concentrated magnetic flux. We reconstruct the solar light curve due to facular brightening using the same technique as in the previous section, but this time using faculae positions and areas from the Greenwich Photo-Heliographic Plate archive, digitized in 1999 by the NOAA Environmental Data Rescue Program, which provides facular positions and areas.
The resulting light curve of excess solar flux due to faculae is shown in Fig. 4. Typical brightening in the Kepler band due to faculae is small (≲20 ppm) compared to the darkening due to sunspots (≲200 ppm). Despite their large relative area coverage compared to starspots, the typical facular intensity contrast (c ∼ 1.05) is relatively small compared with sunspots (c ∼ 0.7), so we expect spots to dominate the rotational modulation of the Sun in the Kepler band, in agreement with Shapiro et al. (2016), for example.
Facular brightening of the Sun from the Greenwich Photo-Heliographic archive spanning three activity cycles. Note that the scale of facular brightening in the Kepler band is small (≲20 ppm) compared to the darkening due to sunspots (≲200 ppm).
The dominance of sunspots over faculae in the rotational light curve of the Sun is not to be confused with the fact that the Sun is considered ‘faculae dominated’ on time-scales of the activity cycle. That is, near solar maximum the Sun is bolometrically brighter than it is at solar minimum (Solanki, Krivova & Haigh 2013). What we refer to as the solar light curve in this work is not the bolometric flux of the Sun, rather it is the flux integrated over a bandpass like those of Kepler, TESS, or Gaia. As such, we choose to ignore the effects of faculae in the remainder of this work, since spots dominate the rotational modulation, which is our primary focus.
2.3 Measuring the solar rotation and activity cycle periods
An astronomer’s first instinct is likely to measure periodicities upon seeing a light curve such as Fig. 2. In this section, we examine the autocorrelation function and Lomb–Scargle (LS) periodogram of the solar light curve to establish benchmark measurements of the rotation and activity cycle periods.
2.3.1 Gaussian process regression
The autocorrelation function of the solar light curve is shown in Fig. 5. There is short-term variation peaking at 26 d – approximately the rotation period of the Sun at the photosphere near the active latitudes (Howe et al. 2000). There is also a long-term decaying cosine-shaped correlation with its first peak at 10.6 yr, corresponding to the magnetic activity cycle period of ∼11 yr (Hathaway 2015). Finally, there is a cosine-shaped correlation with a period of 365 d, corresponding to the orbital period of the Earth. This systematic crops up because the Earth’s orbit is inclined with respect to the solar equator by 7.25° (Meeus 1991), causing starspots to drift slightly towards and away from the solar equator throughout each year, injecting a small correlated signal into the reconstructed light curve.
Upper left: Long-term autocorrelation function of the solar light curve in Fig. 2. The first peak is ∼10.5 yr. Upper right: Short-term signals in the autocorrelation function of the solar light curve, with a peak at 26 d. This estimate of the rotation period is approximately consistent with the rotation period at the active latitudes (Howe et al. 2000). See Section 2.3 for a more robust rotation period measurement using Gaussian process regression with a quasi-periodic kernel. Lower: Posterior distributions for the magnetic activity cycle period Pcyc and the solar rotation period Prot measured with a quasi-periodic Gaussian process regression to the synthetic photometry from 1917 to 1985 (see Fig. 2). We measure Prot = 26.32 ± 0.14 d and Pcyc = 10.61 ± 0.23 yr.
We also measure activity cycle period Pcyc = 10.61 ± 0.23 yr. This is consistent with canonical cycle period measured by taking the dates of the minima of cycle 1 and cycle 23 and dividing by 22, yielding an average cycle period of 10.9 yr (131.7 months, Hathaway 2015).
2.3.2 Lomb–Scargle periodogram
Next, we use the LS periodogram to compare its ability to pick out the quasi-periodic peak – see Fig. 6. The dominant period is 27 d, just longer than the rotation period measured by Gaussian process regression in the previous section. The difference in rotation periods measured with each technique is an artefact of the intrinsically quasi-periodic nature of the Gaussian process kernel in the previous section, and the strict periodicity enforced by the LS periodogram. In addition, the uncertainty in the periodicity measured with the LS periodogram is not well defined, so it is not possible to do a robust comparison between the LS and Gaussian process period measurements. Turning to longer periods, the activity cycle peak is prominent at 10.6 yr – this result is consistent with the Gaussian process regression measurement.
We prefer the value from the quasi-periodic Gaussian process analysis for the apparent rotation period rather than the LS period because (1) we know from high-resolution observations that there’s more than one frequency at play – for example, starspots emerge and decay on time-scales similar to the stellar rotation period and (2) the Gaussian process regression provides us with robust uncertainties on the apparent rotation period. For these reasons, we encourage observers of distant stars to consider using Gaussian process regression over the LS periodogram when searching for the rotation period at the mean active latitudes (VanderPlas 2018).
2.4 Measuring active region evolution time-scales
The autocorrelation function in Fig. 5 has a peak at the rotation period of the star, and smaller peaks at integer multiples of the rotation period with decreasing amplitudes. Giles et al. (2017) developed a technique for measuring active region evolution time-scales by modelling the autocorrelation functions of active stars with a underdamped simple harmonic oscillator (uSHO), which we apply here to the autocorrelation function of the reconstructed solar light curve to estimate active region lifetimes.
We attempt to measure the sunspot lifetimes from the autocorrelation function. However, the form of the autocorrelation function (Fig. 5) is slightly different to those seen in Giles et al. (2017), which typically follow the pattern of a uSHO. In Fig. 5, there is an additional decreasing trend which causes the subsequent peaks to be significantly lower than the central peak, which arises from the much longer time-scale activity cycle pattern. This effect is persistent whether we generate the autocorrelation function for the light curve as a whole, or cut it up into smaller portions and combine the autocorrelation functions.
Although the uSHO fits were unsuccessful, we can still make some qualitative statements from inspection of the autocorrelation function at short lags. The signal of rotation peaking at 26 d has repeated aliases at twice and possibly at three times the rotation period, each with diminished amplitude, before the aliases of the rotation signal appear to decay away at large lags (≳3Prot). This suggests that the typical spot decay time-scale is similar to the rotation period, and only occasional spots survive more than one or two solar rotations. This observation is in agreement with spatial resolved observations which show that the longest-lived sunspots live of order several rotations (Pettit 1951; Howe et al. 2000), but most only survive for less than one rotation (Petrovay & van Driel-Gesztelyi 1997).
3 DIFFERENTIAL ROTATION
Many efforts have been made to quantify differential rotation in Kepler light curves of stars, most notably in Aigrain et al. (2015), where several groups of observers attempted to measure the differential rotation rate in synthetic light curves. The authors found that there was little relation between the injected and recovered differential rotation rates, indicating that Kepler detections of solar-like differential rotation ought to be treated with caution.
In this section, we set out to mimic this perilous exercise by measuring the solar rotation period in consecutive one year bins, using the quasi-periodic Gaussian process technique that we used in Section 2.3 to measure the rotation period of the full light curve. We choose one year bins so that there is sufficiently long baseline to get a fit for the period, but the duration is short compared to the activity cycle period (11 yr). If the rotational modulation contains the signature of differential rotation, we expect to find that the apparent rotation period changes slightly from one year to the next, as spots emerge at different latitudes throughout the activity cycle, and due to differential rotation, the spots rotate with slightly different periods.
Ideally, we would observe that at the beginning of each activity cycle, the spots emerge at high latitudes and therefore the apparent rotation period is long. Then as the activity cycle progresses, spots emerge at lower latitudes, revealing shorter rotation periods.
The rotation period recovered from fitting the quasi-periodic Gaussian process to one year bins of the solar light curve is shown in black points in Fig. 7. The red curve shows the rotation period at the mean area-weighted spot latitude averaged into yearly bins, and shows the small differential rotation signal imparted by the activity cycle which we are attempting to measure. In practice, we observe a spread in measured rotation periods much larger than the variance due to the activity cycle, with similarly large uncertainties. Activity minima can be identified in this figure by the large uncertainties on the rotation period, when few spots are present to drive rotational modulation. In between these points of large uncertainties are intervals where the rotation period is measured more precisely, though it is roughly consistent with a 26.3 d rotation period throughout all phases of the activity cycle. Assuming the rotation period is 26.3 d throughout, the reduced |$\tilde{\chi }^2 = 8$|, indicating that the variance is indeed greater than expected for Gaussian-distributed errors. However, the stochastic nature of the measurements makes it impossible to recover the true differential rotation rate from these rotation measurements. Therefore even at ‘infinite’ signal to noise, we arrive at the same conclusion as Aigrain et al. (2015) – measuring differential rotation shear from rotational modulation alone is a fraught exercise.
Measured solar rotational period inferred from quasi-periodic Gaussian process regression to 1-yr-long bins of the solar light curve (black circles), compared with the rotation period at the mean area-weighted spot latitude averaged in 1 yr bins (red curve). As spots emerge at different latitudes with differential rotation, we hoped to find that the rotation period varied from year to year with the phase of the activity cycle – with spots emerging at high latitudes and long rotation periods, and the rotation period appearing to decrease as spots emerge closer to the solar equator. It appears that due to spots emerging at a broad range of latitudes at all phases of the activity cycle, the apparent rotation period remains largely constant, irrespective of the activity cycle phase.
4 DISCUSSION
The detection of differential rotation from the solar light curve eludes us in Section 3. One reason for this is made clear by the butterfly diagram in Fig. 1 – the distribution of spots within active latitudes of the Sun are broad; spots are distributed within ±8° of the mean ‘active latitude’ at each phase of the activity cycle. The spots at multiple latitudes each contribute to the rotational modulation with their own rotation period, imprinting the mean rotation period on the light curve, rather than the specific rotation period at a high or low latitude. We are encouraged by recent work by Benomar et al. (2018) which may hold the key to measuring differential rotation from stellar photometry via asteroseismology for at least a small subset of stars.
One limitation of this reconstruction approach is that the Howard et al. (1984) spot archive only catalogued spots within 60° longitude of the central solar meridian, meaning that spots on the limb were not logged. If spots near the limb were included in this time series, the overall flux trends might be smoother, and there would be fewer days with flux equal to unity. However, the net effect on the rotation period and activity cycle measurements is likely small, since spots on the limb are geometrically foreshortened, and due to the Wilson depression, they have smaller contrasts than spots at disc centre (Solanki, Walther & Livingston 1993).
5 CONCLUSION
We reconstructed a one day cadence light curve of the Sun using the sunspot archive of Howard et al. (1984). We compared the amplitude of variability due to dark sunspots to the amplitude of brightening from faculae, and found that the dark sunspots dominate the rotational modulation.
With the noise-free light curve, we measured rotation period and activity cycle period of the Sun with both the LS periodogram and a quasi-periodic Gaussian process regression. The rotation periods and activity cycle periods measured with both techniques are consistent with the rotation period at the active latitudes, and the duration of a typical activity cycle. We showed that differential rotation cannot be detected even from this idealized reconstructed light curve.
ACKNOWLEDGEMENTS
We gratefully acknowledge the following software packages which made this work possible: astropy (Astropy Collaboration 2013; The Astropy Collaboration 2018), ipython (Perez & Granger 2007), numpy (Van Der Walt, Colbert & Varoquaux 2011), scipy (Jones et al. 2001), matplotlib (Hunter 2007), celerite (Foreman-Mackey et al. 2017), gatspy (Vanderplas 2015), emcee (Foreman-Mackey et al. 2013). Some of this work has been carried out in the frame of the National Centre for Competence in Research ‘PlanetS’ supported by the Swiss National Science Foundation (SNSF).
JRAD acknowledges support from the DIRAC Institute in the Department of Astronomy at the University of Washington. The DIRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences, and the Washington Research Foundation. This research has made use of NASA’s Astrophysics Data System.
Footnotes
REFERENCES
Author notes
DIRAC Fellow.
