DETERMINING THE MASS OF KEPLER-78b WITH NONPARAMETRIC GAUSSIAN PROCESS ESTIMATION

The Kepler telescope obtained photometry of approximately 150,000 objects in the Kepler field for its four year lifetime from 2009 to 2013. Photometry of Kepler-78 was gathered at a 30 minutes cadence. In this study, we train a GP on the photometric light curve of Kepler-78 to determine the evolution and rotation timescales of the stellar activity, hyperparameters of the GP regression kernel. However, as shown in Equation (2), a matrix inversion is required to calculate the log posterior likelihood of a GP kernel with a given set of hyperparameters, a computationally intensive process with compute time proportional to N³. Therefore, it was necessary to rebin the Kepler photometry to one point every 5 hr (averaging every ten points together) in order to find the best-fit hyperparameters of a full quarter of photometry with our computing resources. While this does marginalize over the planetary signal, we are only using the photometry to estimate the underlying components of the stellar activity which vary on timescales much longer than 5 hr, and the photometric effect of the planet is small compared to that of the stellar activity. This photometry, taken from the 16th quarter of Kepler observations, is not concurrent with the RV measurements; however, since it was taken only ≈55 days before the RV measurements, we can assume that the variation seen in the photometry is at least related to the RV variation, and would exhibit similar structure in the temporal dimension. Figure 1 illustrates the binned photometric measurements as well as the GP regression with the best-fit kernel hyperparameters to the photometry.

Zoom In Zoom Out Reset image size

The RV observations were obtained shortly after the malfunctioning of Kepler's reaction wheels in summer 2013, which prevented the high level of photometric precision achieved in the first era of the spacecraft's lifetime and required that the telescope point at new targets. H13 observed Kepler-78 with HIRES on the Keck I telescope on Mauna Kea. They obtained eight nights of data during each of which approximately nine or more measurements were obtained: three measurements were taken at successive thirty minute intervals at three separate occasions per night. In addition, five individual measurements were taken on separate nights after the rest of the data were collected to provide a longer baseline over which to evaluate the stellar activity. The data set consists of 84 RVs over 45 nights. The P13 team observed Kepler-78 with HARPS-N on the Telescopio Nazionale Galilei in the Canary Islands. They were able to obtain 109 measurements over 97 nights. Most of their observations were made in a consecutive six days period. On subsequent nights only one or two measurements were made per night, giving a total time baseline of approximately three months. Figure 2 illustrates the HIRES and HARPS-N measurements with measurement errors shown, with GP regressions overplotted. Figure 3 shows the phased RVs with the best-fit GP regression subtracted.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

GP is a nonparametric method to describe a data set by evaluating correlations between n data points through a covariance kernel. This kernel describes the relationship of each point in the data set to each other point, and can be expressed as an n × n matrix (subsequently referred to as the covariance matrix). The kernel is a function of hyperparameters. More complicated kernels can have more hyperparameters that characterize different qualities of the correlations in the data, such as various periods, characteristic amplitudes, and length scales, etc.

GP regression is widely used in the field of machine learning (Neal 1997; Herbrich et al. 2003; Quiñonero-Candela & Rasmussen 2005; Wang et al. 2008). Gibson et al. (2012) introduced the technique to the field of exoplanets through analysis of transmission spectroscopy to model correlated noise in the instrumental systematics of HST/NICMOS, as described in Section 1. Concurrently with this work, Haywood et al. (2014) have demonstrated the technique of GP modeling of RV and photometric signals for the CoRoT-7 planetary system, first modeling the photometry with a GP and then using the photometric GP hyperparameters to train the initial RV GP hyperparameters. Haywood et al. (2014) demonstrated that in the case of CoRoT-7b, parametric spot models such as those used by H13 and P13 gave incorrect masses and uncertainties. Thus, it is important to test many time series techniques, and further explore the novel application of GPs in Doppler analysis. We expand upon the GP training technique used in Haywood et al. (2014) here.

Finding the best GP regression requires choosing a kernel and initial hyperparameters, evaluating the likelihood of those hyperparameter values, and then iterating through parameter space until the most likely values are found. The squared exponential kernel, for example, defines a covariance matrix through an operator,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{ij}}=k({t}_{i},{t}_{j})={h}^{2}\mathrm{exp}\left[-{\left(\displaystyle \frac{{t}_{i}-{t}_{j}}{\lambda }\right)}^{2}\right],\end{eqnarray} \tag{ 1 }$$

where h is the covariance amplitude, and λ the covariance length scale. The amplitude observed is described by h, while λ is a characteristic timescale over which the data is going to be correlated.

We discuss other GP kernels and the inferred physical meaning of their hyperparameters in Table 2.

Table 2.  Gaussian Process Kernel Options

Name	Mathematical Expression	Hyperparametersa	Comments
Squared exponential	${h}^{2}\mathrm{exp}\left[-{\left(\frac{{t}_{i}-{t}_{j}}{\lambda }\right)}^{2}\right]$	h, λ	h amplitude of covariance function,
			λ a characteristic timescale
Periodic	${h}^{2}\mathrm{exp}\left[-\frac{{\mathrm{sin}}^{2}[\pi ({t}_{i}-{t}_{j})/\theta ]}{2{w}^{2}}\right]$	h, θ, w	θ equivalent to Prot,
			w represents coherence scale, similar to λ expressed
			as a fraction of θ dependent on recurrent features
Quasi-periodic	${h}^{2}\mathrm{exp}\left[-\frac{{\mathrm{sin}}^{2}[\pi ({t}_{i}-{t}_{j})/\theta ]}{2{w}^{2}}-{\left(\frac{{t}_{i}-{t}_{j}}{\lambda }\right)}^{2}\right]$	h, θ, w, λ	w coherence scale tied to periodic variation,
			while characteristic timescale λ tied to aperiodic variation.

Notes. The name of kernel functions and hyperparameters in this table have been taken from Rasmussen & Williams (2006).

^aEach kernel ${{\boldsymbol{\Sigma }}}_{{ij}}$ can be modified to include an additional hyperparameter, a white noise term ${\sigma }^{2}$ by adding one in quadrature: ${{\boldsymbol{\Sigma }}}_{{ij}}={{\boldsymbol{\Sigma }}}_{{ij}}$+ ${\sigma }^{2}{{\boldsymbol{I}}}_{i}$.

Download table as:  ASCIITypeset image

Table 3.  Planetary Properties

Property	Value
Name	Kepler-78b
Radius, Rpl	1.20 ± 0.09 ${R}_{\oplus }$
Orbital Period ${P}_{\mathrm{orb}}$	0.35500744 ± 0.00000006 days
Doppler Amplitude, K	1.86 ± 0.25 m s−1
Mass, Mpl	${1.87}_{-0.26}^{+0.27}\;{M}_{\oplus }$
Density, ${\rho }_{\mathrm{pl}}$	${6.0}_{-1.4}^{+1.9}$ g cm−3
Iron Fraction	32% ± 26%

Note. The name, radius, and orbital period values in this table have been taken from Howard et al. (2013). All other values have been calculated for this work.

Download table as:  ASCIITypeset image

The logarithm of the posterior likelihood of the GP regression is calculated as

$$\begin{eqnarray}&&\mathrm{log}[{\mathcal{L}}({\boldsymbol{r}})]=-\displaystyle \frac{1}{2}{{\boldsymbol{r}}}^{{\rm{T}}}{{\boldsymbol{\Sigma }}}^{-1}{\boldsymbol{r}}-\displaystyle \frac{1}{2}\mathrm{log}| {\boldsymbol{\Sigma }}| -\displaystyle \frac{n}{2}\mathrm{log}(2\pi ),\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{r}}$ is the vector of residuals after removal of the (optional) mean function, ${\boldsymbol{\Sigma }}$ is the covariance matrix, and n the number of data points. A prior term, ${{\mathcal{L}}}_{\mathrm{prior}}$, can be added to the likelihood to account for any priors placed on the hyperparameters. For example, we apply the Gaussian prior

$$\begin{eqnarray}&&{{\mathcal{L}}}_{\mathrm{prior}}=\mathrm{exp}\left({-\frac{1}{2}\left[{\left(\frac{{\theta }_{\mathrm{true}}-\theta }{{\sigma }_{\theta }}\right)}^{2}\right]}^{}\right),\end{eqnarray} \tag{ 3 }$$

to restrict the hyperparameter θ for the RV data, as the period of the correlated noise signal is equivalent to the stellar rotation period P_rot in both the photometric and RV data sets. We restrict the period hyperparameter θ in the RV regressions within a Gaussian of width ${\sigma }_{\theta }$ determined by the posterior distribution of the period hyperparameter found for our photometric GP regression. Prior knowledge of this hyperparameter helps to ensure convergence of the other hyperparameters of the RV regression, as the RV data is not as well sampled as the photometry.

This likelihood calculation can be used to identify the best fit GP hyperparameters. For more details on how the covariance kernel and parameter boundary conditions were chosen, refer to Section 4.1. For a more complete description of GP regression and posterior likelihood evaluation, see Rasmussen & Williams (2006).

The presence of magnetic features on the stellar surface induces variation in a star's RV that can mimic a planet's orbital RV signal. Starspots are regions of high magnetic fields on the surface of a star that appear relatively cooler and darker than the surrounding photosphere. As the starspots move across the disk of a star (due to stellar rotation), the flux balance between the redshifted and blueshifted halves of the star is broken, producing a shift in the centroid of the spectral lines. This shift corresponds to an apparent RV signal and changes as the spots move across the stellar disk (Dumusque et al. 2011). Starspots, as well as networks of strongly magnetized flux tubes known as faculae, inhibit the convective processes taking place on the stellar surface, thus reducing the net blueshift produced by the up flow of hot, bright granules. This has been shown to be the dominant RV effect on the Sun (Meunier et al. 2010) and on other Sun-like stars (Haywood et al. 2014). These activity-induced RV signals are modulated by the stellar rotation and the surface features evolve with time, resulting in a quasi-periodic, RV signal. The lifetimes of surface features are poorly constrained, making this effect difficult to model parametrically. The nonparametric nature of the GP regression provides an ideal framework for studying stellar surface features in RV. We can choose a GP kernel that reflects the frequency structure of the stellar activity. The resultant GP model is flexible enough to account for the evolution of magnetic features on the stellar surface while keeping a statistical "memory" of the stellar activity patterns.

We fit the photometric data using GP regressions with three different kernels. We test a squared exponential, periodic and quasi-periodic kernel (Table 2). We test kernels with and without an additional white noise term added in quadrature to the likelihood (Equation (2)). The best-fit kernel hyperparameters are found via Markov Chain Monte Carlo analysis powered by the emcee Python package (Foreman-Mackey et al. 2013), described in more detail in the following section. We select the quasi-periodic model with a white noise term based on our assumption that the predominant signal in both data sets is tied to the quasi-periodic variation of the stellar surface, and confirm that the quasi-periodic model is appropriate through a visual check and a reduced ${\chi }^{2}$ analysis. We have plotted the data and GP regression with the best-fit kernel hyperparameters in Figure 1. We ensure that these hyperparameters are robust by performing a GP regression to the photometric data after less stringent binning, and recovering the same best-fit kernel hyperparameter values within errors.

We fit our photometric and RV data sets with the same kernel operator. We took the hyperparameter values from the photometric GP regression and used them as initial fit values for the RV GP kernel hyperparameters, with the exception of the hyperparameters h and σ, as the variation in the photometry and RV related to these quantities is not in equivalent units. We also place a Gaussian prior on θ based on its posterior distribution determined from the photometric GP regression. We then apply a uniform offset to each set of RV measurements to remove RV zero point differences.

The likelihood of the GP regression, given certain kernel hyperparameters and a Doppler amplitude value, determines the quality of the model for those given parameters. This likelihood can be calculated as given in Equation (2), where the residuals ${\boldsymbol{r}}$ can be calculated as

$$\begin{eqnarray}&&{\boldsymbol{r}}={\boldsymbol{v}}-K\mathrm{sin}\left(\displaystyle \frac{2\pi ({\boldsymbol{t}}-{t}_{c})}{{P}_{\mathrm{orb}}}\right),\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{t}}$ is the vector of times of all measurements, ${\boldsymbol{v}}$ the vector of RV measurements, t_c a time of transit, ${P}_{\mathrm{orb}}$ the orbital period, and K the planetary Doppler amplitude. Sanchis-Ojeda et al. (2013) measured t_c to 10⁻⁵ of a day and and ${P}_{\mathrm{orb}}$ to 10⁻⁷ of a day precision, so we can safely treat them as constants for this analysis.

The best-fit GP kernel hyperparameters and the Keplerian Doppler amplitude are found via MCMC exploration of parameter space. We simultaneously fit the Keplerian planetary signal and a GP regression to each RV data set with the Doppler amplitude and kernel hyperparameters as free parameters in the MCMC chain. The emcee package contains an Affine-invariant Monte Carlo Markov Chain Ensemble sampler, which evaluates the likelihood of the GP kernel hyperparameters to the measurement residuals after a Keplerian planetary signal has been subtracted (Foreman-Mackey et al. 2013). This is done for a plethora of steps in free parameter space. We draw the best-fit GP kernel hyperparameters and planetary Doppler amplitude as well as their errors from the posterior distributions generated through this MCMC exploration of parameter space, with 1σ error corresponding to 68% confidence intervals in the MCMC posterior distributions. We plot the data and best-fit GP regression + Keplerian models in Figure 2.

After the MCMC posterior distributions have been created, we calculate a Gelman–Rubin statistic for each parameter to ensure that the parameter chains have converged. Convergence is deemed adequate when the G–R statistic <1.01 (Gelman & Rubin 1992). In addition, we track the average and median log likelihood at each step in the chain. At first, the chains move to higher and higher likelihood parameter space. Once the best-fit solution is found, the chains begin to take more and more unlikely steps, exploring parameter space more completely, and pushing the average likelihood below the median value. The step at which the average log likelihood begins dipping below the median log likelihood indicates the transition from burn-in period to the exploration of parameter space around the maximum likelihood solution, and we remove the steps taken during burn-in from our analysis (e.g., Knutson et al. 2008). To check that the resultant fit is a good description of the data, we compute a ${\chi }_{\mathrm{red}}^{2}$ for the best-fit model of each kernel variant. Finding ${\chi }_{\mathrm{red}}^{2}\approx 1$ confirms that the quasi-periodic GP model and errors describe the data well.

We allow as many of the MCMC parameters as possible to be minimally constrained. We find that whenever the periodic or quasi-periodic kernel is used, and the photometric θ value is not applied as a prior, the RV noise model is able to recover periodicity consistent with the stellar rotation period. This supports our assumption that the stellar surface features are the predominant signal seen in the RV measurements. Applying a Gaussian prior to this value allows us to better constrain all other MCMC parameters. Physical boundary conditions were also placed on some hyperparameters to prevent the MCMC chain from moving into unphysical parameter space (see Table 3 for our adopted planetary parameters and Table 4 for results from all models).

Table 4.  Adopted Quasi-periodic RV Model Parameters

Name	Prior	Value
HIRES Amplitude h	$h\gt 0$	${11.6}_{-2.5}^{+3.7}$ m s−1
HARPS-N Amplitude h2	${h}_{2}\gt 0$	${5.6}_{-1.3}^{+2.0}$ m s−1
Amplitude Ratio a	⋯	${2.0}_{-0.5}^{+0.8}$
Period θ	$\mathrm{exp}\left(-\frac{{(13.26-\theta )}^{2}}{2{(0.12)}^{2}}\right)$	13.12 ${}_{-0.12}^{+0.14}$ days
Coherence Scale w	⋯	${0.28}_{-0.04}^{+0.05}$
Characteristic Timescale λ	$\lambda \gt 0$	26.1 ${}_{-11.0}^{+19.8}$ days
HIRES White Noise ${\sigma }_{\mathrm{jitter}}$	$\sigma \gt 0$	${2.1}_{-0.3}^{+0.3}$ m s−1
HARPS-N White Noise ${\sigma }_{\mathrm{jitter}}$	$\sigma \gt 0$	${1.1}_{-0.5}^{+0.4}$ m s−1
Doppler Amplitude, K	⋯	1.86 ± 0.25 m s−1

Download table as:  ASCIITypeset image

We find a planetary mass within 1σ of our adopted mass result for all models tested. We make the assumption that the predominant signal observed in the photometric and both RV data sets comes from the star, and thus will modulate at the stellar rotation period. We find that when we model the photometric data with a periodic or quasi-periodic GP kernel, we recover a period hyperparameter consistent with previous estimates of the stellar rotation period, supporting our assumption about the predominant photometric signal. Similarly, when we model the RV data sets with a quasi-periodic GP with no priors, we recover a period hyperparameter consistent with the stellar rotation period or its alias. Thus, we conclude that the quasi-periodic GP regression is both well motivated and effective for describing the predominant signal in our data sets.

We calculated Bayesian Information Criterion (BIC) and Akaike Information Criterion (AIC) values for all of our RV models. The BIC and AIC are criteria for model selection among a finite set of models, where the lowest values are preferred. The two criteria are defined as $\mathrm{BIC}=\mathrm{ln}({\mathcal{L}}+k\mathrm{ln}(n)$ and $\mathrm{AIC}=2k-2\mathrm{ln}({\mathcal{L}})$, where ${\mathcal{L}}$ is the maximized value of the likelihood function, k the number of free parameters estimated, and n the number of data points considered (Akaike 1974, Liddle 2007). We test the three kernels of interest both with shared hyperparameters as well as with fully independent hyperparameters in order to ensure that the variation we observe is consistent between data sets and thus astrophysical, rather than local, in origin. Furthermore, we test our squared exponential kernel models with and without a white noise term for each data set. We report these BIC and AIC values in Table 5.

We find that the quasi-periodic model we choose has AIC and BIC values with ΔBIC $\geqslant -0.5$ and ΔAIC $\geqslant 8$ between it and all other models tested. A ΔBIC value of less than 2 suggests two models are indistinguishably likely, whereas a ΔBIC greater than 10 suggests strong evidence for the model with a lower BIC value. Interestingly, the BIC for the quasiperiodic and squared exponential kernels are almost identical, while there is significantly less evidence for the periodic model (BIC of periodic model 18 points higher than SE and quasi-periodic models). Furthermore, the lowest BIC of all was found for the squared exponential kernel with two independent GP regressions, with a ΔBIC = −0.5 as compared to the chosen model. However, the AIC of the quasi-periodic, coupled GP model is significantly lower than that of either of the other models, with a ${\rm{\Delta }}\mathrm{AIC}$ = 8 between it and the next lowest AIC value, that for the quasi-periodic, uncoupled model, and ΔAIC = 14 for the coupled model with a squared exponential kernel. The relative likelihood of one model to another is given by ${e}^{-{\rm{\Delta }}\mathrm{AIC}/2}$, indicating that the uncoupled quasi-periodic model is 1% as likely as the coupled model, and the coupled squared exponential kernel model is 0.09% as likely. The independent GP, squared exponential kernel had a ΔAIC = 19 when compared to the chosen model. We also tested removing the white noise parameter from the squared exponential model, finding that doing so raised the BIC by over 20 points and the AIC by over 15. Thus, we determined that the white noise parameter was justified for all kernels. Strictly speaking, these information criteria require the assumption that the noise in the data is independent and identically distributed, which we know is not the case. However, since we use the same data for all BIC and AIC comparisons, they provide a valuable comparison between the different GP kernels. We report all relevant AIC and BIC values in Table 5.

If we are modeling the stellar activity in both data sets, we would expect the hyperparameters for our simultaneous models to be the same, with the exception of the non-temporal amplitude hyperparameter h and the white noise term σ, which would be different but related between the models due to the different appearance of the disk of the star in different wavelength regimes. We explored the relation of the GP kernel amplitude hyperparameters, whose relationship we evaluated by calculating a, the ratio of the HARPS-N amplitude to the HIRES amplitude. The convergence of a to a best-fit value suggests that there is a steady relation in the ranges of the RV signals observed by HIRES and HARPS-N. Additionally, the value of a is broadly consistent with the results of Desort et al. (2007), who illustrate that the expected radial velocity contribution from starspots on a K dwarf is larger at shorter wavelengths, resulting in a larger signal in HIRES than in HARPS-N, as HARPS-N has significant wavelength coverage in the 600–700 nm regime that HIRES does not. We show the distribution of a along with the other parameters of our MCMC calculation in Figure 4. Any variation due explicitly to the correlated stellar activity should thus cause our models to have a shape (period and phase) that is consistent between the two data sets. We allow the reader to confirm this visually in Figure 2.

Zoom In Zoom Out Reset image size

Distributions for each parameter relative to all other parameters are shown in Figure 4 for the quasi-periodic GP + Keplerian models. We see no clear correlation between Doppler amplitude and any other parameter tested, indicating that the Doppler amplitude observed is a real signal and not a systematic error introduced during our parameter fitting.

We place a Jeffreys prior on the period hyperparameter θ to weight shorter periods more heavily when fitting the quasi-periodic and periodic GP regression kernel hyperparameters to the photometric data (Haywood et al. 2014). When we do this, we are able to recover the rotation period of the star as the period hyperparameter. We then place a Gaussian prior on this hyperparameter in the RV GP + Keplerian model constraining it to the value and errors measured in the photometric data. We also note that although the coherence scale parameter w measured in the photometric data is equal to that measured in the RV data within errors, we do not place a similar Gaussian prior on this parameter, because we observe that it is correlated with both the amplitude and characteristic timescale hyperparameters, which were not necessarily related between the photometric and RV GP regression hyperparameter fits. This correlation between h, h₂, w and λ visible in Figure 4 is hard to interpret physically. We speculate that the positive correlation between w and the h parameters may arise from a connection to starspot size: as starspots grow larger, h grows, any white noise present will become relatively less important, and thus the function will become smoother overall, increasing the w parameter. Similarly, the characteristic timescale λ may also increase at large amplitudes for the same reason—as the uncorrelated noise at small timescales becomes less important, the lengthscale might be weighted more heavily toward larger values. Larger spots also persist longer on the stellar surface, as spots decay on timescales proportional to their size (Bumba 1963). However, the fact that w was consistent for the photometric and RV data sets whereas λ was smaller for the photometric data sets despite their large differences in amplitudes may suggest otherwise. In addition, the amplitude parameters h and h₂ are not directly related to the range of the RV shift observed, and thus difficult to interpret physically. Further tests of these parameters are necessary in order to fully characterize their relationships, such as placing new priors on the coherence scale and characteristic timescale parameters to ensure that they are not biased by noise.