Inferring probabilistic stellar rotation periods using Gaussian processes

TL;DR

This paper introduces a Gaussian Process-based method to accurately infer stellar rotation periods from light curves, outperforming traditional techniques and providing valuable posterior distributions for astrophysical studies.

Contribution

The paper presents a novel Gaussian Process approach with a quasi-periodic kernel for inferring stellar rotation periods, enabling more accurate and probabilistic analysis compared to existing methods.

Findings

Gaussian Process method outperforms sine-fitting and autocorrelation methods

Successfully inferred rotation periods for simulated and real Kepler data

Provides posterior distributions for rotation periods facilitating hierarchical studies

Abstract

Variability in the light curves of spotted, rotating stars is often non-sinusoidal and quasi-periodic --- spots move on the stellar surface and have finite lifetimes, causing stellar flux variations to slowly shift in phase. A strictly periodic sinusoid therefore cannot accurately model a rotationally modulated stellar light curve. Physical models of stellar surfaces have many drawbacks preventing effective inference, such as highly degenerate or high-dimensional parameter spaces. In this work, we test an appropriate effective model: a Gaussian Process with a quasi-periodic covariance kernel function. This highly flexible model allows sampling of the posterior probability density function of the periodic parameter, marginalising over the other kernel hyperparameters using a Markov Chain Monte Carlo approach. To test the effectiveness of this method, we infer rotation periods from 333 simulated stellar light curves, demonstrating that the Gaussian process method produces periods that are more accurate than both a sine-fitting periodogram and an autocorrelation function method. We also demonstrate that it works well on real data, by inferring rotation periods for 275 Kepler stars with previously measured periods. We provide a table of rotation periods for these 1132 Kepler objects of interest and their posterior probability density function samples. Because this method delivers posterior probability density functions, it will enable hierarchical studies involving stellar rotation, particularly those involving population modelling, such as inferring stellar ages, obliquities in exoplanet systems, or characterising star-planet interactions. The code used to implement this method is available online.

Tables

Table 1. Table 1: Ranges and distributions of parameters used to simulate light curves in Aigrain et al. ( 2015 )

Parameter	Range	Distribution
Rotation period, $P_{rot}$	10 - 50 days (90%)	log uniform
	1 - 10 days (10%)	log uniform
Activity cycle length	1 - 10 years	log uniform
Inclination	0 - 90∘	Uniform in ${\sin }^{2}i$
Decay timescale	(1 - 10) $\times P_{rot}$	log uniform

Table 2. Table 2: Priors and bounds on the natural logarithms of the GP model parameters.

Parameter	Prior	Bounds
$\ln A$	$𝒩(-13,5.7)$	(-20, 0)
$\ln l$	$𝒩(7.2,1.2)$	(2, 20)
$\ln \Gamma$	$𝒩(-2.3,1.4)$	(-10, 3)
$\ln \sigma$	$𝒩(-17,5)$	(-20, 0)
$\ln P$	Uniform / ACF-based	($\ln 0.5,\ln 100$)

Table 3. Table 3: Median absolute (MAD), median relative (MRD) deviations and root-mean-squared (RMS) for the LS periodogram, McQuillan et al. ( 2013a ) ACF and GP period recovery methods.

Method	MAD	MRD	RMS
GP (acf prior)	0.40 days	1.97%	0.25 days
GP (uninformative prior)	0.48 days	2.62%	0.44 days
ACF	0.37 days	1.75%	4.58 days
LS periodogram	0.52 days	2.82%	1.04 days

Table 4. Table 4: Priors and bounds on the natural logarithms of the GP model parameters, for Kepler  light curves

Parameter	Prior	Bounds
$\ln A$	$𝒩(-13,5.7)$	(-20, 0)
$\ln l$	$𝒩(5.0,1.2)$	(2, 8)
$\ln \Gamma$	$𝒩(1.9,1.4)$	(0, 3)
$\ln \sigma$	$𝒩(-17,5)$	(-20, 0)
$\ln P$	ACF-based	($\ln 0.5,\ln 100$)

Equations

$$k_{i,j}=A\exp \left(-\frac{(x_i-x_j{)}^2}{2l^2}\right).$$

$$k_{i,j}=A\exp \left[-\frac{(x_i-x_j{)}^2}{2l^2}-{\Gamma }^2{\sin }^2\left(\frac{\pi (x_i-x_j)}{P}\right)\right]+{\sigma }^2{\delta }_{ij}.$$

$$p(\mathbf{\theta }|y)\propto \mathcal{L}(\mathbf{\theta })p(\mathbf{\theta }),$$

$$\ln \mathcal{L}=-\frac{1}{2}\sum _{n=1}^N\left[\frac{(y_n-\mu {)}^2}{{\sigma }_n^2}+\ln (2\pi {\sigma }_n^2)\right],$$

$$\ln \mathcal{L}=-\frac{1}{2}{\mathbf{r}}^{\mathbf{T}}{\mathbf{C}}^{-1}\mathbf{r}-\frac{1}{2}\ln |\mathbf{C}|+\frac{\mathbf{N}}{2}\log 2\pi ,$$

$$\mathbf{C}$$

$$y=A{e}^{-t/\tau }\cos \frac{2\pi t}{T_i},$$

$$Q_i=\left(\frac{{\tau }_i}{T_i}\right)\left(\frac{N_i{h}_i}{R_i}\right),$$

$$\begin{gathered} p(\ln P)= \\ \frac{\sum _i{Q}_i\left(0.9\mathcal{N}(\ln T_i,0.2)+0.05\mathcal{N}(\ln (T_i/2),0.2)+0.05\mathcal{N}(\ln (2T_i),0.2)\right)}{{\sum }_i{Q}_i}. \end{gathered}$$

$${\sigma }_{\nu }=\frac{3\pi {\sigma }_N}{2N^{1/2}TA},$$

$$k_{i,j}=A\exp \left(-\frac{(x_i-x_j{)}^2}{2l^2}\right)\cos \left(\frac{\pi (x_i-x_j)}{P}\right).$$

Full text

††pubyear: 2017††pagerange: Inferring probabilistic stellar rotation periods using Gaussian processes–LABEL:lastpage

1 Introduction

The disk-integrated flux of a spotted, rotating star often varies in a non-sinusoidal and Quasi-Periodic (QP) manner, due to active regions on its surface which rotate in and out of view. Complicated surface spot patterns produce non-sinusoidal variations, and the finite lifetimes of these active regions and differential rotation on the stellar surface produce quasi-periodicity (Dumusque et al., 2011). A strictly periodic sinusoid is therefore not necessarily a good model choice for these time-series. A physically realistic model of the stellar surface would perfectly capture the complexity of shapes within stellar light curves as well as the quasi-periodic nature, allowing for extremely precise probabilistic period recovery when conditioned on the data. However, such physical models require many free parameters in order to accurately represent a stellar surface, and some of these parameters are extremely degenerate (e.g. Russell, 1906; Jeffers and Keller, 2009; Kipping, 2012). In addition to global stellar parameters such as inclination and rotation period, each spot or active region should have (at minimum) a longitude, latitude, size, temperature and lifetime. Considering that stars may have hundreds of spots, the number of free parameters in such a model quickly becomes unwieldy, especially to explore its posterior Probability Density Function (PDF). Simplified spot models, such as the one described in Lanza et al. (2014) where only two spots are modelled, have produced successful results; however, such relatively inflexible models sacrifice precision.

Standard non-inference based methods to measure rotation periods include detecting peaks in a Lomb-Scargle (Lomb, 1976; Scargle, 1982) (LS) periodogram (e.g. Reinhold et al., 2013), Auto-Correlation Functions (ACFs), (e.g. McQuillan et al., 2013b) and wavelet transforms (e.g. García et al., 2014). The precisions of the LS periodogram and wavelet methods are limited by the suitability of the model choice: a sinusoid for the LS periodogram, and a choice of mother wavelet, assumed to describe the data over a range of transpositions and scales (see, e.g. Carter and Winn, 2010), for the wavelet method. In contrast, since it does not rely on a generative model, the ACF method is much better suited to signals that are non-sinusoidal. In fact, as long as the signal is approximately periodic the ACF will display a peak at the rotation period, no matter its shape. A drawback of the ACF method however, is that it requires data to be evenly-spaced. Most ground-based measurements, such as the future Large Synoptic Survey Telescope (LSST), do not nearly satisfy this criterion, and even light curves from the Kepler space telescope, which has provided the richest stellar rotation data to date, can only be approximated as uniformly sampled. An ACF is also an operation performed on the data rather than a generative model of the data, and so is not inherently probabilistic. This means that the effects of the observational uncertainties cannot be formally propagated to constraints on the rotation period.

In this work, we test an effective model for rotationally modulated stellar light curves which captures the salient behaviour but is not physically motivated — although some parameters may indeed be interpreted as physical ones. An ideal effective model should have a small number of non-degenerate parameters and be flexible enough to perfectly capture non-sinusoidal and QP behaviour. A Gaussian process (GP) model fulfills these requirements. We thus use a GP as the generative model at the core of a method to probabilistically infer accurate and precise stellar rotation periods. This enables us to estimate the posterior PDF of the rotation period, thereby producing a justified estimate of its uncertainty. Using a quasi-periodic GP model to infer a rotation period is not a new idea. For example, we used this method in Vanderburg et al. (2015) to measure the rotation period of an exoplanet host. Previously to that, Haywood et al. (2014) used the quasi-periodic GP model to disentangle the correlated noise produced by stellar activity and rotation from the radial velocity signature of an exoplanet. In addition, Littlefair et al. (2017) use a quasi-periodic GP to establish the quasi-periodicity of variability in the light curves of brown dwarfs. What we present here is a test of the GP method and a comparison with alternative rotation period measurement methods. We also provide code allowing others to easily apply this technique themselves.

GPs are commonly used in the machine learning community and increasingly in other scientific fields such as biology, geophysics and cosmology. More recently, GPs have been used in the stellar and exoplanet fields within astronomy, to capture stellar variability or instrumental systematics (see e.g. Gibson et al., 2012; Haywood et al., 2014; Dawson et al., 2014; Barclay et al., 2015; Haywood, 2015; Evans et al., 2015; Rajpaul et al., 2015; Czekala et al., 2015; Vanderburg et al., 2015; Rajpaul et al., 2016; Aigrain et al., 2016; Littlefair et al., 2017). They are useful in regression problems involving any stochastic process, specifically when the probability distribution for the process is a multi-variate Gaussian. If the probability of obtaining a dataset is a Gaussian in $N$ dimensions, where $N$ is the number of data points, a GP can describe that dataset. An in-depth introduction to Gaussian processes is provided in Rasmussen and Williams (2005).

GP models parameterise the covariance between data points by means of a kernel function. As a qualitative demonstration, we present the time-series in figure 1: the Kepler light curve of KIC 5809890. This is a relatively active star that rotates once every $\sim$ 30.5 days, with stochastic variability typical of Kepler FGK stars. Clearly, data points in this light curve are correlated. Points close together in time are tightly correlated, and points more widely separated are loosely correlated. A GP models this variation in correlation as a function of the separation between data points; that is, it models the covariance structure rather than the data directly. This lends GPs their flexibility—they can model any time series with a similar covariance structure. In addition, a very simple function can usually capture the covariance structure of a light curve, whereas modelling the time series itself might require much more complexity. Figure 1 demonstrates how a GP model fits the light curve of KIC 5809890.

A range of covariance models, or kernel functions, can describe stellar variability. For example, the most commonly used kernel function, the ‘Squared Exponential’ (SE), defined as follows, could adequately fit the KIC 5809890 light curve:

[TABLE]

Here $A>0$ is the amplitude of covariance, $l$ is the length scale of exponential decay, and $x_{i}-x_{j}$ is the separation between data points. The SE kernel function has the advantage of being very simple, with just two parameters, $A$ and $l$. If $l$ is large, two data points far apart in $x$ will be tightly correlated, and if small they will be loosely correlated. Another property of the SE kernel function is that it produces functions that are infinitely differentiable, making it possible to model a data set and its derivatives simultaneously. However, The SE kernel function does not well describe the covariance in stellar light curves, nor is it useful for the problem of rotation period inference because it does not capture periodic behaviour. Inferring rotation periods thus requires a periodic kernel function. For this reason, we use the ‘Quasi-Periodic’ kernel. Rasmussen and Williams (2005) model QP variability in CO2 concentration on the summit of the Mauna Loa volcano in Hawaii (data from Keeling and Whorf, 2004) using a kernel which is the product of a periodic and a SE kernel: the QP kernel. This kernel is defined as

[TABLE]

It is the product of the SE kernel function, which describes the overall covariance decay, and an exponentiated, squared, sinusoidal kernel function that describes the periodic covariance structure. $P$ can be interpreted as the rotation period of the star, and $\Gamma$ controls the amplitude of the $\sin^{2}$ term. If $\Gamma$ is very large, only points almost exactly one period away are tightly correlated and points that are slightly more or less than one period away are very loosely correlated. If $\Gamma$ is small, points separated by one period are tightly correlated, and points separated by slightly more or less are still highly correlated, although less so. In other words, large values of $\Gamma$ lead to periodic variations with increasingly complex harmonic content. This kernel function allows two data points that are separated in time by one rotation period to be tightly correlated, while also allowing points separated by half a period to be weakly correlated. The additional parameter $\sigma$ captures white noise by adding a term to the diagonal of the covariance matrix. This can be interpreted to represent underestimation of observational uncertainties — if the uncertainties reported on the data are too small, it will be non-zero — or it can capture any remaining “jitter,” or residuals not captured by the effective GP model. We use this QP kernel function (Equation 2) to produce the GP model that fits the Kepler light curve in figure 1.

There are many ways to construct a QP kernel function, involving a range of choices for both the periodic and aperiodic components of the model. The kernel function presented above reproduces the behaviour of stellar light curves, but other model choices can do so as well. We do not attempt to test any other models in this paper, noting only that this kernel function provides an adequate fit to the data. We leave formal comparisons with other kernel function choices to a future publication.

To infer a stellar rotation period $P$ from a light curve, we fit this QP-kernel GP model to the data. As with any model-fitting exercise, the likelihood of the model could be maximised to find the maximum-likelihood value for $P$. In this study, however, we explore the full posterior PDFs using a Markov Chain Monte Carlo (MCMC) procedure. While this approach comes at a computational cost, such posterior exploration importantly provides a justified uncertainty estimate.

Performing inference with GPs is computationally expensive. Is it worth spending the extra computation power to get slightly more accurate, probabilistic rotation periods for a large number of stars? A growing number of astronomers are interested in hierarchical probabilistic modelling, particularly in the stellar and exoplanet communities. These individuals are mindful of carefully inferring probabilistic parameters wherever possible. Probabilistic inference allows one to marginalise over nuisance parameters while formally taking any uncertainty into account when measuring the quantity of interest – the star’s rotation period, in this case. This becomes especially important as the signal-to-noise ratio decreases and the rotation period is no longer tightly constrained, and as the sampling becomes uneven. Throughout the astronomical literature, hierarchical Bayesian inference has been demonstrated to be a robust method of population inference (e.g. Hogg et al., 2010; Foreman-Mackey et al., 2014a; Wolfgang et al., 2015; Rogers, 2015). Probabilistic measurements with full, rigorous propagation of uncertainties are a crucial ingredient for any hierarchical inference where the target is the population level distribution. Inferring the ages of stars in the galaxy via gyrochronology (age-rotation relations) (e.g. Skumanich, 1972; Kawaler, 1989; Barnes, 2003, 2007) is an example of a hierarchical inference problem using stellar rotation. Hierarchical inference allows you to marginalise over individual rotation periods and infer age distribution parameters directly from the light curves. Probabilistic rotation period inference is highly relevant to a large number of galactic archaeology and exoplanet population studies which require rotation period (or age) distributions.

This paper is laid out as follows. The GP method is described in §2. Its performance is demonstrated and compared with literature methods in in §3. In §4 we apply the GP method to real Kepler data, and the results are discussed in §5.

2 GP Rotation Period Inference

In order to recover a stellar rotation period from a light curve using a quasi-periodic Gaussian process (QP-GP), we sample the following posterior PDF:

[TABLE]

where $y$ are the light curve flux data, $\bm{\theta}$ are the hyperparameters of the kernel described in Equation 2, $\mathcal{L}$ is the QP-GP likelihood function, and $p({\bm{\theta}})$ is the prior on the hyperparameters. Sampling this posterior presents several challenges:

• •

The likelihood evaluation is computationally expensive;

• •

The GP model is very flexible, sometimes at the expense of reliable recovery of the period parameter; and

• •

The posterior may often be multimodal.

This Section discusses how we address these challenges through the implementation details of the likelihood (Section 2.1), priors (Section 2.2), and sampling method (Section 2.3).

2.1 Likelihood

The GP likelihood is similar to the simple Gaussian likelihood function where the uncertainties are Gaussian and uncorrelated. The latter can be written

[TABLE]

where $y_{n}$ are the data, $\mu$ is the mean model and $\sigma_{n}$ are the Gaussian uncertainties on the data. The equivalent equation in matrix notation is

[TABLE]

where $\bf{r}$ is the vector of residuals and $\bf{C}$ is the covariance matrix,

[TABLE]

In the case where the uncertainties are uncorrelated, the noise is ‘white’, (an assumption frequently made by astronomers, sometimes justified) and the off-diagonal elements of the covariance matrix are zero. However, in the case where there is evidence for correlated ‘noise’111In our case the ‘noise’ is actually the model! Incidentally, this approach is the reverse of the regression techniques usually employed by astronomers. In most problems in astronomy one tries to infer the parameters that describe the mean model and, if correlated noise is present, to marginalise over that noise. Here, the parameters describing the correlated noise are what we are interested in and the mean model is simply a straight line at $y=0$., as in the case of Kepler light curves, those off-diagonal elements are non-zero. With GP regression, a covariance matrix generated by the kernel function ${\bf K}$ replaces ${\bf C}$ in the above equation (for our purposes, the QP kernel of Equation 2 generates ${\bf K}$).

Evaluating this likelihood for a large number of points can be computationally expensive. For example, evaluating $\mathcal{L}$ for an entire Kepler lightcurve ($\sim$40,000 points) takes about $\sim$5 s— too slow to perform inference on large numbers of light curves222All computational times cited in this section are based on evaluations on a single core of a 2015 Macbook Pro, 3.1 GHz Intel Core i7.. The matrix operations necessary to evaluate the GP likelihood scale as $N\ln(N)^{2}$, where $N$ is the number of data points in the light curve, using the fast matrix solver HODLR (Ambikasaran et al., 2014), implemented in the george (Foreman-Mackey et al., 2014b) python package.

We accelerate the likelihood calculation using two complementary strategies: subsampling the data and splitting the light curve into independent sections. To subsample Kepler data, for example, we randomly select 1/30th of the points in the full light curve (an average of $\sim$1.5 points per day). This decreases the likelihood evaluation time by a factor of about 50, down to about 100 ms. We then split the light curve into equal-sized chunks containing approximately 300 points per section (corresponding to about 200 days), and evaluate the log-likelihood as the sum of the log-likelihoods of the individual sections (all using the same parameters ${\bm{\theta}}$). This reduces computation time because the section-based likelihood evaluation scales as $mn\ln(n)^{2}$, where $n$ is the number of data points per section and $m$ is the number of sections. This method further reduces computation time for a typical light curve (subsampled by a factor of 30) by about a factor of two, to about 50 ms.

2.2 Priors

2.2.1 Non-period Hyperparameters

The flexibility of this GP model allows for posterior multimodality and “over-fitting”–like behavior. For example, if $l$ is small, the non-periodic factor in the covariance kernel may dominate, allowing for a good fit to the data without requiring any periodic covariance structure—even if clear periodic structure is present. Additionally, for large values of $\Gamma$ the GP model becomes extremely flexible and can fit the data without varying the period. Managing this flexibility to reliably retrieve the period parameter requires imposing informative priors on the non-period GP parameters. In particular, we find it necessary to avoid large values of $A$ and $\Gamma$, and small values of $l$; though the exact details of what works well may differ among datasets.

2.2.2 Period

We use an informed period prior, based on the autocorrelation function (ACF) of the light curve. The ACF has proven to be very useful for measuring stellar rotation periods (McQuillan et al., 2012, 2013b, 2014); however, the method has several shortcomings, most notably the inability to deliver uncertainties, but also the necessity of several heuristic choices, such as a timescale on which to smooth the ACF, how to define a peak, whether the first or second peak gets selected, and what constitutes a secure detection. While this paper presents a rotation period inference method that avoids these shortcomings, it seems prudent to still use information available from the ACF. We thus use the ACF to define a prior on period, which can help the posterior sampling converge on the true period.

We do not attempt to decide which single peak in the ACF best represents the true rotation period, but rather we identify several candidate periods and define a weighting scheme in order to create a noncommital, though useful, multimodal prior. While this procedure does not avoid heuristic choices, we soften the potential impact of these choices because we are simply creating a prior for probabilistically justified inference rather attempting to identify a single correct period.

As another innovation beyond what ACF methods in the literature have presented, we also bandpass filter the light curves (using a 5th order Butterworth filter, as implemented in SciPy) before calculating the autocorrelation. This suppresses power on timescales shorter than a chosen minimum period $P_{\rm min}$ and longer than a chosen maximum $P_{\rm max}$, producing a cleaner autocorrelation signal than an unfiltered light curve.

We use the following procedure to construct a prior for rotation period given a light curve:

For each value of $P_{i}$, where $i=\{1,2,4,8,16,32,64,128\}$ d, we apply a bandpass filter to the light curve using $P_{\rm min}=0.1$ d and $P_{\rm max}=P_{i}$. We then calculate the ACF of the filtered light curve out to a maximum lag of $2P_{i}$ and smooth it with a boxcar filter of width $P_{i}/10$. 2. 2.

We identify the time lag corresponding to the first peak of each of these ACFs, as well as the first peak’s trough-to-peak height, creating a set of candidate periods $T_{i}$ and heights $h_{i}$. 3. 3.

We assign a quality metric $Q_{i}$ to each of these candidate periods, as follows. First, we model the ACF as a function of lag time $t$ as a damped oscillator with fixed period $T_{i}$:

[TABLE]

and find the best-fitting parameters $A_{i}$ and $\tau_{i}$ by a non-linear least squares minimization procedure–though enforcing a maximum value of $\tau=20*P_{i}/T_{i}$ to avoid unreasonably large solutions of $\tau$. We then define the following heuristic quality metric:

[TABLE]

where $h_{i}$ is the height of the ACF peak at $T_{i}$, $N_{i}$ is the length of the lag vector in the ACF (directly proportional to the maximum allowed period $P_{i}$), and $R_{i}$ is the sum of squared residuals between the damped oscillator model and the actual ACF data. The idea behind this quality metric is to give a candidate period a higher score if

1. [(a)]
2. (a)

it has many regular sinusoidal peaks, such that the decay time $\tau_{i}$ is long compared to the oscillation period $T_{i}$, 3. (b)

the ACF peak height is high, and 4. (c)

the damped oscillator model is a good fit (in a $\chi^{2}$ sense) to the ACF, with extra bonus for being a good fit over more points (larger $N_{i}$, or longer $P_{i}$). 4. 4.

Given this set of candidate periods $T_{i}$ and quality metrics $Q_{i}$, we finally construct a multimodal prior on the $P$ parameter of the GP model as a weighted mixture of Gaussians:

[TABLE]

That is, in addition to taking the candidate periods themselves as mixture components, we also mix in twice and half each candidate period at a lower level, which compensates for the possibility that the first peak in the ACF may actually represent half or twice the actual rotation period. The period width of 0.2 in log space (corresponding to roughly 20% uncertainty) is again a heuristic choice, balancing a healthy specificity with the desire to not have the results of the inference overly determined by the ACF prior.

Incidentally, while we use the procedure described here to create a prior on $P$ which we use while inferring the parameters of the quasi-periodic GP model, this same procedure may also be used in the service of a rotation-period estimating procedure all on its own, perhaps being even more robust and accurate than the traditional ACF method. We leave exploration of this possibility to future work.

2.3 Sampling

To sample the posterior in a way that is sensitive to potential multimodality, we use the emcee3333https://github.com/dfm/emcee3 MCMC sampler. emcee3 is the successor to the emcee project (Foreman-Mackey et al., 2013) that includes a suite of ensemble MCMC proposals that can be combined to efficiently sample more distributions than the stretch move (Goodman, J. and Weare, 2010) in emcee. For this project, we use a weighted mixture of three proposals. First, we include a proposal based on the kombine package444https://github.com/bfarr/kombine (Farr & Farr, in prep.) where a kernel density estimate (KDE) of the density represented by the complementary ensemble is used as the proposal for the other walkers. The other two proposals are a “Differential Evolution (DE) MCMC” proposal (Ter Braak, 2006; Nelson et al., 2014) and the “snooker” extension of DE (ter Braak and Vrugt, 2008).

We initialise 500 walkers with random samples from the prior and use a weighted mixture of the KDE, DE, and snooker proposals with weights of 0.4, 0.4, and 0.2 respectively. We run 50 steps of the sampler at a time, checking for convergence after each iteration, up to a maximum of 50 iterations. We declare convergence if the total effective chain length is at least 8$\times$ the maximum autocorrelation time. When convergence is achieved, we discard the first two autocorrelation lengths in the chain as a burn-in, and randomly choose 5000 samples as representative of the posterior. This fitting process takes several hours for a typical simulated light curve, though in some cases it can take 12 hours or longer to converge.

3 Performance and Comparison to Literature: Simulated Data

In order to test this new rotation period recovery method, we apply it to a set of simulated light curves and compare to the performance of established literature methods. Section 3.1 describes the simulated data we use; Section 3.2 demonstrates the performance of the QP-GP method; and Section 3.3 compares it to the performance of the Lomb-Scargle periodogram and autocorrelation function methods.

3.1 Simulated light curves

We take our test data set from the Aigrain et al. (2015) ‘hare and hounds’ rotation period recovery experiment. These light curves result from placing dark, circular spots with slowly evolving size on the surface of bright, rotating spheres, ignoring limb-darkening effects. Aigrain et al. (2015) simulated one thousand such light curves to test the ability of participating teams to recover both stellar rotation periods and rotational shear (the amplitude of surface differential rotation). However, in this work, in order to focus on demonstrating reliable period recovery, we select only the 333 light curves without differential rotation, as differential rotation may produce additional scatter in the measured rotation periods.

Each of these light curve simulations uses a real Kepler long-cadence time array: one data point every thirty minutes over a four year duration. 90% of the rotation periods of the simulations come from a log-uniform distribution between 10 and 50 days, and 10% from a log-uniform distribution between 1 and 10 days. Figure 2 shows the distribution of solid-body rotation periods. The simulations also have a range of stellar inclination angles, activity levels, spot lifetimes and more (see Table 1). In order to preserve Kepler noise properties, Aigrain et al. (2015) add real Kepler light curves with no obvious astrophysical variability to the theoretical rotationally modulated light curves. Figure 3 shows an example of a simulated light curve with a period of 20.8days.

3.2 Method Performance

We apply the QP-GP inference method described in Section 2 to each of these 333 simulated light curves. As discussed in Section 2.2.1, reliable inference requires defining a useful set of priors on the non-period hyperparameters. For this simulation dataset, we determined these by first running the method using very broad priors on all the non-period parameters (log-flat between -20 and 20) and then inspecting the distribution of their posteriors for those cases that successfully recovered the true period. We also experimented with constraining the allowed ranges of the parameters after discovering that some regions of parameter space (such as large values of $A$ and $\Gamma$ and small values of $l$) tended to allow fits that ignored the desired periodicity. We list the final priors and bounds this process led us to adopt in Table 2. For the period prior, we tried two different methods: an uninformed (log-flat) prior between 0.5 d and 100 d, and an ACF-informed prior (Section 2.2.2).

Figures 4 and 5 summarise our results compared to the injected ‘true’ stellar rotation periods, for the uninformed and ACF-informed priors on period, respectively. The periods inferred using the GP method show good agreement with the true periods. To assess the performance of the QP-GP and other period recovery methods, we compute Median Absolute Deviations (MADs) of the results, relative to the input periods. We also compute the Median Relative Deviation (MRD), as a percentage and the root-mean-squared (RMS). These metrics are presented for the three different methods tested in this paper in table 3. The informative and uninformative prior versions of the GP method have MRDs of 1.97% and 2.62% respectively. The marginal posterior distributions of the QP kernel hyperparameters, for the example simulated light curve in figure 3, are shown in figure 6.

3.3 Comparison with literature methods

3.3.1 ACF

The ACF method has proven extremely useful for measuring rotation periods. The catalogue of rotation periods of Kepler stars provided in McQuillan et al. (2013a) has been widely used by the community and has provided ground-breaking results for stellar and exoplanetary science. The method also performed well in the Aigrain et al. (2015) recovery experiment, producing a large number of accurate rotation period measurements (see, e.g., their figure 8). Another advantage is its fast implementation speed. However, because the ACF method is non-probabilistic, ACF-estimated rotation period uncertainties are poorly defined — a clear disadvantage. It also requires evenly-spaced data, therefore is not applicable to ground-based light curves. Figure 7 contains an example ACF of the light curve in Figure 3.

We compare our results to those of the Tel Aviv team in Aigrain et al. (2015) who use the ACF method described in McQuillan et al. (2014) and Aigrain et al. (2015). This method involves the calculation of an autocorrelation function, a measure of the self-similarity of the light curve over a range of lags, which is then smoothed with a Gaussian kernel. Periodic variability in the light curve produces a peak in the ACF at the period of the signal and a decaying sequence of peaks at integer multiples of the period. McQuillan et al. (2014) first pre-processed the light curves by removing very long-term variations. They then calculated an ACF for each star and measured the positions of up to the first four peaks in the ACF. They fit a straight line to the periods of these peaks as a function of their integer multiple and adopted the gradient of that line as the rotation period. The uncertainty on the ACF period is the uncertainty on the best-fit slope.

Figure 8 shows ACF-measured versus true rotation periods, with the $2n$ and $\frac{1}{2}n$ harmonic lines as dashed lines. The injected and recovered rotation periods generally agree well: the MRD of their periods, relative to the true rotation periods, is 1.75% — slightly smaller than the results produced by the GP method. However, their results contain a larger number of extreme outliers — they predict periods that are significantly different from the true period in more cases than the GP method. This is revealed by comparing the RMS of the two sets of results since the RMS is more sensitive to outliers. The RMS of the Tel Aviv method results is 4.58 days and the RMS of GP method results is 0.25 days.

3.3.2 LS periodogram

The Aigrain et al. (2015) Tel Aviv team were the only group to report periods for all 333 simulated stars, and are therefore the only team with which we compare results. The other participating groups chose to report only the period measurements in which they were confident and are likely to have omitted large outliers. In order therefore to test our results against another commonly used technique, we implemented a simple LS periodogram method.

We first applied a high-pass filter to the light curves, removing long-term trends that could produce large peaks at long periods in the periodogram. We used the 3rd order Butterworth filter in scipy with a cut-off of 35 days, attenuating signals with periods greater than this threshold. This filter was applied in frequency space, rather than period space, and full attenuation was attained at zero frequency. In other words, the level of attenuation increased smoothly with period, for periods above 35 days but no periods were fully suppressed. Experimentally, we found the 35 day cut off removed the largest number of long period outliers while preserving shorter period signals. We also tried 50, 45, 40 and 30 day cut-offs but the 35 day cut-off value minimised the RMS and MAD of the LS periodogram results. For each simulated light curve, we computed a LS periodogram555LS periodograms were calculated using the gatspy Python module: https://github.com/astroML/gatspy/tree/master/gatspy/periodic. over a grid of 10,000 periods, evenly spaced in frequency, between 1 and 100 days. We adopted the period of the highest peak in the periodogram as the measured rotation period. The uncertainties on the rotation periods were calculated using the following equation for the standard deviation of the frequency (Horne and Baliunas, 1986; Kovacs, 1981):

[TABLE]

where $A$ is the amplitude of the signal of highest power, $\sigma_{N}$ is the variance of the time series, with the signal of highest power removed, $N$ is the number of observations and $T$ is the timespan of the data. These formal uncertainties are only valid in the case that the noise is white, the data are evenly sampled and there is only one signal present. Since there are multiple signals present in these light curves, this formal uncertainty is an underestimate of the true uncertainty.

Figure 9 shows the resulting recovered rotation periods as a function of true period. The MRD of the periodogram-recovered periods is 2.82%, slightly worse than the ACF and GP methods (see table 3 for a side-by-side comparison with the other methods). The filter suppresses signals with periods greater than 35 days but doesn’t eliminate them altogether, thus allowing high amplitude frequencies to be detectable. For this reason, some erroneous long periods are recovered by the periodogram method. These results still show a marked improvement in comparison to periodograms calculated from unfiltered data. Simply adopting the highest peak of a periodogram calculated from the unfiltered simulated data resulted in an RMS of 18.98 and relative MAD of 4.70.

4 Real Kepler data

In order to test our rotation period inference method on real data, we apply it to a set of 275 Kepler Object of Interest (KOI) host stars, that had rotation periods previously measured by the ACF method by McQuillan et al. (2013a). We use the pipeline-corrected flux (pdcsap_flux column in the Kepler light curve table), median-normalized and unit-subtracted, and mask out all known transiting planet candidate signals. As with the simulated light curves, we randomly subsample each light curve by a factor of 30 and split it into segments of about 300 points for the purposes of evaluating the likelihood. We also follow the same MCMC fitting procedure as with the simulated data, using the ACF-based prior as before.

Initially, we also use the same priors on the hyperparameters for the KOIs as for the simulated light curves (Table 2). However, we found that $\ln l$ and $\ln\Gamma$ tended toward slightly different values than the simulations. We also found that the allowed hyperparameter range that we used for the simulations was too large for the KOI population, as maybe $\sim$15% of the fits tended toward corners in the hyperparameter space, resulting in poor period measurements. As a result, after this initial test, we subsequently adjusted the priors and re-fit all the KOIs. The final priors and parameter ranges that we used are in Table 4. The posterior samples of the model parameters for the Kepler objects of interest are available online: https://zenodo.org/record/292340#.WKWpiBIrJE4.

Figure 10 compares the periods inferred with the GP method to the ACF-based periods from McQuillan et al. (2013a) for the 275 overlapping KOIs. This comparison shows generally very good agreement, demonstrating that this method works not only on simulated data, but also on real data—with the caveat that for any particular data set, some care is needed regarding the setting the priors and ranges for the GP hyperparameters. The new GP rotation periods calculated for 1102 KOIs are provided in an online table accompanying this paper. Samples of the rotation period posterior PDFs are available online: https://doi.org/10.5281/zenodo.804440

5 Discussion

The majority of long period outliers seen in the periodogram comparison figure (figure 9) are produced when the rotation period signal is highly non-sinusoidal. This can be understood intuitively as the LS periodogram fits a series of sine waves to the data; if the data are not sinusoidal, a sine wave will not produce a good fit. It is difficult to assess the reason for the ACF incorrect recoveries because we do not have the individual ACFs computed by the Aigrain et al. (2015) Tel Aviv team, however it seems that the ACF method is more sensitive than the LS periodogram method to noise. This particular implementation of the LS periodogram method is capable of detecting periodic signals at lower amplitudes that the ACF. We do not suggest that the LS periodogram method is better than the ACF method in general. It was made clear in Aigrain et al. (2015) that the performance of each of these rotation period detection algorithms depends sensitively on the implementation technique. The exact steps performed in the pre-processing of the light curves and the post-processing of the periodogram or ACF seem to affect the results strongly. This is also true for the GP method. Using a different MCMC sampler, different priors, different subsampling, or a different kernel would all alter the results. Clearly, none of these rotation period measurement methods is perfect and they all require a lot of background engineering and tuning. This is simply due to inherent difficulty of inferring periods, or any parameter that produces a highly multimodal probability distribution, from noisy data.

The QP-GP inference method we test in this paper produces slightly more accurate rotation periods than the ACF or periodogram methods. Additionally, because it explores the posterior PDF with MCMC, the GP method produces probabilistically justified uncertainties. Unfortunately, these uncertainties still appear to be underestimated. Of the rotation periods recovered from the simulated light curves using the ACF-informed prior, only 25% of the measured periods lie within 1$\sigma$ of the true period, 50% lie within 2$\sigma$, and 66% within 3$\sigma$. The largest outlier is 114$\sigma$ away from the true value. The median uncertainty is underestimated by around a factor of two. These numbers are similar for the uninformative prior results. We attribute these underestimated uncertainties to the model choice. Although more appropriate than a sinusoid, the QP-GP is still only an effective model. A perfect physical model of the star would produce more accurate uncertainties.

The QP kernel function represents a simple effective model of a stellar light curve. It can describe a wide range of quasi-periodic behavior and is relatively simple, with only a few hyperparameters. Nevertheless, it is still a somewhat arbitrary choice. Another valid choice would be a squared cosine function multiplied by a squared exponential, used by Brewer and Stello (2009) to model asteroseismic pulsations,

[TABLE]

This function produces a positive semi-definite matrix and has the $P$ parameter of interest, but differs qualitatively from the QP function by allowing negative covariances. Is it realistic to allow negative covariances? In practice, the ACFs of Kepler light curves often go negative. However, many stars have two active regions on opposite hemispheres that produce two brightness dips per rotation, and it may be difficult to model such stars with a kernel that forces anti-correlation of points $\frac{1}{2}$ a period apart. It would be worthwhile to test this assumption and this alternative kernel function (and others) in the future.

Not all Kepler light curves show evidence of stellar rotation. In some cases a star may have few or no active regions, be rotating pole-on, or be rotating so slowly that the Kepler data detrending pipeline removes any signal. In other cases there may be another source of variability present in the light curve, generating a false period detection. These sources may be physical: e.g. binary star interactions, intra-pixel contamination from other astrophysical objects, pulsating variable stars, asteroseismic oscillations in giants and even stellar activity cycles. Identifying many of these astrophysical false positives falls outside the scope of this GP method (e.g. applying colour cuts to avoid giant contamination); however, for some, such as variable stars, they may have distinctive hyperparameters (e.g. long coherence timescales) that identify them. Testing this may be an interesting follow-up study. As well as astrophysical contamination, instrumental sources may contribute to contaminating variability, e.g. temperature variations or pointing shifts of the Kepler spacecraft. These are unlikely to be periodic and, again, may produce unusual combinations of hyperparameters. We also hope to test this in the future.

In addition, we are continuing to develop several other aspects of this GP method:

• •

To use an updated method for calculating the GP likelihood and drastically improve computational efficiency. celerite, a new method for rapid computation of the ‘solve’ and determinant calculations of covariance matrices necessary for GP inference has recently been developed (Foreman-Mackey et al., 2017). celerite scales with the number of data points, $\mathcal{O}(N)$ rather than the $N\log(N)^{2}$ scaling of the HODLR algorithm used in this analysis. It requires a slightly different functional form for the kernel as the increased speed is obtained by using mathematical properties of exponential functions. Kernels must be written as the sums of exponentials in order to be implemented in celerite, however a close approximation to the kernel function used here is demonstrated to recover rotation periods in Foreman-Mackey et al. (2017).

• •

To investigate more physically motivated priors. In this paper we have only explored the physical interpretation of one parameter in the kernel function ($P$), which we related to stellar rotation periods, but others also warrant physically motivated priors. For example, the harmonic complexity parameter, $\Gamma$, defines the typical number of inflection points per function period; suitable priors may be informed by the typical properties of observed light curves. For the evolution time-scale, $l$, it is sensible to enforce $l>P$, such that functions that have clear periodicity. Intuitively, if $l<P$, functions may evolve significantly over time scales shorter than one period, so a function $f(t+P)$ would look significantly different to $f(t)$, in which case any claims of periodicity become dubious. More stringent criteria may be suggested by stellar astrophysical considerations, e.g. knowledge of typical active region evolution time-scales or persistence lifetimes (e.g. $l\sim 2P$). Finally, the amplitude hyper-parameter $A$ may be related to active region coverage. Unfortunately, because these hyper-parameters control the covariance structure of functions, their interpretation is not always straightforward; interpretation is further complicated by significant degeneracies between the hyperparameters. Therefore we defer such considerations, along with strategies for efficient optimization or marginalization of hyperparameters, to a separate paper (currently in prep.).

• •

To build in a noise model for Kepler data. Because it is a generative model of the data, the GP method models the rotation period at the same time as systematic noise features. One can then marginalise over the parameters of the noise model. This approach would be extremely advantageous for Kepler data since long-term trends are often removed by the Kepler detrending pipeline. Marginalising over the noise model at the same time as inferring the parameters of interest will insure that the periodic signal is preserved.

6 Conclusion

We have attempted to recover the rotation periods of 333 simulated Kepler-like light curves for solid-body rotators (Aigrain et al., 2015) using three different methods: a Gaussian process method, an autocorrelation function method, and a Lomb-Scargle periodogram method. We demonstrate that the GP method produces the most accurate rotation periods of the three techniques. Comparing our results with the results of the Aigrain et al. (2015) Tel Aviv team who implement the ACF method of McQuillan et al. (2013a), we find their results to be slightly more precise in a median sense than the GP method results, but with a larger number of significant outliers. Additionally, we measure the rotation periods of 275 Kepler objects of interest using the GP method, and find that these results compare well to those measured previously by McQuillan et al. (2013a). The good agreement between these two sets of results demonstrates that the GP method works well on real Kepler data. Samples of the full set of 1102 KOI rotation period posterior PDFs can be found at https://doi.org/10.5281/zenodo.804440.

Unlike the ACF and LS periodogram methods, the GP method provides posterior PDF samples which can be used to estimate rotation periods uncertainties. In general it provides more accurate rotation periods than either method and produces significantly fewer outliers. The GP method can also be applied to non-uniformly sampled light curves, unlike the ACF method.

Although an improvement on competing methods, the uncertainties are still underestimated by the GP method in many cases because it is an approximate model. A quasi-periodic Gaussian process is clearly a good model, as demonstrated here, however it is still only an effective model, not an accurate physical one. It captures the posterior PDF of the periodic component of the covariance matrix of a time series, not the actual rotation period of a physical star. The marginalised posterior PDF does not therefore exactly reflect the probability of a rotation period, given the data. Similarly, the 16th and 84th percentile ranges calculated from the posterior PDF of the QP-GP model, are not exactly the uncertainty of a physical rotation period. Although underestimated uncertainties are a drawback of this method, they are at least well motivated and well defined which is not the case for uncertainties calculated from the ACF and LS periodogram methods.

Although GP inference is slow as it requires a ‘solve’ and a matrix determinant calculation at every step of an MCMC, a new GP solver, celerite (Foreman-Mackey et al., 2017) shows promise for increasing the computional speed of GP rotation period inference by several orders of magnitude. We intend to integrate celerite into our code in the near future. How much computational time is too much is a question that depends on the available resources of each individual and their scientific priorities. GP period inference is trivially parallelisable, so if one has access to thousands of computer cores, one can infer rotation periods of thousands of stars in a few hours (perhaps just a few minutes with celerite). For single stars GP rotation inference takes a reasonable amount of time and has already been used on exoplanet hosts (e.g. Haywood et al., 2014; Vanderburg et al., 2015).

The code used in this project is available at https://github.com/RuthAngus/GProtation/.

This research was funded by the Simons Foundation and the Leverhulme Trust. TDM is supported by NASA grant NNX14AE11G, and acknowledges the hospitality of the Institute for Advanced Study, where part of this work was completed. VR thanks Merton College and the National Research Foundation of South Africa for financial support. We would like to thank the anonymous referee who’s comments significantly improved this paper. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate.