Determining the Mass of Kepler-78b With Nonparametric Gaussian Process Estimation

1.2 MB·English

by  Samuel K. Grunblatt

Draft version June 19, 2015
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
DETERMINING THE MASS OF KEPLER-78b WITH NONPARAMETRIC GAUSSIAN PROCESS ESTIMATION
Samuel K. Grunblatt and Andrew W. Howard
InstituteforAstronomy,UniversityofHawaii,2680WoodlawnDrive,Honolulu,HI96822
Raphae¨lle D. Haywood
SUPASchoolofPhysicsandAstronomy,UniversityofStAndrews,NorthHaugh,StAndrews,FifeKY169SS,UnitedKingdom
5
Draft version June 19, 2015
1
0 ABSTRACT
2
Kepler-78b is a transiting planet that is 1.2 times the radius of Earth and orbits a young, active K

n dwarfevery8hours. ThemassofKepler-78bhasbeenindependentlyreportedbytwoteamsbasedon
u radialvelocitymeasurementsusingtheHIRESandHARPS-Nspectrographs. Duetotheactivenature
J of the host star, a stellar activity model is required to distinguish and isolate the planetary signal in

8 radial velocity data. Whereas previous studies tested parametric stellar activity models, we modeled
1 this system using nonparametric Gaussian process (GP) regression. We produced a GP regression of
relevant Kepler photometry. We then use the posterior parameter distribution for our photometric

] fit as a prior for our simultaneous GP + Keplerian orbit models of the radial velocity datasets. We
P
testedthreesimplekernelfunctionsforourGPregressions. BasedonaBayesianlikelihoodanalysis,we
E
selectedaquasi-periodickernelmodelwithGPhyperparameterscoupledbetweenthetwoRVdatasets,
h. giving a Doppler amplitude of 1.86 ± 0.25 m s−1and supporting our belief that the correlated noise
p we are modeling is astrophysical. The corresponding mass of 1.87+0.27 M is consistent with that
−0.26 ⊕
- measured in previous studies, and more robust due to our nonparametric signal estimation. Based
o
on our mass and the radius measurement from transit photometry, Kepler-78b has a bulk density
r
t of 6.0+1.9 g cm−3. We estimate that Kepler-78b is 32±26% iron using a two-component rock-iron
s −1.4
a model. This is consistent with an Earth-like composition, with uncertainty spanning Moon-like to
[ Mercury-like compositions.

2
v 1. INTRODUCTION dius of the planet and the age of the star established by
9 Sanchis-Ojedaetal.(2013),theexpectedDopplerampli-
6 Measuring the radial velocity (RV) of a planet’s host tudewas1-2ms−1,whiletheexpectedspot-inducedRV
3 star is the most common method to measure the mass signal was approximately an order of magnitude larger
0 of a planet. The planetary RV signal, or Doppler ampli- (≈10 m s−1). Howard et al. (2013, H13) and Pepe et al.
0 tude, is directly related to the planet mass. If we know (2013,P13)independentlyobservedandmodeledtheRV
.
1 the mass of the star, we can find the mass of the planet signal of Kepler-78 parametrically to remove any contri-
0 from its Doppler amplitude. If the planet radius is also butionfromcorrelatednoiseactivitywhilemeasuringthe
5 known,wecanthencalculatetheplanet’sdensityandes- planetary Doppler signal. H13 used a sum of sinusoids
1 timate a composition. Precise density and composition atthestellarrotationperiodanditsharmonicstomodel
: information is available only for a handful of transiting the correlated noise activity, accounting for rotational
v
i rocky planets. variability while neglecting any component of the pre-
X dominant signal that was not periodic at precisely the
However, stellar activity can produce spurious RV sig-
r nals larger than some planetary signals. In young, spot- rotation period of the star or its aliases.
a
ted stars, this activity can cause RV variations much P13 noted the fact that the planetary and stellar RV
larger than most observed planetary RV signals (Hillen- signal timescales differ by over an order of magnitude:
brandetal.2015),andeveninrelativelyquietstarssuch theorbitalperiodoftheplanetisonly8.5hours,whereas
as the Sun, these variations are at least one order of the stellar activity-induced signals are modulated by
magnitude larger than the RV signal of Earth (Meunier the stellar rotation period of 12.8 days. This allowed
et al. 2010). Before we can have any hope of confirming P13 to test a floating chunk offset model in which the
the discovery of an Earthlike planet around any star, we free parameters were the planetary Doppler amplitude
must account for the contribution of stellar activity. and an RV zero point to represent the noise signal, as-
sumed constant for each night of observations. How-
Discovered in 2013, the radius of Kepler-78b was mea-
ever, P13 found that the evidence ratio for their float-
sured to be 1.16 ± 0.16 R (Sanchis-Ojeda et al. 2013).
⊕ ing chunk offset model was significantly lower than the
Subsequently, the mass of the planet was measured and
sum of sinusoids model, given the large number of free
reported simultaneously by two teams. Based on the ra-
parameters, and used the parametric, sinusoidal model
Email: [email protected] for their reported mass measurement. Howard et al.
2 Grunblatt et al.
(2013)used Keck/HIRES and measured the mass of Ke-
TABLE 1
pler 78b to be 1.69 ± 0.41M , while Pepe et al. (2013)
⊕ Stellar Properties
used TNG/HARPS-N and measured a mass of 1.86+0.38
−0.25
M⊕. Despite the difference in their methods and obser- Property Value
vational techniques, the H13 and P13 planet mass mea-
Name Kepler-78/KIC8435766/Tycho3147-188-1
surements are consistent, attesting to the robustness of
Age 625±150millionyears
the results. Vsin(i) 2.6±0.5kms−1
Further exploration into nonparametric estimation is Mass,Mstar 0.83±0.05M(cid:12)
Inclination,i 75.2+2.6 degrees
justified by the limited scope of the purely parametric −2.1
RotationperiodProt ≈12.5days
modelstestedpreviously. Wedrawattentiontothecom-
parable case of HST transmission spectroscopy of hot
Note. —ThevaluesinTable1arehavebeentakenfromHowardet
Jupiter HD 189733, studied by Swain et al. (2008), Gib- al.(2013).
son et al. (2011) and Gibson et al. (2012). Gibson et al.
(2012) note that both previous studies of the transmis- 2.1. Kepler photometry
sion spectroscopy used linear basis functions to account
forsystematicerrors,andarguethatthisisnotsufficient The Kepler telescope obtained photometry of approx-
to account for instrumental systematics, and therefore imately 150,000 objects in the Kepler field for its four
provides an unrealistic treatment of the uncertainties. year lifetime from 2009 to 2013. Photometry of Kepler-
Gibson et al. (2012) reanalyzed the spectroscopic data 78 was gathered at a 30-minute cadence. In this study,
with a Gaussian process to marginalize ignorance of the we train a GP on the photometric light curve of Kepler-
functionalformofthesystematics,givinglargerbutmore 78 to determine the evolution and rotation timescales of
robusterrors. Similarly,wehavereanalyzedtheRVdata the stellar activity, hyperparameters of the GP regres-
of Kepler-78b with several Gaussian process estimators sion kernel. However, as shown in Equation 2, a matrix
to find a model of the RVs with the highest evidence, inversion is required to calculate the log posterior likeli-
resulting in more robust uncertainties of the mass, den- hoodofaGPkernelwithagivensetofhyperparameters,
sity, and composition of Kepler-78b. In this case, we are a computationally intensive process with compute time
usingGaussianprocesstodescribeasignalofastrophys- proportional to N3. Therefore, it was necessary to rebin
ical origin, rather than instrumental, which we confirm theKepler photometrytoonepointevery5hours(aver-
through tests described in our Results section. agingeverytenpointstogether)inordertofindthebest-
fit hyperparameters of a full quarter of photometry with
In this study, we combine the RV measurements from
our computing resources. While this does marginalize
H13 and P13. We describe each radial velocity dataset
over the planetary signal, we are only using the photom-
with a Gaussian process (GP) regression combined with
etry to estimate the underlying components of the stel-
a Keplerian orbit signal at the orbital period and phase
laractivitywhichvaryontimescalesmuchlongerthan5
oftheplanet(bothknownpreciselyfromthetransitpho-
hours, and the photometric effect of the planet is small
tometry). GP regression is a nonparametric method for
compared to that of the stellar activity. This photome-
modelingcorrelatedRVnoise,andcanprovideamorero-
try, taken from the 16th quarter of Kepler observations,
bustnoisemodelbecauseofitsflexibilityoverparametric
is not concurrent with the RV measurements; however,
models. The observations analyzed are described in §2.
since it was taken only ≈55 days before the RV mea-
Our analysis of the combined RV dataset is described in
surements, we can assume that the variation seen in the
depth in §3, starting with a review of Gaussian process
photometry is at least related to the RV variation, and
regression and the benefit of a nonparametric model in
would exhibit similar structure in the temporal dimen-
§ 3.1-3.3, and a focus on how different GP models were
sion. Figure 1 illustrates the binned photometric mea-
tested and chosen for the data in sections § 3.4 and 3.5.
surements as well as the GP regression with the best-fit
In §4, we report our results for the planet Doppler am-
kernel hyperparameters to the photometry.
plitude and mass, and discuss the selection of the cho-
sen analysis model. Using the known properties of the
hoststar(summarizedinTable1)andtheplanetradius, 2.2. HIRES and HARPS-N RV data
we recalculate the planetary density. In §5, we compare
The RV observations were obtained shortly after the
our mass measurement and density calculation to pre-
malfunctioning of Kepler’s reaction wheels in summer
vious results and discuss possible composition scenarios
2013, which prevented the high level of photometric pre-
for Kepler-78b. We explore the applications of this tech-
cision achieved in the first era of the spacecraft’s life-
nique to other RV datasets in §6.
time and required that the telescope point at new tar-
gets. H13 observed Kepler-78 with HIRES on the Keck
I telescope on Mauna Kea. They obtained eight nights
2. OBSERVATIONS of data during each of which approximately 9 or more
measurements were obtained: three measurements were
InordertodeterminetheDoppleramplitudeofKepler- takenatsuccessivethirtyminuteintervalsatthreesepa-
78b, it was necessary to first confirm the general struc- rateoccasionspernight. Inaddition,fiveindividualmea-
ture of the stellar activity by modeling the Kepler pho- surementsweretakenonseparatenightsaftertherestof
tometry. Thesemodelparameterswerethenusedtopro- the data were collected to provide a longer baseline over
vide reasonable initial hyperparameters for the RV anal- which to evaluate the stellar activity. The dataset con-
ysis. The observations taken are described below. sists of 84 RVs over 45 nights. The P13 team observed
Kepler-78b Mass Measurement with Gaussian Process 3
Kepler-78 with HARPS-N on the Telescopio Nazionale spot models such as those used by H13 and P13 gave in-
Galilei in the Canary Islands. They were able to obtain correct masses and uncertainties. Thus, it is important
109 measurements over 97 nights. Most of their obser- to test many time series techniques, and further explore
vations were made in a consecutive six-day period. On the novel application of GPs in Doppler analysis. Weex-
subsequent nights only one or two measurements were pand upon the GP training technique used in Haywood
made per night, giving a total time baseline of approx- et al. (2014) here.
imately three months. Figure 2 illustrates the HIRES
and HARPS-N measurements with measurement errors
3.2. Gaussian Process Covariance Kernel Choice
shown, with GP regressions overplotted.
Finding the best GP regression requires choosing a
3. METHODS kernel and initial hyperparameters, evaluating the likeli-
hood of those hyperparameter values, and then iterating
We model the photometric variations observed by Ke- through parameter space until the most likely values are
pler with a GP regression with a quasi-periodic kernel. found. The squared exponential kernel, for example, de-
We then model the radial velocity measurements with a fines a covariance matrix through an operator,
quasi-periodic kernel GP regression, using the kernel pe-
riod hyperparameter, θphot, to train the RV GP period (cid:20) (cid:16)t −t (cid:17)2(cid:21)
hnyapnterspiagrnaaml eotbesrerθvReVd. inThbiosthassduamtaessettshactomthees fprroemdomthie- Σij =k(ti,tj)=h2exp − i λ j , (1)
star, and thus will modulate at the stellar rotation pe-
wherehisthecovarianceamplitude,andλthecovariance
riod. Whenwemodelthephotometricdata,werecovera
length scale. The amplitude observed is described by h,
periodhyperparameterthatagreeswithotherautocorre-
while λ is a characteristic timescale over which the data
lationanalysesoftheKepler photometry(Sanchis-Ojeda
is going to be correlated.
etal.2013),supportingthisassumption. WecreateaGP
regression of each RV dataset with a quasi-periodic ker- We discuss other GP kernels and the inferred physical
nel, where the period hyperparameter has been trained meaning of their hyperparameters in Table 2.
ontheperiodhyperparameterofthephotometricregres-
The logarithm of the posterior likelihood of the GP
sion. We then fit all the RV GP hyperparameters, h,
regression is calculated as
θ, w, and λ in the quasi-periodic case, σ, a white noise
parameter, as well as the Doppler amplitude K of a Ke- 1 1 n
plerian signal at the known orbital period and phase of log[L(r)] = − rTΣ−1r − log|Σ| − log(2π), (2)
2 2 2
the planet, to the RV data simultaneously. We find that
all of our GP + Keplerian models are able to recover where r is the vector of residuals after removal of the
theplanetaryDoppleramplitudeK atavalueconsistent (optional) mean function, Σ is the covariance matrix,
with previous work. and n the number of data points. A prior term, Lprior,
can be added to the likelihood to account for any priors
placed on the hyperparameters. For example, we apply
3.1. Gaussian Process Regression: Concept
the Gaussian prior
Gaussian process is a nonparametric method to de- −1(cid:104)(cid:16)θtrue−θ(cid:17)2(cid:105)
scribe a dataset by evaluating correlations between n L =e 2 σθ , (3)
prior
data points through a covariance kernel. This kernel
describes the relationship of each point in the dataset to restrict the hyperparameter θ for the RV data, as the
to each other point, and can be expressed as an n×n period of the correlated noise signal is equivalent to the
matrix (subsequently referred to as the covariance ma- stellar rotation period Prot in both the photometric and
trix). Thekernelisafunctionofhyperparameters. More RV datasets. We restrict the period hyperparameter θ
complicatedkernelscanhavemorehyperparametersthat in the RV regressions within a Gaussian of width σθ de-
characterize different qualities of the correlations in the termined by the posterior distribution of the period hy-
data, such as various periods, characteristic amplitudes perparameter found for our photometric GP regression.
and length scales, etc. Prior knowledge of this hyperparameter helps to ensure
convergence of the other hyperparameters of the RV re-
Gaussian process regression is widely used in the field gression, as the RV data is not as well sampled as the
of machine learning (Neal 1997; Herbrich et al. 2003; photometry.
Quin˜onero-Candela & Rasmussen 2005; Wang et al.
2008). Gibson et al. (2012) introduced the technique to This likelihood calculation can be used to identify the
the field of exoplanets through analysis of transmission best fit GP hyperparameters. For more details on how
spectroscopytomodelcorrelatednoiseintheinstrumen- thecovariancekernelandparameterboundaryconditions
tal systematics of HST/NICMOS, as described in §1. werechosen, referto§4.1. Foramorecompletedescrip-
Concurrentlywiththiswork,Haywoodetal.(2014)have tion of Gaussian process regression and posterior likeli-
demonstrated the technique of GP modeling of RV and hood evaluation, see Rasmussen & Williams (2006).
photometric signals for the CoRoT-7 planetary system,
first modeling the photometry with a GP and then us- 3.3. Gaussian Process and Stellar Activity
ing the photometric GP hyperparameters to train the
initial RV GP hyperparameters. Haywood et al. (2014) The presence of magnetic features on the stellar sur-
demonstrated that in the case of CoRoT-7b, parametric face induces variation in a star’s RV that can mimic a
4 Grunblatt et al.
Fig. 1.—FluxofKepler-78versustime,asrecordedbytheKepler spacecraftduringits16thquarterofobservations. Thedataareplotted
asblackpoints,withphotometricerrorsshown. ThebluelinecorrespondstotheGaussianprocessregressionofthephotometrywiththe
best-fitkernelhyperparameters,andtheshadedregionscorrespondto1-σ and2-σ uncertainties,asdefinedbytheposteriordistributions
of the kernel hyperparameters. The Kepler data has been shown as fitted with a single quasi-periodic kernel function and a white-noise
stellarjitterparameter.
TABLE 2
Gaussian Process Kernel Options
Name Mathematicalexpression Hyperparametersa Comments
Squaredexponential h2exp(cid:20)−(cid:16)ti−tj(cid:17)2(cid:21) h,λ hamplitudeofcovariancefunction,
λ
λacharacteristictimescale
(cid:20) (cid:21)
Periodic h2exp − sin2[π(2twi−2tj)/θ] h,θ,w θ equivalenttoProt,
w representscoherencescale,similartoλexpressed
asafractionofθ dependentonrecurrentfeatures
Quasi-Periodic h2exp(cid:20)− sin2[π(ti−tj)/θ] −(cid:16)ti−tj(cid:17)2(cid:21) h,θ,w,λ w coherencescaletiedtoperiodicvariation
2w2 λ
whilecharacteristictimescaleλtiedtoaperiodicvariation.
Note. —ThenameofkernelfunctionsandhyperparametersinTable2aretakenfromRasmussen&Williams(2006).
a EachkernelΣij canbemodifiedtoincludeanadditionalhyperparameter,awhitenoisetermσ2 byaddingoneinquadrature: Σij=Σij+σ2Ii.
planet’s orbital RV signal. Starspots are regions of high on the stellar surface while keeping a statistical “mem-
magnetic fields on the surface of a star that appear rel- ory” of the stellar activity patterns.
atively cooler and darker than the surrounding photo-
sphere. As the starspots move across the disk of a star
(due to stellar rotation), the flux balance between the 3.4. Photometric Model
redshifted and blueshifted halves of the star is broken,
producing a shift in the centroid of the spectral lines. WefitthephotometricdatausingGPregressionswith
This shift corresponds to an apparent RV signal and three different kernels. We test a squared exponential,
change as the spots move across the stellar disk (Du- periodic and quasi-periodic kernel (Table 2). We test
musque et al. 2011). Starspots, as well as networks of kernels with and without an additional white noise term
strongly magnetizedflux tubes knownas faculae, inhibit addedinquadraturetothelikelihood(Eq. 2). Thebest-
the convective processes taking place on the stellar sur- fit kernel hyperparameters are found via Markov Chain
face, thus reducing the net blueshift produced by the up Monte Carlo analysis powered by the emcee Python
flow of hot, bright granules. This has been shown to be package(Foreman-Mackeyetal.2013),describedinmore
thedominantRVeffectontheSun(Meunieretal.2010) detail in the following section. We select the quasi-
andonotherSun-likestars(Haywoodetal.2014). These periodic model with a white noise term based on our as-
activity-induced RV signals are modulated by the stellar sumptionthatthepredominantsignalinbothdatasetsis
rotation and the surface features evolve with time, re- tied to the quasi-periodic variation of the stellar surface,
sulting in a quasi-periodic, RV signal. The lifetimes of andconfirmthatthequasi-periodicmodelisappropriate
surface features are poorly constrained, making this ef- through a visual check and a reduced χ2 analysis. We
fect difficult to model parametrically. The nonparamet- have plotted the data and GP regression with the best-
ric nature of the GP regression provides an ideal frame- fit kernel hyperparameters in Figure 1. We ensure that
work for studying stellar surface features in RV. We can these hyperparameters are robust by performing a GP
choose a GP kernel that reflects the frequency structure regression to the photometric data after less stringent
ofthestellaractivity. TheresultantGPmodelisflexible binning, and recovering the same best-fit kernel hyper-
enoughtoaccountfortheevolutionofmagneticfeatures parameter values within errors.
Kepler-78b Mass Measurement with Gaussian Process 5
3.5. RV Model ouranalysis(e.g.Knutsonetal.2008). Tocheckthatthe
resultant fit is a good description of the data, we com-
WefitourphotometricandRVdatasetswiththesame pute a χ2 for the best-fit model of each kernel variant.
red
kernel operator. We took the hyperparameter values Finding χ2 ≈ 1 confirms that the quasi-periodic GP
red
from the photometric GP regression and used them as model and errors describe the data well.
initial fit values for the RV GP kernel hyperparameters,
with the exception of the hyperparameters h and σ, as WeallowasmanyoftheMCMCparametersaspossible
the variation in the photometry and RV related to these beminimallyconstrained. Wefindthatwheneverthepe-
quantities is not in equivalent units. We also place a riodicorquasi-periodickernelisused,andthephotomet-
Gaussianprioronθbasedonitsposteriordistributionde- ric θ value is not applied as a prior, the RV noise model
termined from the photometric GP regression. We then is able to recover periodicity consistent with the stellar
apply a uniform offset to each set of RV measurements rotation period. This supports our assumption that the
to remove RV zero point differences. stellarsurfacefeaturesarethepredominantsignalseenin
theRVmeasurements. ApplyingaGaussianpriortothis
The likelihood of the GP regression, given certain ker- value allows us to better constrain all other MCMC pa-
nelhyperparametersandaDoppleramplitudevalue, de- rameters. Physicalboundaryconditionswerealsoplaced
termines the quality of the model for those given pa- on some hyperparameters to prevent the MCMC chain
rameters. This likelihood can be calculated as given in from moving into unphysical parameter space (see Table
Equation 2, where the residuals r can be calculated as 4 for details).
(cid:16)2π(t−t )(cid:17)
r=v−Ksin c , (4)
P 4. RESULTS
orb
where t is the vector of times of all measurements, v the We find a Doppler amplitude of Kepler-78b of 1.86 ±
vectorofRVmeasurements,tc atimeoftransit,Porb the 0.25 m s−1, corresponding to a mass of 1.87+0.27 M
orbital period, and K the planetary Doppler amplitude. −0.26 ⊕
using the quasi-periodic kernel with common temporal
Sanchis-Ojedaetal.(2013)measuredt to10−5 ofaday
c hyperparameters, excepting amplitudes and white noise
andandPorb to10−7 ofadayprecision, sowecansafely terms, for the two RV datasets. We find that the quasi-
treat them as constants for this analysis. periodic kernel provides the best fit to the Kepler pho-
tometrybasedonvisualsignalinspection, andthevalid-
The best-fit GP kernel hyperparameters and the Ke-
ityofthefitisconfirmedthroughχ2 calculation. Based
plerian Doppler amplitude are found via MCMC explo- red
onaBayesiananalysis,wefindthestrongestevidencefor
ration of parameter space. We simultaneously fit the
a quasi-periodic kernel regression of the radial velocity
Keplerian planetary signal and a GP regression to each
data over simpler kernels as well, and that stronger evi-
RV dataset with the Doppler amplitude and kernel hy-
dence exists for one set of temporal hyperparameters to
perparameters as free parameters in the MCMC chain.
describebothradialvelocitydatasetsasopposedtotwo,
The emcee package contains an Affine-invariant Monte
suggestingthatthedominantsignalinbothRVdatasets
Carlo Markov Chain Ensemble sampler, which evaluates
is indeed astrophysical. This supports the findings of
the likelihood of the GP kernel hyperparameters to the
Haywood et al. (2014). We provide the adopted model
measurement residuals after a Keplerian planetary sig-
parameters for the best-fit quasi-periodic RV GP regres-
nal has been subtracted (Foreman-Mackey et al. 2013).
sion + Keplerian orbit models chosen in Table 4.
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 4.1. GP kernel selection
from the posterior distributions generated through this
MCMC exploration of parameter space, with 1-σ error We find a planetary mass within 1-σ of our adopted
correspondingto68%confidenceintervalsintheMCMC mass result for all models tested. We make the assump-
posteriordistributions. Weplotthedataandbest-fitGP tion that the predominant signal observed in the photo-
regression + Keplerian models in Figure 2. metric and both RV datasets comes from the star, and
thus will modulate at the stellar rotation period. We
AftertheMCMCposteriordistributionshavebeencre-
find that when we model the photometric data with a
ated, we calculate a Gelman-Rubin statistic for each pa-
periodic or quasi-periodic GP kernel, we recover a pe-
rameter to ensure that the parameter chains have con-
riod hyperparameter consistent with previous estimates
verged. Convergence is deemed adequate when the G-R
ofthestellarrotationperiod,supportingourassumption
statistic < 1.01 (Gelman & Rubin 1992). In addition,
about the predominant photometric signal. Similarly,
we track the average and median log likelihood at each
when we model the RV datasets with a quasi-periodic
stepinthechain. Atfirst,thechainsmovetohigherand
GP with no priors, we recover a period hyperparameter
higherlikelihoodparameterspace. Oncethebest-fitsolu-
consistent with the stellar rotation period or its alias.
tionisfound,thechainsbegintotakemoreandmoreun-
Thus,weconcludethatthequasi-periodicGPregression
likely steps, exploring parameter space more completely,
is both well motivated and effective for describing the
and pushing the average likelihood below the median
predominant signal in our datasets.
value. The step at which the average log likelihood be-
gins dipping below the median log likelihood indicates We calculated BIC and AIC values for all of our RV
the transition from burn-in period to the exploration of models. We test the three kernels of interest both with
parameter space around the maximum likelihood solu- shared hyperparameters as well as with fully indepen-
tion,andweremovethestepstakenduringburn-infrom dent hyperparameters in order to ensure that the varia-
6 Grunblatt et al.
Fig. 2.— RV of Kepler-78 versus time, measured by Keck-HIRES and HARPS-N. The HIRES and HARPS-N data are plotted as
blue and red points respectively, with errors in RV shown. A Keplerian orbit signal with the calculated best-fit Doppler amplitude has
been subtracted from the data. The colored lines correspond to the Gaussian process regressions with best-fit kernel parameters of the
correspondingly colored RV measurements, where the shaded regions correspond to 1-σ and 2-σ uncertainties. Both datasets have been
fitted with a single quasi-periodic kernel operator with common period, roughness, and lengthscale (θ, w, and λ) hyperparameters, but
separatecovarianceamplitudeandwhitenoise(handσ)parameters.
ther of the other models, with a ∆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−∆AIC/2, indicating
that the uncoupled quasi-periodic model is 1% as likely
as the coupled model, and the coupled squared exponen-
tial 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 ker-
nels. Strictlyspeaking,theseinformationcriteriarequire
theassumptionthatthenoiseinthedataisindependent
and identically distributed, which we know is not the
Fig. 3.—TheresidualsoftheHIRESandHARPS-Ndataafter case. However, since we use the same data for all BIC
the quasi-periodic GP regression to the data with the Keplerian
and AIC comparisons, they provide a valuable compar-
signal subtracted is removed, phase-folded at the known orbital
periodoftheplanet. HIRESdataisshowninblue,andHARPS-N ison between the different GP kernels. We report all
dataisshowninred. Theplanetarysignalmodelisshownbythe relevant AIC and BIC values in Table 5.
dottedblackline,andbinneddataareshownbyblackpoints.
Ifwearemodelingthestellaractivityinbothdatasets,
tion we observe is consistent between datasets and thus we would expect the hyperparameters for our simulta-
astrophysical, rather than local, in origin. Furthermore, neous models to be the same, with the exception of the
we test our squared exponential kernel models with and non-temporalamplitudehyperparameterhandthewhite
without a white noise term for each dataset. We report noise term σ, which would be different but related be-
these BIC and AIC values in Table 5. tween the models due to the different appearance of the
disk of the star in different wavelength regimes. We ex-
We find that the quasi-periodic model we choose has
plored the relation of the GP kernel amplitude hyperpa-
AIC and BIC values with ∆BIC ≥−0.5 and ∆AIC ≥8
rameters,whoserelationshipweevaluatedbycalculating
between it and all other models tested. A ∆BIC value
a, the ratio of the HARPS-N amplitude to the HIRES
of less than 2 suggests two models are indistinguishably
amplitude. The convergence of a to a best-fit value sug-
likely, whereas a ∆BIC greater than 10 suggests strong
gests that there is a steady relation in the ranges of the
evidence for the model with a lower BIC value. Interest-
RV signals observed by HIRES and HARPS-N. Addi-
ingly, the BIC for the quasiperiodic and squared expo-
tionally, the value of a is broadly consistent with the
nential kernels are almost identical, while there is signif-
results of Desort et al. (2007), who illustrate that the
icantlylessevidencefortheperiodicmodel(BICofperi-
expected radial velocity contribution from starspots on
odic model 18 points higher than SE and quasi-periodic
a K dwarf is larger at shorter wavelengths, resulting in
models). Furthermore, the lowest BIC of all was found
a larger signal in HIRES than in HARPS-N, as HARPS-
forthesquaredexponentialkernelwithtwoindependent
N has significant wavelength coverage in the 600-700nm
GPregressions,witha∆BIC=−0.5ascomparedtothe
regime that HIRES does not. We show the distribution
chosen model. However, the AIC of the quasi-periodic,
of a along with the other parameters of our MCMC cal-
coupled GP model is significantly lower than that of ei-
Kepler-78b Mass Measurement with Gaussian Process 7
TABLE 3 TABLE 4
Planetary Properties Adopted Quasi-periodic RV Model Parameters
Property Value Name Prior Value
Name Kepler-78b HIRESAmplitudeh h>0 11.6+3.7 ms−1
OrbRitaaldipuesr,ioRdplPorb 0.35500714.240±±00.0.0090R00⊕006days HARPS-NAmplitudeh2 h2>0 5.6+−−212..0.35ms−1
Doppleramplitude,K 1.86±0.25ms−1 Amplituderatioa - 2.0+0.8
−0.5
DMenasssit,yM,ρpl 61.0.8+71+−.900..g2276cmM−⊕3 Periodθ e−(123(.02.61−2)θ2)2 13.12+−00..1142 days
Ironfractiopnl −321.±4 26% Coherencescalew - 0.28+−00..0054
Characteristictimescaleλ λ>0 26.1+19.8 days
−11.0
Note. —Thename,radius,andorbitalperiodvaluesinTable3are HIRESwhitenoiseσjitter σ>0 2.1+−00..33 ms−1
takenfromHowardetal.(2013). Allothervalueshavebeencalculated HARPS-Nwhitenoiseσjitter σ>0 1.1+−00..45 ms−1
forthiswork. Doppleramplitude,K - 1.86±0.25ms−1
culation in Figure 4. Any variation due explicitly to the
correlated stellar activity should thus cause our models to the range of the RV shift observed, and thus difficult
tohaveashape(periodandphase)thatisconsistentbe- tointerpretphysically. Furthertestsoftheseparameters
tween the two datasets. We allow the reader to confirm arenecessaryinordertofullycharacterizetheirrelation-
this visually in Figure 2. ships, such as placing new priors on the coherence scale
and characteristic timescale parameters to ensure that
Distributions for each parameter relative to all other
they are not biased by noise.
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 5. INTERPRETATION AND DISCUSSION
tested, indicating that the Doppler amplitude observed
is a real signal and not a systematic error introduced Withourmeasuredmassof1.87+0.27 M andaradius
−0.26 ⊕
during our parameter fitting. of 1.20 ± 0.09 R (H13), Kepler-78b has a bulk density
⊕
of6.0+1.9 gcm−3,suggestingarockycompositionsimilar
We place a Jeffreys prior on the period hyperparame- −1.4
terθ toweightshorterperiodsmoreheavilywhenfitting toEarth(ρ⊕ =5.52gcm−3). Usingthetwo-component,
the quasi-periodic and periodic GP regression kernel hy- rock-ironmodelsfromFortneyetal.(2007), weestimate
perparameters to the photometric data (Haywood et al. anironfractionof32±26%,consistentwithanEarth-like
2014). When we do this, we are able to recover the ro- composition (iron mass fraction of 0.319, McDonough &
tation period of the star as the period hyperparameter. Sun 1995). These simplified models consider Kepler-78b
We then place a Gaussian prior on this hyperparame- as an iron core surrounded by a rocky mantel, and ac-
ter in the RV GP + Keplerian model constraining it to countforcompressionthatisimportantforhigherplanet
the value and errors measured in the photometric data. masses. Weignoretheeffectofanatmosphereonthera-
We also note that although the coherence scale parame- dius of Kepler-78b due to its equilibrium temperature of
ter w measured in the photometric data is equal to that 1500-3000 K (Sanchis-Ojeda et al. 2013).
measured in the RV data within errors, we do not place
Wenowspeculateabouttheimplicationsofthesemea-
a similar Gaussian prior on this parameter, because we
surementstothecompositionofKepler-78b. InFigure5,
observethatitiscorrelatedwithboththeamplitudeand
we illustrate that the best-fit composition of Kepler-78b
characteristictimescalehyperparameters,whichwerenot
isindistinguishablefromthecompositionofEarth. How-
necessarily related between the photometric and RV GP
ever, Kepler-78b’s range in iron mass fraction stretches
regressionhyperparameterfits. Thiscorrelationbetween
from almost purely rock (a Moon-like composition) to
h, h , w and λ visible in Figure 4 is hard to interpret
2 60% iron (a Mercury-like composition). This is consis-
physically. Wespeculatethatthepositivecorrelationbe-
tent with a range of rocky solar system planets but dis-
tween w and the h parameters may arise from a connec-
tinguishable from other solid celestial bodies (such as
tion to starspot size: as starspots grow larger, h grows,
M-type asteroids, which are almost 100% iron by mass,
anywhitenoisepresentwillbecomerelativelylessimpor-
or comets, made up of rock and ice with no iron at all).
tant, and thus the function will become smoother over-
Since its best-fit composition and mass are more simi-
all,increasingthew parameter. Similarly,thecharacter-
lar to Earth than any other solar system body, Kepler-
istic timescale λ may also increase at large amplitudes
78bmighthavehadaformationprocessthatwassimilar
for the same reason–as the uncorrelated noise at small
to Earth’s. The proximity of Kepler-78 to its host star
timescalesbecomeslessimportant,thelengthscalemight
(0.009AU)makesitunlikelythatitformedinitscurrent
be weighted more heavily toward larger values. Larger
location, and migration to its current orbit is likely.
spots also persist longer on the stellar surface, as spots
decay on timescales proportional to their size (Bumba Theimprovementinmassdeterminationofthismodel
1963). However, the fact that w was consistent for the relative to H13 is shown in Figure 5. Despite the fact
photometric and RV datasets whereas λ was smaller for that the mass measured by this method is a 6.5-σ mea-
the photometric datasets despite their large differences surementcomparedtoH13’s4-σandP13’s6-σdetection,
in amplitudes may suggest otherwise. In addition, the the new measurement of Kepler-78b’s density is as pre-
amplitude parameters h and h are not directly related cise as those obtained by either of the two competing
2
8 Grunblatt et al.
Fig. 4.—ParameterdistributionsfromourchosenGP+KeplerianmodeloftheRVdataplottedagainsteachother. Themedianvalue
ofeachparameter(bluelines)and1-σ ranges(dashedlines)areshown. Theparametera=h/h2.
teams in 2013, who estimated its density to be 5.3+2.0 stellar mass.
−1.6
(H13) and 5.57+3.0 g cm−3 (P13). Relative to P13, the
−1.3 Hatzes (2014) tested several traditional models on the
errors on the density in this study are 15% smaller, al- same RV data analyzed here to check the robustness of
thoughonlymarginallysmallerthanthatofH13. Thisis the planetary detection and verify whether the planet
becauseeventhoughthenewmassmeasurementissome- could be found without previous transit knowledge. A
whatmoreprecise,theerroronthedensityisdriventhree range of Doppler amplitudes of the planet from 1.31 to
times more strongly by the error on the radius than the 1.96ms−1 arereported,whichisbroadlyconsistentwith
mass. Highercadencephotometryoftheplanetarytran- the result of this work as well as that of H13 and P13.
sitwouldallowforingressandegressofthesystemtobe Hatzes (2014) reports the planet can be identified with-
observed, breaking the degeneracy between the impact out prior knowledge of the system using a modified ver-
parameter and transit depth. Estimates of the planet’s sion of the parametric method originally tested by P13,
density could also be improved with more accurate mea- thefloatingchunkoffsetmodel,toanalyzebothdatasets,
surementsofthestellarmass,asconvertingRVmeasure- andreportsamassof1.31±0.24M ,inconsistentatthe
⊕
mentsintoaplanetarymassisdirectlydependentonthe 1-σ levelwiththemasscalculatedinboththisstudyand
Kepler-78b Mass Measurement with Gaussian Process 9
TABLE 5
Model MCMC Fit Diagnostics and Results
Kernelfunction K(ms−1) BICvalue(2GPs)b AICvalue(2GPs)b Comment
Quasi-periodic(adopted)a 1.86±0.25ms−1 1065.3(1079.6) 1039.2(1047.0) withwhitenoise,commonw,θ,λ
butdifferenth,σ parameters
Periodica 1.82±0.29ms−1 1083.8(1090.6) 1067.5(1061.2) withwhitenoise,commonw,θ,
butdifferenth,σ parameters
Squaredexponential 1.92±0.27ms−1 1066.4(1064.8) 1053.4(1058.4) withwhitenoise,commonλbut
butdifferenthparameters
a Fortheperiodicandquasi-periodickernelmodels,theperiodhyperparameterθ isconstrainedbyaGaussianprior. Thispriorisshapedbythe
posteriordistributionofthecorrespondingphotometricGPhyperparameter. Whenthispriorwasnotapplied,thebest-fitperiodfoundintheRV
datawasconsistentwiththestellarrotationperiod,butDoppleramplitudeerrorswerelarger.
b WemeasuredtheAICandBICvaluesforallmodelswithcommonRVGPhyperparametersandwith2independenthyperparametersets,which
wecall2GPshere. WealsomeasuredBICandAICofthe2GPsquaredexponentialkernelwithoutawhitenoisetermbutdonotincludethemin
thistable.
Kepler-78 system, and can obtain a mass measurement
for Kepler-78b that is consistent with previous work.
6. SUMMARY
We performed a combined analysis of the photometric
and RV data of Kepler-78 in order to better extract the
RVsignaloftheplanetKepler-78b. Wefitthedatausing
Gaussianprocessregression,andtestthreesimplekernel
configurations. Aftertestingmultiplemodels,wefindthe
strongest evidence for squared exponential kernel mod-
els, and a quasi-periodic kernel model in which the hy-
perparametersfortheGPregressionforeachRVdataset
are coupled (aside from the amplitude hyperparameter
and the white noise term). We prefer the quasi-periodic
model due to evidence indicated by AIC calculation and
because it supports our understanding of the physical
origins of the observed signals. We measure the Doppler
amplitude and therefore the mass of Kepler-78b to 6.5-
σ significance, comparable to or better than all previous
massmeasurementsofthisplanet. Weconstraintheiron
massfractionoftheplanetto32±26%, illustratingthat
Fig. 5.— A graph of mass versus radius for rocky planets. The
Kepler-78b is most likely Earth-like in composition.
constraints on Kepler-78b are shown by the red point. The es-
timation of Kepler-78b by H13 is shown as the light gray error
Theanalysisdoneinthiswork(andpreviousstudiesof
bars. Compositioncurvesrangingfromapurewater(blue)topure
ironcomposition(black)havebeenplotted. Earth-like(67%rock, this system) is possible because the orbital period of the
33% iron) and Mercury-like (67% iron, 33% rock) compositions planet P is an order of magnitude smaller than and
orb
are denoted by green and brown curves, respectively. Solar sys-
not a harmonic of the rotation period P of Kepler-78.
templanets are shownas greensquares. Otherwell-characterized rot
exoplanets are plotted as black points. Exoplanet masses, radii, Thismadeseparationofthesignalsrelatedtothestellar
andtheirassociatederrorsarefromtheExoplanetOrbitDatabase rotation and the signal related to the planetary orbital
(http://exoplanets.org; downloaded on 23 August 2014). Planets periodpossible. Iftheplanetaryperiodwasalargerfrac-
withfractionalmassuncertaintiesofover50%arenotshown.
tion or a harmonic of the stellar rotation period, decon-
volving the signal due to the star and the signal due to
the previous studies. The spread in the Doppler ampli-
the planet would be much more difficult.
tudes reported by Hatzes (2014) illustrates that slight
differences in the choice of noise model has a signifi- The true benefit of this analysis technique comes from
cant effect on the planetary signal extracted, and that its nonparametric nature. This analysis is particularly
the error on the planetary mass measurement has likely powerful in that even if the actual nature of the noise
been underestimated. Thus, by starting with the sim- being modeled is unclear, a GP model can still be used
plestpossiblenonparametricdescriptionsofthedata,we to explore the nature of the noise and identify the most
find a description of the data that is minimally complex evident components. Due to the success of extracting a
yet strongly supported by Bayesian information criteria. quasi-periodic signal at the rotation period of the star
Furthermore, the structure of the best fit GP kernel as in independent datasets, we conclude that the predom-
well as the success of sharing kernel parameters to de- inant noise signal in the RV datasets likely comes from
scribe both datasets supports our belief that we have the stellar activity of Kepler-78. The variability of this
robustly measured the astrophysical variability of the noise, due to the growth and decay of active regions on
10 Grunblatt et al.
thestellarsurface,cannotbefullyparametrizedwiththe exponentialsignalscanbemodeledwiththeGP,indicat-
information we have. Thus, testing the nonparametric ingthatitcouldbeparticularlyusefulindescribingother
Gaussian process noise model is a valuable exploration noise modes seen in radial velocity studies as well as the
afterearlierparametricandnonparametricmodelsofthe already well-understood stellar surface signals. We hope
RVactivity. Furthermore,theuseoftheKepler photom- to usethis technique ona systemwith contemporaneous
etry as a prior on our estimate of the RV activity makes photometryandspectroscopytoexploretherelationship
our Gaussian process analysis especially valuable when between photometric and spectroscopic signals of exo-
photometric data is more readily available for a system planetary systems.
than spectroscopic data, as is often the case.
We plan to further test the robustness of this tech-
nique by analyzing other RV datasets of exoplanetary We thank the anonymous referee for their suggestions.
systems. This method is particularly useful for analyz- WealsothankSuzanneAigrain,B.J.Fulton,ConorMc-
ing stellar RV datasets over long time baselines which Partland, and Maxwell Service for helpful discussions.
cannot be modeled easily because of the spot evolution A.W.H. acknowledges NASA grant NNX12AJ23G. We
on a timescale less than an order of magnitude larger gratefully acknowledge the efforts and dedication of the
than the stellar rotation period (as determined by the Keck Observatory staff and extend special thanks to
autocorrelation function analysis done by P13). Such those of Hawai‘ian ancestry on whose sacred mountain
scenarios can be easily described with a quasi-periodic of Mauna Kea we are privileged to be guests. With-
GP.Inaddition, sincekernelfunctionscanbecombined, outtheirgeneroushospitality,theKeckobservationspre-
anysortofphysicalcombinationofperiodicandlinearor sented herein would not have been possible.
REFERENCES
Bumba,V.1963,BulletinoftheAstronomicalInstitutesof Hillenbrand,L.,Isaacson,H.,Marcy,G.,etal.2015,in
Czechoslovakia,14,91 CambridgeWorkshoponCoolStars,StellarSystems,andthe
Desort,M.,Lagrange,A.-M.,Galland,F.,Udry,S.,&Mayor,M. Sun,Vol.18,CambridgeWorkshoponCoolStars,Stellar
2007,A&A,473,983 Systems,andtheSun,ed.G.T.vanBelle&H.C.Harris,
Dumusque,X.,Santos,N.C.,Udry,S.,Lovis,C.,&Bonfils,X. 759–766
2011,A&A,527,A82 Howardetal.,A.W.2013,Nature,503,381
Foreman-Mackey,D.,Hogg,D.W.,Lang,D.,&Goodman,J. Knutson,H.A.,Charbonneau,D.,Allen,L.E.,Burrows,A.,&
2013,PASP,125,306 Megeath,S.T.2008,ApJ,673,526
Fortney,J.J.,Marley,M.S.,&Barnes,J.W.2007,ApJ,659, McDonough,W.,&Sun,S.-S.1995,ChemicalGeology,120,223
1661 Meunier,N.,Desort,M.,&Lagrange,A.-M.2010,A&A,512,A39
Gelman,A.,&Rubin,D.B.1992,StatisticalScience,7,457 Neal,R.M.1997,arXivpreprintphysics/9701026
Gibson,N.P.,Aigrain,S.,Roberts,S.,etal.2012,MNRAS,419, Pepe,F.,Cameron,A.C.,Latham,D.W.,etal.2013,Nature,
2683 503,377
Gibson,N.P.,Pont,F.,&Aigrain,S.2011,MNRAS,411,2199 Quin˜onero-Candela,J.,&Rasmussen,C.E.2005,TheJournalof
Hatzes,A.P.2014,A&A,568,A84 MachineLearningResearch,6,1939
Haywood,R.D.,CollierCameron,A.,Queloz,D.,etal.2014, Rasmussen,C.E.,&Williams,C.K.I.2006,GaussianProcesses
MNRAS,443,2517 forMachineLearning(AdaptiveComputationandMachine
Herbrich,R.,Lawrence,N.D.,&Seeger,M.2003,inAdvancesin Learning)(TheMITPress)
NeuralInformationProcessingSystems15,ed.S.Becker, Sanchis-Ojeda,R.,Rappaport,S.,Winn,J.N.,etal.2013,The
S.Thrun,&K.Obermayer(MITPress),625–632 AstrophysicalJournal,774,54
Swain,M.R.,Vasisht,G.,&Tinetti,G.2008,Nature,452,329
Wang,J.M.,Fleet,D.J.,&Hertzmann,A.2008,Pattern
AnalysisandMachineIntelligence,IEEETransactionson,30,
283