RADIAL VELOCITY OBSERVATIONS AND LIGHT CURVE NOISE MODELING CONFIRM THAT KEPLER-91b IS A GIANT PLANET ORBITING A GIANT STAR^*

The first discovered exoplanets orbiting non-degenerate stars were Jupiter-sized (Campbell et al. 1988) and most orbited just a few stellar radii from their host star (Mayor & Queloz 1995; Marcy & Butler 1996). This was a milestone moment that demonstrated planetary systems need not resemble our own. While data from the Kepler spacecraft have revealed that planetary systems come in many flavors (e.g., Lissauer et al. 2011; Carter et al. 2012; Barclay et al. 2013), hot Jupiters remain a key area of interest because the large sizes of these planets and their short orbital periods yield the highest signal-to-noise (S/N) light curve data with which otherwise undetectable effects can be observed. Examples of this include the detection of in-homogenous clouds (Demory et al. 2013) and the determination of planet masses from Doppler boosting (Shporer et al. 2011; Barclay et al. 2012).

Additionally, the formation of hot Jupiters is still highly debated, with competing theories including migration through the protoplanetary disk (Lin et al. 1996) and high eccentricity migration triggered by dynamical events such as planet–planet scattering (Nagasawa et al. 2008). The confirmation and characterization of new systems, in particular for evolved stars for which only a handful of hot Jupiters are known, is important to shed light on the origin of these planets.

Kepler-91 was designated KIC 8219268 in the Kepler Input Catalog (Brown et al. 2011). The star was observed for the entire four year duration of the Kepler mission in long cadence mode. The star exhibits stochastic oscillations which confirm that Kepler-91 is an ascending red-giant branch star with a mass of 1.3 M_☉ and radius of 6.3 R_☉ (Huber et al. 2013; Lillo-Box et al. 2014b).

A transiting planet candidate with an orbital period of 6.3 days was detected by the Kepler pipeline (Jenkins et al. 2010; Tenenbaum et al. 2013) and assigned Kepler Object of Interest (KOI) number 2133.01 (Batalha et al. 2013). Lillo-Box et al. (2014b) confirmed the planetary nature of the apparently transiting body. However, the status of this planet has recently been called into question by both Esteves et al. (2013), who use phase variations to deduce that the occulting body is self-luminous, and Sliski & Kipping (2014) who find that the stellar density derived from a transit model differs significantly from the density calculated using asteroseismic techniques. In this paper we present the results of a light curve model combined with radial velocity (RV) observations obtained from the ground and find that the transit-signal is unambiguously caused by a Jupiter-sized planet orbiting the red giant target star.⁷

2.1. Stellar Properties

2.2. Kepler Data

In this work we utilized the full set of Kepler long cadence (29.4 minute) observations from the Kepler spacecraft obtained over 4 yr. These data consist of 17 observational Quarters (Q1–Q17) where all but the first and last Quarter consist of around 90 days of nearly continuous data. Q1 lasted 40 days and Q17 consists of 31 days of data after which Kepler suffered the failure of a reaction wheel.

We used data that have undergone Presearch Data Conditioning (Stumpe et al. 2012; Smith et al. 2012) using the multi-scale maximum a priori (MS-MAP) method (Stumpe et al. 2014). This preprocessing removes signals related to the spacecraft while retaining much of the variability of astrophysical origin. However, the MS-MAP algorithm does not entirely retain astrophysical signals (Thompson et al. 2013). While short timescale events such as transits are largely unaffected, signals on timescales of the orbital period of Kepler-91 (6.3 d) are on average damped by about 5%. The full photometric time series data used in this work is shown in Figure 1.

Zoom In Zoom Out Reset image size

2.3. Radial Velocity Data

We obtained precise RV measurements using the High-Resolution-Spectrograph (Tull 1998) instrumental setup for Kepler follow-up observations and data reduction algorithms as described in Endl et al. (2011). Nine spectra were obtained using an iodine (I₂) cell with an exposure time of 900 s and are shown plotted as a function of time in Figure 1. The data have a resolving power of R = λ/δλ = 30, 000 and were sky background subtracted. We also obtained two template spectra (without the I₂ cell) of Kepler-91, one at R = 30, 000 and a second spectrum with R = 60, 000. The second template was obtained by combining four 900 s exposures and yielded better RV precision. The RV data reported here were obtained using the higher resolution template.

We also measured bisectors and bisector velocity spans (BVS) for the nine spectra used for the RV computation. We measured the bisector and BVS of the cross-correlation-function (CCF) for 11 orders that do not contain significant I₂ lines. We cross-correlated each spectrum with the R = 30, 000 template to search for variability of the BVS that could indicate a false positive and computed the BVS as the velocity difference of two arbitrary points on the CCF bisector at flux values of 0.4 and 0.84, following Hatzes et al. (1998). The BVS results are displayed in Figure 2. Excluding the poor quality measurement from the lowest S/N spectrum, the remaining data have a total rms-scatter of 15 m s⁻¹ with a mean uncertainty of 35 m s⁻¹. The BVS results are consistent with no variability and they also do not correlate with either the orbital phase or the RV measurements. The RV and BVS measurements with uncertainties are listed in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Radial Velocities

Time	Velocitya	Uncertainty	BVS	Uncertainty
(BJD-2454833)	(m s−1)	(m s−1)	(m s−1)	(m s−1)
1208.86670891	114	24	−43	31
1266.71041653	21	18	−15	31
1267.70865078	−25	20	1	32
1268.70698350	−46	17	−31	32
1271.68968502	67	18	−18	28
1275.69264679	−26	25	−3	56
1300.86443801	−9	27	−198b	86
1358.70858740	96	18	−1	32
1382.63282948	0	15	−20	36

Notes. ^aThe values presented here have had an arbitrary offset subtracted to enforce a median of zero. ^bThis outlier comes from the lowest S/N spectrum and is mostly caused by the BVS result from just 2 of the 11 orders used.

Download table as:  ASCIITypeset image

To provide a self-consistent model of both the light curve and the RV observations we chose to model both data sets simultaneously using the orbital model described in Rowe et al. (2014).

Significant planet-induced variability outside of the transit of Kepler-91b has been noted in previous work (Lillo-Box et al. 2014b; Esteves et al. 2013). We chose to model these light curve variations rather than remove them via filtering. We include five physical components in our model of the light curve: a transit, an occultation, ellipsoidal modulation, Doppler boosting, and reflection from the planet. We additionally include the RV data as an additional component in the model and finally we include a model for the correlated noise.

3.1. Parameterization

We used a limb darkened transit model (Mandel & Agol 2002) following a quadratic limb darkening law, and a uniform disk model for the occultation. The ellipsoidal variations, Doppler boosting, and reflection of the planet were modeled in the manner described by Lillo-Box et al. (2014b), but we parameterized the Doppler beaming in terms of K, the RV semi-amplitude, to retain a consistent solution between the light curve data and the spectroscopic RVs. The scaling between the RV semi-amplitude and the Doppler beaming amplitude is proportional to a "beaming factor" B such that

$$\begin{equation} A_b = B \frac{K}{c} (e\cos {\omega }) \end{equation} \tag{ 1 }$$

where A_b is the semi-amplitude of the Doppler beaming signal, c is the speed of light, e is the eccentricity and ω is the argument of periastron. We calculated B in the manner described by Bloemen et al. (2011) and found a value of 5.46 which we kept fixed.

We parameterized the combined model in terms of ρ the mean stellar density, zp a photometric zero point nuisance parameter, linear (γ₁) and quadratic (γ₂) limb darkening coefficients, T₀ the mid-point of transit, P the orbital period of the planet, b the impact parameter, R_p/R_⋆ the planet-to-star radius ratio, eccentricity vectors ecos ω and esin ω where e is eccentricity and ω is the argument of periastron, the amplitude of the ellipsoid variations A_e, the amplitude of the reflection from the planet A_r, the occultation depth F_e, RV semi-amplitude K, and V an RV zero-point. We account for the 29.4 minute integration time of the exposures by subsampling the model 11 times per observation then integrating over these subsamples (Kipping 2010).

Because we are using disparate data sets is it useful to include a parameter for both the light curve and RV data that is an additional noise term which is added in quadrature with the formal uncertainty (σ_lc and σ_rv). These account for missing physics on our model in addition to dealing with underestimation of the reported uncertainties. This creates a flexible model that enables us to scale the two data sets appropriately. We sample in the natural log of σ_lc an σ_rv which is equivalent to using a Jeffery's prior on these parameters.

3.2. Gaussian Process Noise Model

Red giants are known to show significant correlated noise on timescales of hours to weeks due to stellar granulation (Mathur et al. 2012), which is also clearly apparent in Kepler-91. Not accounting for this correlated noise component can bias the observed planet parameters (Carter & Winn 2009).

A Gaussian process is a general framework for modeling correlated noise (Rasmussen & Williams 2006) and it has been applied to transmission spectroscopy (Gibson et al. 2012; Evans et al. 2013) and Kepler light curves (Ambikasaran et al. 2014; Dawson et al., 2014). The basic idea is to model the light curve as being drawn from a large N_data-dimensional Gaussian where the covariance matrix K is modeled as

$$\begin{eqnarray} K_{ij} &=& \sigma _i^2\,\delta _{ij} + k(t_i,\,t_j) \end{eqnarray} \tag{ 2 }$$

where σ_i is the observational uncertainty, δ_ij is the Kronecker delta, and k(t_i, t_j) is the so-called covariance kernel function that quantifies the correlations between data points. There is a lot of freedom in how we choose this covariance function but here we choose a very simple and commonly used model called the squared-exponential or radial basis function kernel

$$\begin{eqnarray} k(t_i,\,t_j) &=& A^2\,\exp \left(-\frac{[t_i-t_j]^2}{2\,\tau ^2}\right) \end{eqnarray} \tag{ 3 }$$

where the covariance amplitude A_GP is measured in flux units and the length scale τ_GP is measured in days. We chose the squared-exponential kernel because it exhibits a number of traits that make it suitable for modeling granulation in red giant light curves: (1) it is smooth with no sharp discontinuities and (2) the covariance between observed points is only a function of the distance (in time) between points. When performing our analysis, we first subtract a realization of the physical light curve model, then calculate the likelihood of the residuals given the covariance matrix defined by the Gaussian Process model. We marginalize over the parameters A_GP and τ_GP with uniform priors in our MCMC analysis.

3.3. Priors

Priors on our model parameters are shown in Table 3. Of particular note is ρ, where we enforce a normal prior constrained from the probability density found from asteroseismology. We enforce a 1/e prior because without this our parameterization in terms of esin ω and ecos ω would bias e high (Eastman et al. 2013). We use a normal prior on limb darkening with the expectation obtained through interpolation of model limb darkening in the Kepler bandpass with T_eff, log g and [Fe/H] fixed at the values shown in Table 3, and with a standard deviation of 0.1. Finally, we set a number of prior constraints on linear combinations of γ₁ and γ₂ that prevent them from taking unphysical values (Burke et al. 2008).

Table 3. Model Parameters

Property	Prior
ρ (g/cc)	$\mathcal {N}(0.0073;0.0001)$
zp	$\mathcal {U}($−0.0002; 0.0002)
γ1	$\mathcal {N}(0.67;0.6)$
γ2	$\mathcal {N}(0.09;0.6)$
T0 (BKJD)	$\mathcal {U}(136.3;136.5)$
P (days)	$\mathcal {U}(6.2;6.3)$
b	$\mathcal {U}(0;1+R_{p}/R_{\star })$
Rp/R⋆	$\mathcal {U}(0;1)$
ecos ω	$\mathcal {U}(-1;1)$
ecos ω	$\mathcal {U}(-1;1)$
e	1/e
Ae (ppm)	$\mathcal {U}(-500;500)$
Ar (ppm)	$\mathcal {U}(-500;500)$
Fe (ppm)	$\mathcal {U}(-500;500)$
V (m s−1)	$\mathcal {U}(-300;300)$
K (m s−1)	$\mathcal {U}(-1000;1000)$
log σlc	$\mathcal {U}(-50;0.0)$
log σrv (m s−1)	$\mathcal {U}(-15;6)$
AGP	$\mathcal {U}(0;0.001$
τGP	$\mathcal {U}(0;0.5)$

Download table as:  ASCIITypeset image

3.4. Markov-Chain Monte Carlo Modeling

We numerically integrated the posterior probability using an efficient affine invariant Markov-Chain Monte Carlo (MCMC) algorithm (Goodman & Weare 2010; Foreman-Mackey et al. 2013). This method utilizes many walkers to reduce autocorrelation time; we opted to use 700 walkers each taking 15,000 steps for a total of 10.5 × 10⁶ samples of the posterior probability. However, we discard the first 5000 samples in each walker as burn-in which leaves 7 × 10⁶ samples used to calculate posterior distributions.

In Figure 3 we show 10 transits seen in the Kepler data, the mean noise model, the light curve model and the 1σ uncertainty on the combined light curve and noise model. The combination of our noise model and transit model does a good job of matching the correlated noise and planet signals seen in the data.

Zoom In Zoom Out Reset image size

We were able to produce a self consistent model that described the data well. Parameters from the modeling are reported in Table 4, we give the median and central 68% bounds of the marginalized posterior distribution for each parameter. We also list the sample that—in our chains—obtained the highest probability. This so-called "best fit" is not the result of an optimization but it does have a self-consistent set of parameters and if you need a single model, this is the one to use. However, if one needs to draw conclusions about the physical properties of the planet then we advise using the median with reported uncertainty.

In Figure 4 we show in the upper left panel the observed data in black (this data still contains correlated noise) with the best-fit model in red and binned observed data in blue. The data have been folded on the orbital period of the planet. The lower left shows the RV data folded on the orbital period of the planet with the same phase as the top right panel. The upper right panel shows the data with a mean Gaussian Process noise model subtracted. Subtracting the mean model gives an idea of what the Gaussian Process is modeling and what noise level remains in the data. We stress that the mean noise model subtracted data should never be used for inference and is for demonstration purposes only. The lower right panel is the same as the upper right panel but zoomed in to better show the ellipsoidal variations and the secondary eclipse.

Zoom In Zoom Out Reset image size

Our parameter estimates are largely consistent with Lillo-Box et al. (2014b). While we obtain a slightly lower estimate for R_p/R_⋆, we are still consistent with Lillo-Box et al. (2014b) at the level of <1.5σ. Our estimate of R_p/R_⋆ is inconsistent with that found by Sliski & Kipping (2014), as shown in Figure 5. The probability that our R_p/R_⋆ estimate agrees with the value found by Sliski & Kipping (2014) is 1%.

Zoom In Zoom Out Reset image size

We derive a planet mass of 0.73 ± 0.13 M_Jup which, combined with a planet radius of $1.308^{+0.061}_{-0.074}$ R_Jup, yields a density of $0.40^{+0.10}_{-0.09}$ g cm³ implying that Kepler-91b is somewhat inflated (Baraffe et al. 2010). The consistency between the planet mass we obtain with RV observations and the mass Lillo-Box et al. (2014b) obtain from the phase curve alone is remarkably good (Lillo-Box et al. determined a mass for Kepler-91b of $0.88^{+0.17}_{-0.33}$ M_Jup). Kepler-91b is one of a growing number of examples demonstrating that hot Jupiter densities can be estimated from Kepler data alone (e.g., Barclay et al. 2012; Quintana et al. 2013; Esteves et al. 2013).

Lillo-Box et al. (2014b) suggest that Kepler-91b may be on an eccentric orbit. In Figure 6 we show our posterior distribution for eccentricity derived from our esin ω and ecos ω samples. While the data prefers an eccentric model, a circular orbit cannot be ruled out.

Zoom In Zoom Out Reset image size

We note that, as mentioned in Section 2, signals on timescales of more than a few days are damped by a few percent owing to the detrending performed within the Kepler pipeline. The only low frequency signal that we significantly detected in these data was from stellar ellipsoidal variations. With an uncertainty of 10% on the amplitude of the ellipsoidal signal, a 5% damping of the signal would not significantly change our results. It should be kept in mind that signals with longer timescales (more than 15–20 days) are significantly damped by the detrending method used within the Kepler pipeline.

The additional uncertainty we report on the RV observations, σ_rv, accounts for astrophysical noise sources that we do not include in our model. One of the sources of astrophysical noise that σ_rv accounts for is the RV jitter induced by the stellar oscillations. The oscillation signal expected is a few meters per second (Kjeldsen & Bedding 1995) and given its short timescale relative to the planetary orbital period acts as an additional noise source on the data. However, photon noise is still the dominant source of uncertainty.

We obtained good fits to both the RV and light curve data with our model as shown in Figure 4. Our Gaussian process model is intended to describe the surface granulation on the red giant. In Figure 7 we show the power spectral density of the original light curve in gray and in purple we show the power spectral density for the light curve with the mean noise model and best fitting transit model subtracted out. We can clearly see that the low frequency noise is well fit by the model. We note that the red giant oscillations are not removed by this model and can be seen as a Lorentzian peak around 100 μHz. In Figure 8 we show the histogram of our flux values from the light curve (with median of zero) and the residuals of flux time series with the mean nose model and transiting planet model subtracted. It is clear here that the residuals with the noise model removed are significantly smaller than in the standard deviation of original time series.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The RV observations phase well with the orbital period defined by the transit and therefore almost certainly are caused by a reflex motion of the star that the planet orbits. The transit model fits the data well with low eccentricity. If the density constraint from asteroseismology was found to be a poor choice of prior we would need a model with high eccentricity, which is not seen. The mass we derive is consistent with a planetary body. We therefore conclude Kepler-91b is a bona fide planet.

The validation presented in this work raises questions as to why Kepler-91 was previously suspected to be a false-positive. Esteves et al. (2013) concluded that Kepler-91b is a false positive based on a phase curve fit where they find a night-side temperature inconsistent with a body that is not self luminous. We do not draw the same conclusions. Using the same equation as Esteves et al. for geometric albedo, we derived a value of $0.52^{+0.22}_{-0.20}$ compared with $2.49^{+0.55}_{-0.60}$ from Esteves et al.

When modeling the correlated noise we find a transit depth that is inconsistent with the transit depth found by both Esteves et al. (2013) and Sliski & Kipping (2014) who both report shallower transits. The method we use to calculate the geometric albedo is

$$\begin{equation} A_g = F_e \left(\frac{a}{R_p}\right)^{2} = F_e \left(\frac{a}{R_{\star }}\frac{R_{\star }}{R_{p}}\right)^{2}. \end{equation} \tag{ 4 }$$

Therefore, underestimating R_p/R_⋆ will overestimate the albedo A_g and lead to a conclusion that the occulting body is self-luminous.

Sliski & Kipping (2014) find that the mean stellar density derived purely from a transit model is not consistent with the asteroseismic density of the star. In Figure 9 we show our estimate of the stellar density and the estimate of Sliski & Kipping (2014), the two estimates differ very significantly. They hypothesize that the observed transit is not caused by a body occulting the target star which is inconsistent with our results.

Zoom In Zoom Out Reset image size

We have hypothesized that the different conclusions regarding the planetary nature of Kepler-91b are owing to correlated noise in the Kepler data. With this in mind, we performed an additional MCMC simulation, but this time did not include our Gaussian Process noise model, and removed the asteroseismic prior on the stellar density. We found that without the Gaussian Process model our posteriors mirrored the results found by Sliski & Kipping (2014) with ρ = $0.043^{+0.001}_{-0.003}$ g cm⁻³ and R_p/R_⋆ = 0.01955 ± 0.0002, in tension to our original parameters when using the noise model (a comparison is shown in Figure 10). This test shows that it is the Gaussian Process noise model that causes the differing results from Sliski & Kipping (2014).

Zoom In Zoom Out Reset image size

We believe it is reasonable to draw the conclusion that not including correlated noise in the model for this star is what led to the erroneous results of the aforementioned teams. However, it is important to note that it is only because we were in possession of the RV data, which raised doubts on the interpretation of Esteves et al. and Sliski & Kipping, that we considered the more sophisticated noise model described here.

We combined Kepler data presented in previous studies of Kepler-91 with RV data obtained from the Hobby–Eberly Telescope. We find that these data can be fit in a self consistent manner and the results lead us to conclude that Kepler-91b is a planet, validating the work of Lillo-Box et al. (2014b). We hypothesize that recent work claiming this as a false positive erred owing to the challenge of modeling a star with a strong correlated noise component and relatively high amplitude out-of-transit variations. We use a Gaussian Process to model the correlated noise present in the light curve of this red giant star and we derive parameters that are in agreement with the planetary nature of Kepler-91b. We suggest that Gaussian Processes may prove to be very useful in future studies of transiting exoplanets orbiting stars exhibiting correlated noise properties.

The validation of Kepler-91 presented in this work confirms that hot Jupiters exist in close-in orbits around red giant stars. Future studies of Kepler-91 and similar planets will allow insights on the internal composition of hot Jupiters as a function of host star evolution (Spiegel & Madhusudhan 2012), and measurements of the spin-orbit inclination of such systems may allow us to test theories of obliquity damping that have been proposed for hot Jupiter populations around less evolved stars (Winn et al. 2010).

We note that after the submission of our paper, independent RV observations confirming the planetary nature of Kepler-91b were published by Lillo-Box et al. (2014a).

This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission Directorate. We would like to express our gratitude to all those who have worked on the Kepler pipeline over the many years of the Kepler mission. Some Kepler data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI). 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. Our MCMC sampling was performed on the Pleiades supercomputer of the NASA Advanced Supercomputing Division at NASA's Ames Research Center. We used data obtained from The Hobby–Eberly Telescope (HET), a joint project of the University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximilians-Universität München, and Georg-August-Universität Göttingen. The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly. We thank the Statistical and Applied Mathematical Sciences Institute (SAMSI) for hosting the Modern Statistical and Computational Methods for Analysis of Kepler Data workshop where the idea for this work came about. The work performed at SAMSI was partially supported by the National Science Foundation under grant DMS-1127914 to the Statistical and Applied Mathematical Sciences Institute. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. We thank Ruth Angus (University of Oxford) for valuable discussions on noise sources in RV observations. We thank Jeffrey C. Smith, Joe Catanzarite (both SETI Inst.), and David Kipping (Harvard/CfA) for suggestions on how to improve the manuscript. E.V. Quintana is supported by a NASA Senior Fellowship at the Ames Research Center, administered by Oak Ridge Associated Universities through a contract with NASA. D. Huber acknowledges support by NASA under grant NNX14AB92G issued through the Kepler Participating Scientist Program.

In this section we have included a number of additional plots to allow the reader to better assess the convergence of our MCMC chains. In Figure 11 we show a trace of each of the 19 parameters included in our model, while Figure 12 shows histograms representing marginalized posteriors of the modeled parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size