Additional Keplerian Signals in the HARPS data for Gliese 667C from a Bayesian Re-analysis
aa r X i v : . [ a s t r o - ph . S R ] D ec Mon. Not. R. Astron. Soc. , 1– ?? (2002) Printed 18 April 2018 (MN L A TEX style file v2.2)
Additional Keplerian Signals in the HARPS data forGliese 667C from a Bayesian Re-analysis
Philip C. Gregory ⋆ Physics and Astronomy Department, University of British Columbia, 6224 Agricultural Rd., Vancouver, BC V6T 1Z1, Canada
Submitted to MNRAS Dec. 16, 2012
ABSTRACT
A re-analysis of Gliese 667C HARPS precision radial velocity data was carried out witha Bayesian multi-planet Kepler periodogram (from 0 to 7 planets) based on a fusionMarkov chain Monte Carlo algorithm. The most probable number of signals detectedis 6 with a Bayesian false alarm probability of 0.012. The residuals are shown to beconsistent with white noise. The 6 signals detected include two previously reportedwith periods of 7.198 (b) and 28.14 (c) days, plus additional periods of 30.82 (d), 38.82(e), 53.22, and 91.3 (f) days. The 53 day signal is probably the second harmonic ofthe stellar rotation period and is likely the result of surface activity. The existence ofthe additonal Keplerian signals suggest the possibilty of further planets, two of which(d and e) could join Gl 667Cc in the central region of the habitable zone. N-bodysimulations are required to determine which of these signals are consistent with astable planetary system. M sin i values corresponding to signals b, c, d, e, and f are ∼ ⊕ , respectively. Key words: stars: planetary systems; methods: statistical; methods: data analysis;techniques: radial velocities.
The HARPS spectrograph mounted on the ESO/3.6-m LaSilla telescope is yielding a rich bounty of information re-garding planets around M dwarfs (eg., Bonfils et al. 2011).There is a lot of interest in searches for planets about lowmass M dwarfs. Firstly, a planet of given mass and orbitalseparation induces a larger stellar radial velocity variationaround a lower mass star. Secondly, the low luminosity ofM dwarfs moves their habitable zone much nearer to thestar. For these reasons a habitable zone (HZ) planet arounda 0.3 M ⊙ M dwarf produces a 7 times larger radial velocitywobble than the same planet orbiting a solar-mass G dwarf(Delfosse et al. 2012).One M dwarf of particular interest is Gliese 667C (Gl667C, also known as GJ 667C), an isolated component of ahierarchical triple system. The two others components, Gl667AB, are a closer couple of K dwarfs. Gl 667AB has asemi-major axis of 1.82 AU (period of 42.15 years) and a to-tal mass of 1.27 M ⊙ (dynamically determined, Soderhjelm1999). Gl 667C is the lightest component with a mass of0 . ± .
019 M ⊙ (Anglada-Escud´e et al. 2012a). It is at aprojected distance of 32.4” from Gl 667AB, giving an ex-pected semi-major axis of ∼
300 AU (for a distance of 7.23 ⋆ E-mail: [email protected] pc and a factor of 1.26 between expected and projected semi-major axis, Fischer and Marcy 1992).In 2011, Bonfils et al. (2011) reported a super-earth Gl667Cb ( M sin i = 5 . ⊕ ) with a period of 7.2 d andevidence for two other interesting periods of 90 and 28d. The possibility of a 28 d period planet was particu-larly interesting because it would fall in the HZ. Two re-cent papers (Anglada-Escud´e et al. 2012a & Delfosse et al.2012), have confirmed planet Gl 667Cb and a 28 d pe-riod planet Gl 667Cc ( M sin i = 4 . ⊕ ) in the HZ. TheAnglada-Escud´e et al. (2012a) results are not fully inde-pendent as they depend on the 143 Bonfils et al. (2011)RV observations which they reduced using the HARPS-TERRA (Template Enhanced Radial velocity Re-analysisApplication) algorithm (Anglada-Escud´e et al. 2012b), sup-ported by observations from two other spectrographs. TheDelfosse et al. (2012) data includes an additional 36 HARPSradial velocity measurements.The excitement generated by this and many other exo-planetary discoveries has spurred a significant effort to im-prove the statistical tools for analyzing data in this field(e.g., Loredo & Chernoff 2003, Loredo 2004, Cumming 2004,Gregory 2005a & b, Ford 2005 & 2006, Ford & Gregory2007, Tuomi & Kotiranta 2009, Dawson & Fabrycky 2010,Cumming & Dragomir 2010). Much of this work has high-lighted a Bayesian MCMC approach as a way to better un- c (cid:13) P. C. Gregory derstand parameter uncertainties and degeneracies and tocompute model probabilities.Gregory (2011, 2011b, 2012) developed a Bayesian fu-sion MCMC algorithm that incorporates parallel tempering(PT), simulated annealing and a genetic crossover operationto facilitate the detection of a global minimum in χ . Thisenables the efficient exploration of a large model parameterspace starting from a random location. When implementedwith a multi-planet Kepler model , it is able to identifyany significant periodic signal component in the data thatsatisfies Kepler’s laws and is able to function as a multi-planet Kepler periodogram . The fusion MCMC algorithmincludes an innovative adaptive control system that auto-mates the selection of efficient parameter proposal distribu-tions even if the parameters are highly correlated. One appli-cation of the algorithm (Gregory & Fischer 2010) confirmedthe existence of a disputed second planet (Fischer et al.2002) in 47 Ursae Majoris (47 UMa) and provided orbitalconstraints on a possible additional long period planet witha period ∼ The adaptive fusion MCMC (FMCMC) is a very generalBayesian nonlinear model fitting program. After specify-ing the model, M i , the data, D , and priors, I , Bayes theo-rem dictates the target joint probability distribution for themodel parameters which is given by p ( ~X | D, M i , I ) = C p ( ~X | M i , I ) × p ( D | M i , ~X, I ) . (1)where C is the normalization constant and ~X represent theset of model parameters. The first term on the RHS of theequation, p ( ~X | M i , I ), is the prior probability distribution of ~X , prior to the consideration of the current data D . Thesecond term, p ( D | ~X, M i , I ), is called the likelihood and it isthe probability that we would have obtained the measureddata D for this particular choice of parameter vector ~X ,model M i , and prior information I . At the very least, theprior information, I , must specify the class of alternativemodels (hypotheses) being considered (hypothesis space ofinterest) and the relationship between the models and thedata (how to compute the likelihood). In some simple casesthe log of the likelihood is simply proportional to the familiar For multiple planet models, there is no analytic expression forthe exact radial velocity perturbation. In many cases, the radialvelocity perturbation can be well modeled as the sum of multipleindependent Keplerian orbits which is what has been assumed inthis paper. Following on from the pioneering work on Bayesian peri-odograms by Jaynes (1987) and Bretthorst (1988) χ statistic. For further details of the likelihood function forthis type of problem see Gregory (2005b).To compute the marginals for any subset of theparameters it is necessary to integrate the joint probabilitydistribution over the remaining parameters. For example,the marginal probability distribution (a probability densityfunction) of the orbital period in a one planet radial velocitymodel fit is given by p ( P | D, M , I ) = Z dK Z dV Z de Z dχ × Z dω Z dslope Z ds × p ( P, K, V, e, χ, ω, slope, s | D, M , I ) ∝ p ( P | M , I ) Z dK · · · Z ds × p ( K, V, e, χ, ω, slope, s | M , I ) × p ( D | M , P, K, V, e, χ, ω, slope, s, I ) , (2)where P, K, V, e, χ, ω are the radial velocity model parame-ters and s is an extra noise parameter which is discussed inSection 3. slope is a parameter that accounts for the longperiod orbital motion induced in Gl 667C by Gl 667AB. p ( P | M , I ) is the prior for the orbital period parameter, p ( K, V, e, χ, ω, slope, s | M , I ) is the joint prior for the otherparameters, and p ( D | M , P, K, V, e, χ, ω, slope, s, I ) is thelikelihood. For a seven planet model fit we need to inte-grate over 37 parameters to obtain the marginal probabilitydistribution for each of the seven period parameters. Inte-gration is more difficult than maximization, however, theBayesian solution provides the most accurate informationabout the parameter errors and correlations without theneed for any additional calculations, i.e., Monte Carlo sim-ulations. Bayesian model selection requires integrating overall the model parameters.In high dimensions, the principle tool for carrying outthe integrals is Markov chain Monte Carlo based on theMetropolis algorithm. The greater efficiency of an MCMCstems from its ability, after an initial burn-in period, to gen-erate samples in parameter space in direct proportion to thejoint target probability distribution. In contrast, straightMonte Carlo integration randomly samples the parameterspace and wastes most of its time sampling regions of verylow probability.MCMC algorithms avoid the requirement for com-pletely independent samples, by constructing a kind of ran-dom walk in the model parameter space such that the num-ber of samples in a particular region of this space is pro-portional to a target posterior density for that region. Therandom walk is accomplished using a Markov chain, wherebythe new sample, ~X t +1 , depends on previous sample ~X t ac-cording to a time independent entity called the transitionkernel , p ( ~X t +1 | ~X t ). The remarkable property of p ( ~X t +1 | ~X t )is that after an initial burn-in period (which is discarded)it generates samples of ~X with a probability density pro-portional to the desired posterior p ( ~X | D, M , I ) (e.g., seeChapter 12 of Gregory (2005a) for details). c (cid:13) , 1– ?? Bayesian Re-analysis of the Gliese 667C HARPS Data The joint posterior probability distribution in model pa-rameter space is typically highly multi-modal for exoplanetradial velocity (RV) analysis. The Metropolis algorithm canbecome stuck in the vicinity of a local probability maxi-mum. To avoid this fusion MCMC (FMCMC) incorporatesthree other algorithms each of which is designed to facilitatethe detection of a global minimum in χ (or a maximum inprobability). They are parallel tempering (PT) (Geyer 1991and re-invented by Hukushima & Nemoto 1996), simulatedannealing and a genetic crossover operation. All three areimplement in each FMCMC run. The combination greatlyfacilitates the identification of a global maximum in prob-ability. This was made possible through the developmentof an adaptive control system which has been described indetail most recently by Gregory (2011b, 2012). Further re-finements of the control system are ongoing.At each iteration of the FMCMC, a single joint pro-posal is made to jump to a new location in parameter spacefrom the current location. The key to an efficient MCMCis an efficient method of proposing new jumps especiallywhen there are correlations between the parameters. Thisis further complicated in PT because multiple MCMC areexecuted in parallel each exploring a different probabilitydistribution. In FMCMC this difficult step is automated bythe control system saving a great deal of time and effort. In this section we describe the model fitting equations andthe selection of priors for the model parameters. We haveinvestigated the Gl 667C data using models ranging from0 to 7 planets. For a one planet model the predicted radialvelocity is given by v ( t i ) = V + K [cos { θ ( t i + χP ) + ω } + e cos ω ] , (3)and involves the 6 unknown parameters V = a constant velocity. K = velocity semi-amplitude. P = the orbital period. e = the orbital eccentricity. ω = the longitude of periastron. χ = the fraction of an orbit, prior to the start of datataking, that periastron occurred at. Thus, χP = the numberof days prior to t i = 0 that the star was at periastron, foran orbital period of P days. θ ( t i + χP ) = the true anomaly, the angle of the star inits orbit relative to periastron at time t i .We utilize this form of the equation because we obtainthe dependence of θ on t i by solving the conservation ofangular momentum equation dθdt − π [1 + e cos θ ( t i + χ P )] P (1 − e ) / = 0 . (4)Our algorithm is implemented in Mathematica and it provesfaster for
Mathematica to solve this differential equationthan solve the equations relating the true anomaly to themean anomaly via the eccentric anomaly.
Mathematica gen-erates an accurate interpolating function between t and θ so the differential equation does not need to be solved sepa-rately for each t i . Evaluating the interpolating function for each t i is very fast compared to solving the differential equa-tion. Details on how equation 4 is implemented and the ac-curacy of the interpolation as a function of eccentricity aregiven in the Appendix A of Gregory (2011b).As described in more detail in Gregory 2007, we em-ployed a re-parameterization of χ and ω to improve theMCMC convergence speed motivated by the work of Ford(2006). The two new parameters are ψ = 2 πχ + ω and φ = 2 πχ − ω . Parameter ψ is well determined for all ec-centricities. Although φ is not well determined for low ec-centricities, it is at least orthogonal to the ψ parameter.We use a uniform prior for ψ in the interval 0 to 4 π anduniform prior for φ in the interval − π to +2 π . This in-sures that a prior that is wraparound continuous in ( χ, ω )maps into a wraparound continuous distribution in ( ψ, φ ).To account for the Jacobian of this re-parameterization itis necessary to multiply the Bayesian integrals by a factorof (4 π ) − nplan , where nplan = the number of planets in themodel. In this work we utilized the orthogonal combination( ψ, φ ) as well as the FMCMC correlated proposal schemedescribed in Gregory (2012).In a Bayesian analysis we need to specify a suitable priorfor each parameter. These are tabulated in Table 1. For thecurrent problem, the prior given in equation 2 is the productof the individual parameter priors. Detailed arguments forthe choice of each prior were given in Gregory (2007) andGregory & Fischer (2010).All of the models considered in this paper incorporatean extra noise parameter, s , that can allow for any addi-tional noise beyond the known measurement uncertainties.We adopt an independent Gaussian distribution with a vari-ance s where is becomes another parameter in the modelfit. Thus, the combination of the known errors and extranoise has a Gaussian distribution with variance = σ i + s ,where σ i is the known measurement uncertainty for i th datapoint. In general, nature is more complicated than our modeland known measurement errors. Marginalizing s has the de-sirable effect of treating anything in the data that can’t beexplained by the model and known measurement errors asnoise, leading to conservative estimates of orbital parame-ters . See Sections 9.2.3 and 9.2.4 of Gregory (2005a) fora tutorial demonstration of this point. If there is no extranoise then the posterior probability distribution for s willpeak at s = 0. The upper limit on s was set equal to K max .We employed a modified scale invariant prior for s with aknee, s = 1 m s − . There is an additional benefit for incorporating the extra noiseparameter s . When the χ of the fit is very large, the BayesianMarkov chain automatically inflates s to include anything in thedata that cannot be accounted for by the model with the currentset of parameters and the known measurement errors. This resultsin a smoothing out of the detailed structure in the χ surfaceand allows the Markov chain to explore the large scale structurein parameter space more quickly. This is similar to simulatedannealing, but does not require choosing a cooling scheme. Thechain begins to decrease the value of the extra noise as it settlesin near the best-fit parameters. This is demonstrated in Fig. 2 ofGregory (2011b) and in Section 2.1 of Gregory (2012).c (cid:13) , 1– ?? P. C. Gregory
Table 1.
Prior parameter probability distributions.Parameter Prior Lower bound Upper boundOrbital frequency p (ln f , ln f , · · · ln f n | M n , I ) = n ![ln( f H /f L )] n a ( n =number of planets)Velocity K i Modified scale invariant b = 1) K max (cid:0) P min P i (cid:1) / p − e i (m s − ) ( K + K ) − ln (cid:20) K max K (cid:0) P min Pi (cid:1) / p − e i (cid:21) K max = 2129V (m s − ) Uniform − K max K max e i Eccentricity Ecc. noise bias correction filter 0 0.99 χ orbit fraction Uniform 0 1 ω i Longitude of Uniform 0 2 π periastron slope Uniform c − s Extra noise (m s − ) ( s + s ) − ln (cid:0) s max s (cid:1) = 1) K max a Gl 667C is part of a triple star system with Gl 667AB, a much closer binary. We adopted at an upperlimit of 95 yr by setting the gravitational pull on the planet due to Gl 667AB = 1% of the pull from Gl667C. b Since the prior lower limits for K and s include zero, we used a modified scale invariant prior of the form p ( X | M, I ) = 1 X + X (cid:0) X max X (cid:1) (5)For X ≪ X , p ( X | M, I ) behaves like a uniform prior and for X ≫ X it behaves like a scale invariant prior.The ln (cid:0) X max X (cid:1) term in the denominator ensures that the prior is normalized in the interval 0 to X max . c Since this parameter is common to all models including the 0 planet model, the exact upper and lowerbounds are not critical for model selection. The range simply needs to be large enough so as not to effectthe parameter estimates.
In Gregory & Fischer (2010), the velocities of model fitresiduals were randomized in multiple trials and processedusing a one planet version of the FMCMC Kepler peri-odogram. In this situation periodogram probability peaksare purely the result of the effective noise. The orbits corre-sponding to these noise induced periodogram peaks exhib-ited a well defined statistical bias towards high eccentric-ity. We offered the following explanation for this effect. Tomimic a circular velocity orbit the noise points need to becorrelated over a larger fraction of the orbit than they doto mimic a highly eccentric orbit. For this reason it is morelikely that noise will give rise to spurious highly eccentricorbits than low eccentricity orbits. The bias found using multiple sets of randomized residuals froma 5 planet fit to 55 Cancri combined Lick and Keck data agreedclosely with the bias found for multiple sets of randomized resid-uals from both 2 and 3 planet fits to 47 UMa Lick data.
Gregory & Fischer (2010) characterized this eccentric-ity bias and designed a correction filter that can be used asan alternate prior for eccentricity to enhance the detectionof planetary orbits of low or moderate eccentricity. On thebasis of our understanding of the mechanism underlying theeccentricity bias, we expect the eccentricity prior filter tobe generally applicable to searches for low amplitude orbitalsignals in precision radial velocity data sets. The probabilitydensity function for this filter is given by pdf ( e ) = 1 . − . e + 0 . e − . e − . . (6)Fig. 11 of Gregory (2011b) demonstrates that the effectof this noise eccentricity bias correction filter on the finalmarginal eccentricity distributions for the low and moder-ate eccentricity orbits of Gl 581b, c, & d is insignificant.In a related study, Shen and Turner (2008) exploredleast- χ Keplerian fits to synthetic radial velocity data sets.They found that the best fit eccentricities for low signal-to-noise ratio
K/σ c (cid:13) , 1– ?? Bayesian Re-analysis of the Gliese 667C HARPS Data - -
500 0 500 1000 - - - H - L V e l o c i t y H m s - L Figure 1.
The HARPS data for Gl 667C after correcting for thesecular acceleration. N obs
60, were systematically biased to higher values, lead-ing to a suppression of the number of nearly circular orbits.More recently, Zakamska et al. (2011) found that eccentric-ities of planets on nearly circular orbits are preferentiallyoverestimated, with typical bias of one to two times the me-dian eccentricity uncertainty in a survey, e.g., 0.04 in theButler et al. catalogue (Butler et al. 2006). When perform-ing population analysis, they recommend using the mode ofthe marginalized posterior eccentricity distribution to min-imize potential biases.In the analysis of the Gl 667C data we used the eccen-tricity noise bias correction filter as the eccentricity prior onfits of all the models.
For all the analysis we used the HARPS data given byDelfosse et al. (2012). We subtracted a mean velocity of6.5474477 km s − and converted to units of m s − . Follow-ing Delfosse et al. (2012) we also subtracted a componentdue to the secular acceleration of 0.21 m s − yr − arisingfrom the proper motion. This small correction resulted in achange in radial velocity over the span of the data amount-ing to 1.527 m s − . Fig. 1 shows the corrected HARPS datafor Gl 667C which exhibits a large positive slope indicativeof a period much longer than the data. A likely explana-tion is the orbital motion of Gl 667C relative to the centerof mass of the Gl 667ABC system. The Gl 667ABC sys-tem parameters suggested an orbital period of ∼ ∼ − and a maximum rateof change in radial velocity of ∼ − yr − . Multiplyingby the data duration of 7.27 yr yields a maximum expectedvelocity change of 20 m s − , which is comparable to whatis observed. Our zero reference time is the mean time ofthe HARPS observations which corresponds to a Julian daynumber = 2 , , . slope parameter in all ourmodels.The models considered range from a zero planet modelto a seven planet model. As mentioned all models include a - Iterations H ‰ L Log @ P r i o r ‰ L i ke D Iterations H ‰ L P e r i od s Figure 2.
The upper panel is a plot of the Log [Prior × Like-lihood] versus iteration for a 2 planet Kepler periodogram of theHARPS data. The lower shows the values of the two unknownperiod parameters versus iteration number. The two starting pe-riods of 7.2 and 15 d are shown on the left hand side of the plotat a negative iteration number. slope parameter and extra noise ( s ) parameter. The currentsection deals with model parameter estimation while Sec-tion 5 deals with model selection. The next 6 subsectionsshow the FMCMC based Kepler periodograms for the 2 - 7planet cases. We don’t show Kepler periodogram for the oneplanet model as its parameters are well understood. Five 2 planet Kepler periodograms were computed. In allcases signals were detected at 7.2 and 28 d but the 28 dsignal was never the dominant second peak. Other periodswhich occured in individual periodograms included 106, 184,249, 383 d. Fig. 2 shows the results for the periodogramwhich achieved the highest peak Log [Prior × Likelihood].The upper plot shows the Log [Prior × Likelihood] versusiteration. The lower plot shows the two period parametersversus iteration which shows a stable 7.2 d signal while thesecond period transitions over many peaks thanks to theparallel tempering feature of the FMCMC algorithm. Therelative peak probabilities of the second period options isillustrated in Fig. 3 which shows a plot of the two periodparameter values versus a normalized value of Log [Prior × Likelihood]. Fig. 4 shows a plot of the eccentricity parametervalues versus period. Only the 7.2 and 28.1 d periods areconsistent with low eccentricity orbits. The median value ofthe extra noise parameter s = 1 . − . Six 3 planet Kepler periodograms were computed for thedata. In all cases signals were detected at 7.2 and 28.1 d. c (cid:13) , 1– ?? P. C. Gregory P = d P = d - - - - - Period H d L Log @ P r i o r ‰ L i ke D Figure 3.
A plot of the 2 period FMCMC parameter samplesversus a normalized value of Log [Prior × Likelihood] for the 2planet Kepler periodogram. P = d Periods E cce n t r i c i t y Figure 4.
A plot of eccentricity versus period for the FMCMCparameter samples from the 2 planet Kepler periodogram.
Other periods which occured in individual periodogramswere 38.8, 106, 184, 368 d. In three cases the third periodwas 184 d. Fig. 5 shows a plot of eccentricity versus pe-riod for the periodogram which achieved the highest peakLog [Prior × Likelihood]. The 7.2 and 28.1 d signals areconsistent with low eccentricity orbits. The third period in-variably had a high eccentricity as is the case for the 184d period in the figure. The median value of the extra noiseparameter s = 1 . − . Periods E cce n t r i c i t y Figure 5.
A plot of eccentricity versus period for the FMCMCparameter samples from the 3 planet Kepler periodogram.
Periods E cce n t r i c i t y Figure 6.
A plot of eccentricity versus period for the FMCMCparameter samples from the 4 planet Kepler periodogram.. P = . d P = . d - - - - - - - Periods
Log @ P r i o r ‰ L i ke D Figure 7.
A plot of the five period parameter values versus anormalized value of Log [Prior × Likelihood] for the 5 planetKepler periodogram.
Three 4 planet Kepler periodograms were computed for thedata. In all cases signals were detected at 7.2, 28.1, 106, &190 days. Additional periods that were detected in differentruns were 38.8 and 91d. Fig. 6 shows a plot of eccentricityversus period for the periodogram which achieved the high-est peak Log [Prior × Likelihood]. The median value of theextra noise parameter s = 1 . − . Four 5 planet Kepler periodograms were computed for thedata. In all cases signals were detected at 7.2, 28.1, 53.2, & 91d. Other periods which occured in individual periodogramsincluded 30.8, 38.8, 87, 106, 128, 184 d. In some cases thefifth period exhibited a variety of periods. Fig. 7 shows theperiodogram which achieved the highest peak Log [Prior × Likelihood]. In this case the extra period exhibits a dominantpeak at 86.6 d and a weaker peak at 30.8 d. Fig. 8 is a plot ofeccentricity versus period. The 86.6 d signal exhibits a higheccentricity while the 30.8 d signal has a low eccentricity.The median value of the extra noise parameter s = 0 .
79 ms − . c (cid:13) , 1– ?? Bayesian Re-analysis of the Gliese 667C HARPS Data P = . d P = . d Period H d L E cce n t r i c i t y Figure 8.
A plot of eccentricity versus period for the FMCMCparameter samples from the 5 planet Kepler periodogram..
We next computed a 6 planet Kepler periodogram analysisof the data. An initial blind search with starting periodsof 5, 10, 15, 20, 50, 100 d yielded 5 well defined periodsat 7.2, 28.1, 30.8, 53.2, 91.3 d. The remaining period wassplit between a high eccentricity 14.4 d period and a loweccentricity 35 d period. A second run starting from the 5well defined periods plus the 35 d low eccentricity option isshown in Fig. 9. After a temporary stay at the 35 d period,the black trace transitions to a stable 38. 8 d peak . Asdiscussed in Section 5 on model selection, the Bayes factorfavors the 6 planet model by a factor of 137 compared tothe next leading contender. The median value of the extranoise parameter s = 0 .
54 m s − .Fig. 9 shows a plot of the 6 period parameter values (in-cluding burn-in points to show the 35 d signal) versus a nor-malized value of Log [Prior × Likelihood] for the 6 planetKepler periodogram. The fourth period, shown in black, ex-hibits two peaks, a weak one at 35 d and the other at 38.8d. The 35 d period coincides with a one year alias of thestronger 38.8 d period. The spectral window function for theHARPS data exhibits two peaks, 1 d and one year. We cangain further insight into the relationship of the 35 and 38.8d periods from a 6 planet Kepler periodogram of a subset ofthe HARPS data, i.e., the first 143 data points. Again dur-ing the burn-in, the fourth period locks on to the 35 d peakbefore transitioning to the 38.8 d peak. Fig. 11 shows the6 period parameter values (including burn-in points) versusa normalized value of Log [Prior × Likelihood] for the 143d sub-sample. In this case, the alias at 35 d actually has ahigher peak probability density by a factor of ∼
10, eventhough the much larger number of samples for the 38.8 dpeak indicates that there is more probable associated withthe 38.8 d peak. A narrow peak in the joint parameter spaceof high probability density can contain much less total prob-ability than a broader region of lower probability density.Note that the other one year alias of the 38.8 d period at Another 6 planet Kepler periodogram was computed startingfrom the best 6 planet parameters set with the exception that53.2 d period which was replaced by twice this period which cor-responds roughly with the assumed rotation period of the star of105 d (Delfosse et al. 2012, based on a periodogram of a stellaractivity diagnostic). This solution with periods of 7.2, 28.1, 30.8,38.8, 91 and 106.5 d had a peak probability 2600 times lower.
Iterations H ‰ L Log @ P r i o r ‰ L i ke D Iterations H ‰ L P e r i od s Figure 9.
The upper panel is a plot of the Log [Prior × Like-lihood] versus iteration for the six planet fit of the HARPS data.The lower shows the values of the six unknown period parametersversus iteration number. The six starting periods are shown onthe left hand side of the plot at a negative iteration number. . d d
3. 7.2 28.1 38.8 53.2 91.3 300 - - - - - Period H d L Log H P r i o r ‰ L i ke li hood L Figure 10.
Full HARPS data set: A plot of the 6 period parame-ter values versus a normalized value of Log [Prior × Likelihood]for the 6 planet Kepler periodogram (includes burn-in points toshow the 35 d signal). The fourth period, shown in black, exhibitstwo peaks one at 35 d and the other at 38.8 d. The 35 d periodis a one year alias of the stronger 38.8 d period. ∼ .It was not until the 6 planet model that the 30.8 and38.8 d signals became the dominant third and fourth periods.Both are consistent with low eccentricity orbits as shown inFig. 12 which is a plot of eccentricity versus period for thepost-burn-in samples for the full HARPS data set.The closeness of the 28.1 and 30.8 day periods raisesquestions about a possible relationship between these twosignals. Could the 30.8 d signal be an alias of the 28.1 d c (cid:13) , 1– ?? P. C. Gregory d d . d
3. 7.2 16 28 38.8 54 91 - - - - - Period H d L Log H P r i o r ‰ L i ke li hood L Figure 11.
Subset of HARPS data (first 143 points): A plotof the 6 period parameter values versus a normalized value ofLog [Prior × Likelihood] for the 6 planet Kepler periodogram(includes burn-in points to show the 35 d signal). The fourthperiod, shown in black, exhibits two peaks one at 35 d and theother at 38.8 d. For the partial data set the 35 d alias has a higherpeak probability density. The other one year alias at 43.4 d is justdiscernible. . d
3. 7.2 28.1 38.8 53.2 91.30.00.20.40.60.81.0
Periods E cce n t r i c i t y Figure 12.
Full HARPS data set: A plot of eccentricity versusperiod for the 6 planet Kepler periodogram for the post-burn-inFMCMC samples. signal? An alias results from a convolution of the spectralwindow function of the data sample times with the spec-trum of the real signal. If the real signal is removed then analias of that signal should not be present in the residuals.In our case the multi-planet model is requiring both to bepresent. Nevertheless, we transformed the marginal distri-bution of the second period, P (28.14 d signal), by the rela-tionship P alias = 1 / (1 /P − / .
25) to obtain the pridictedmarginal distribution the one year alias. This is comparedit to the marginal for P (30.82 d signal) in the top panelof Fig. 13. Clearly, the predicted alias (shown dashed) doesnot overlap the marginal for P . For comparison, the bottompanel shows the predicted marginal distribution of the lowerone year alias (dashed) of the 38.8 d signal to the marginaldistribution obtained for the 35 d signal (solid) which wasobserved during the burn-in phase of the full HARPS dataset. In this case the distributions coincide.Fig. 14 shows a plot of a subset of the FMCMC param-eter marginal distributions for the 6 planet FMCMC fit ofthe data. The eccentricities of the 4 lowest periods are con-sistent with low eccentricity orbits while the 53.2 and 91.3d periods show peak eccentricities of 0.41 and 0.36, respec- P H d L P D F P H d L P D F Figure 13.
The upper panel shows a comparison of the predictedmarginal distribution (dashed) of the closest one year alias of the28.14 d signal to the marginal distribution for the 30.82 d signal(solid). For comparison, the bottom panel shows the predictedmarginal distribution of the lower one year alias (dashed) of the38.8 d signal to the marginal distribution obtained for the 35 dsignal (solid) which was observed during the burn-in phase of thefull HARPS data set. tively. The median value of extra noise parameter s = 0 .
54m s − .Phase plots for the 6 planet model are shown in Fig. 15.The top left panel shows the data and model fit versus 7.2d period phase after removing the effects of the five otherorbital periods plus V and slope parameter. The FMCMCoutput for each iteration is a vector of the 6 planet orbitalparameter set plus V , slope , and the extra noise parameters.The 7.2 d period phase plot is constructed from a sample ofFMCMC iterations (typically 300) and for each iteration wecompute the predicted velocity points for that realizationof the 5 planet plus V and slope parameters. We then con-struct the mean of these model predictions and subtract themean prediction from the data points. These residuals forthe set of observation times are converted to residuals ver-sus phase using the mode of the marginal distribution for the7.2 d period parameter. A period phase model velocity fitis then computed at 100 phase points for each realization ofthe 7.2 d planet parameter set obtained from the same sam-ple of 300 FMCMC iterations. At each of these 100 phasepoints we construct the mean model velocity fit and mean ± ± slope parameter after removing the effects of the 6 orbitalperiods plus V .Fig. 16 shows a generalized Lomb-Scargle (GLS) peri- c (cid:13) , 1– ?? Bayesian Re-analysis of the Gliese 667C HARPS Data P H d L P D F e P D F K H m s - L P D F P H d L P D F e P D F
2. 2.50.00.51.01.52.0 K H m s - L P D F P H d L P D F e P D F
1. 2.0.00.51.01.52.0 K H m s - L P D F P H d L P D F e P D F
1. 1.50.00.51.01.52.0 K H m s - L P D F P H d L P D F e P D F
1. 2.0.00.51.01.52.0 K H m s - L P D F P H d L P D F e P D F K H m s - L P D F H m (cid:144) s (cid:144) yr L P D F H m s - L P D F Figure 14.
A plot of a subset of the FMCMC parameter marginaldistributions for the 6 planet Kepler periodogram. odogram (Zechmeister & K¨urster 2009) for the maximum a posteriori (MAP) parameter values of the 6 planet fitresiduals. The GLS allows for a floating offset and weights.The dashed horizontal lines correspond to peak periodogramlevels for which the false alarm probability (FAP) would= 0 . . . There is no evidence of any significant peaksor red noise in these residuals.Another way of exploring the residuals is to computethe autocorrelation function, ρ ( j ). Fig. 17 shows ρ ( j ) of theresiduals for the 6, 3 and 1 planet fits computed from equa-tion (7). ρ ( j ) = P overlap [( x i − x ) ( x i + j − x )] qP overlap ( x i − x ) × qP overlap ( x i + j − x ) , (7)where x i is the i th residual, j is the lag and x is the meanof the samples in the overlap region. Because the data arenot uniformly sampled, for each lag all sample pairs thatdiffered in time by this lag ± . The interesting region of power is where the frequentist p-valueis small ( ≪ ≈ M ∗ (1 − z ) ( N − / (Zechmeister & K¨urster 2009), where z = maximum periodogram power, M = the number of indepen-dent frequencies and N is the number of data points. Clearly,the FAP value for the actual highest peak in the periodogramis much larger than 0.1. Cumming (2004) recommends setting M = ∆ f/δf , where ∆ f is the frequency range examined ≈ f max ,and δf = the resolution of the periodogram ≈ /T where T isthe duration of the data. P = - P Orbital phase V e l o c i t y H m s - L P = - P Orbital phase V e l o c i t y H m s - L P = - - - P Orbital phase V e l o c i t y H m s - L P = - P Orbital phase V e l o c i t y H m s - L P = - - - P Orbital phase V e l o c i t y H m s - L P = - P Orbital phase V e l o c i t y H m s - L Slope = (cid:144) s (cid:144) yr - -
500 0 500 1000 - - H d L V e l o c i t y H m s - L Figure 15.
The top left panel shows the data and model fit versus7.2 d period phase after removing the effects of the 5 other orbitalperiods plus V and slope parameters. The upper and lower curvesare the mean FMCMC model fit ± slope parameter after removing the effectsof the 6 orbital periods plus V . FAP = = frequency H (cid:144) d L P o w e r Figure 16.
A periodogram of the 6 planet fit residuals. right panel of Fig. 17 shows a plot of the number of suchsample pairs as a function of the lag. Clearly for large lagsthe uncertainty in computed ρ ( j ) is expected to be larger.Clearly, the 6 planet model autocorrelation function isconsistent with white noise. The thick solid line in the 1planet residuals panel is the average autocorrelation func-tion generated from 400 simulated data sets of a 5 planetmodel (28.1, 30.8, 38.8, 53.2, and 91.4 d periods) togetherwith the measurement errors. The 5 planet model param-eter values are the MAP values derived from the 6 planetKepler periodogram of the real data. The good agreementbetween the simulation and actual ρ ( j ) argues that in thecase of Gl 667C any colored noise evident in the 1 planet fit c (cid:13) , 1– ?? P. C. Gregory ææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ = - - Lag H d L A u t o c o rr e l a t i on ææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ = - - Lag H d L A u t o c o rr e l a t i on æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ ææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ = - - Lag H d L A u t o c o rr e l a t i on æææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ Lag H d L N u m b e r o f po i n t s Figure 17.
The top panel shows the autocorrelation function, ρ ( j ), of the HARPS 6 planet fit residuals computed as explainedin the text. The next two panels show ρ ( j ) for the 3 planet resid-uals and the 1 planet residuals (together with a simulation thicksolid line). The bottom panel shows the number of HARPS samplepairs available for computing the autocorrelation function versuslag. is nicely accounted for by additional signals (mostly planets)not included in the model. In spite of the absence of any obvious signal from the 6planet residual periodogram, we computed three 7 planetKepler periodogram to see if there were any new surprises. . d Periods E cce n t r i c i t y Figure 18.
A plot of eccentricity versus period for the 7 planetKepler periodogram of the full HARPS data.
Each trial used a starting period set of 7.2, 28.1, 30.8, 35.2,53.2, 91.3 d and a different choice for the seventh period.Encouragingly, all trials recovered the 6 periods found inthe 6 planet fit which indicates they are stable features.Fig. 18 shows the plot of eccentricity versus period for thetrial which achieved the highest peak Log [Prior × Likeli-hood]. The additional 3.5 d period is a weak high eccentric-ity feature typical of noise (Gregory & Fischer 2010). Themedian value of the extra noise parameter s = 0 .
23 m s − . One of the great strengths of Bayesian analysis is the built-in Occam’s razor. More complicated models contain largernumbers of parameters and thus incur a larger Occampenalty, which is automatically incorporated in a Bayesianmodel selection analysis in a quantitative fashion (see forexample, Gregory (2005a), p. 45). The analysis yields therelative probability of each of the models explored.To compare the posterior probability of the i th planetmodel to the 4 planet model we need to evaluate the oddsratio, O i = p ( M i | D, I ) /p ( M | D, I ), the ratio of the poste-rior probability of model M i to model M . Application ofBayes theorem leads to, O i = p ( M i | I ) p ( M | I ) p ( D | M i , I ) p ( D | M , I ) ≡ p ( M i | I ) p ( M | I ) B i (8)where the first factor is the prior odds ratio, and the secondfactor is called the Bayes factor , B i . The Bayes factor isthe ratio of the marginal (global) likelihoods of the models.The marginal likelihood for model M i is given by p ( D | M i , I ) = Z d ~Xp ( ~X | M i , I ) × p ( D | ~X, M i , I ) . (9)Thus Bayesian model selection relies on the ratio of marginallikelihoods, not maximum likelihoods. The marginal likeli-hood is the weighted average of the conditional likelihood,weighted by the prior probability distribution of the modelparameters.The marginal likelihood can be expressed as the prod-uct of the maximum likelihood and the Occam penalty (seeGregory (2005a), page 48). The Bayes factor will favor themore complicated model only if the maximum likelihood ra-tio is large enough to overcome this penalty. In the simple c (cid:13) , 1– ?? Bayesian Re-analysis of the Gliese 667C HARPS Data case of a single parameter with a uniform prior of width∆ X , and a centrally peaked likelihood function with char-acteristic width δX , the Occam factor is ≈ δX/ ∆ X . If thedata is useful then generally δX ≪ ∆ X . For a model with m parameters, each parameter will contribute a term to theoverall Occam penalty. The Occam penalty depends not onlyon the number of parameters but also on the prior range of each parameter (prior to the current data set, D ), as sym-bolized in this simplified discussion by ∆ X . If two modelshave some parameters in common then the prior ranges forthese parameters will cancel in the calculation of the Bayesfactor.To make good use of Bayesian model selection, we needto fully specify priors that are independent of the currentdata D . In most instances we are not particularly interestedin the Occam factor itself, but only in the relative probabil-ities of the competing models as expressed by the Bayes fac-tors. Because the Occam factor arises automatically in themarginalization procedure, its effect will be present in anymodel selection calculation. Note: no Occam factors arise inparameter estimation problems. Parameter estimation canbe viewed as model selection where the competing modelshave the same complexity so the Occam penalties are iden-tical and cancel out.The MCMC algorithm produces samples which are inproportion to the posterior probability distribution whichis fine for parameter estimation but one needs the pro-portionality constant for estimating the model marginallikelihood. Clyde et al. (2007) reviewed the state of tech-niques for model selection from a statistical perspective andFord & Gregory (2007) have evaluated the performance ofa variety of marginal likelihood estimators in the exoplanetcontext.Estimating the marginal likelihood is a very big chal-lenge for models with large numbers of parameters, e.g., ourseven planet model has 38 parameters. In this work we em-ploy the nested restricted Monte Carlo (NRMC) methodintroduced in Gregory & Fischer (2010) and described inmore detail in Gregory (2012) to estimate the marginal likeli-hoods. Monte Carlo (MC) integration can be very inefficientin exploring the whole prior parameter range because it ran-domly samples the whole volume. The fraction of the priorvolume of parameter space containing significant probabilityrapidly declines as the number of dimensions increase. Forexample, if the fractional volume with significant probabil-ity is 0.1 in one dimension then in 38 dimensions the frac-tion might be of order 10 − . In restricted MC integration(RMC) this is much less of a problem because the volumeof parameter space sampled is greatly restricted to a regiondelineated by the outer borders of the marginal distributionsof the parameters for the particular model.In RMC, most of the random samples occur close tothe outer borders of the restricted region because the con-tribution to the volume of parameter space is greatest there.In NRMC integration, multiple boundaries are constructedbased on credible regions ranging from 30% to > The more surprising the result the stronger the evidence re-quired to overcome our skepticism. contributions. For example, for the interval between the 30%and 60% credible regions, we generate random parametersamples within the 60% region and reject any sample thatfalls within the 30% region. Using the remaining samples wecan compute the contribution to the NRMC integral fromthat interval.The left panel of Fig. 19 shows the contributions fromthe individual intervals for 5 repeats of the NRMC eval-uation for the 6 planet model. Note the large range inparameter volume on the abscissa. The right panel showsthe summation of the individual contributions versus thevolume of the credible region. The credible region listedas 9995% is defined as follows. Let X U and X L corre-spond to the upper and lower boundaries of the 99% cred-ible region, respectively, for any of the parameters. Simi-larly, X U and X L are the upper and lower boundaries ofthe 95% credible region for the parameter. Then X U = X U + ( X U − X U ) and X L = X L + ( X L − X L ).Similarly, X U = X U + ( X U − X U ). For each credi-ble region interval approximately 320,000 MC samples wereused.The mean value of the prior × likelihood within the30% credible region is a factor of 4 . × larger than themean in the shell between the 97 and 99% credible regions.However, the volume of parameter space in the shell betweenthe 97 and 99% credible regions is a factor of 5 . × larger than the volume within the 30% credible region sothe contribution from the latter to the marginal likelihoodis negligible. For further details on the NRMC method seeGregory (2012).The NRMC method is expected to underestimate themarginal likelihood in high dimensions and this underes-timate is expected to become worse the larger the num-ber of model parameters, i.e. increasing number of plan-ets (Gregory 2012). When we conclude, as we do, that theNRMC computed odds in favor of the six planet model com-pared to the four planet model is 137 we mean that the trueodds is > m planets the FAP m is the probability that there areactually fewer than m planets, i.e., m − m = m − X i =0 (prob . of i planets) (10)If we assume a priori that all models under considera-tion are equally likely, then the probability of each model isrelated to the Bayes factors by p ( M i | D, I ) = B i P Nj =0 B j (11)where N is the maximum number of planets in the hypoth-esis space under consideration, and of course B = 1. Forthe purpose of computing FAP m we set N = m . Substitut-ing Bayes factors from Table 2 into equation 10 givesFAP = ( B + B + B + B + B + B ) P j =0 B j = 0 . c (cid:13) , 1– ?? P. C. Gregory
Table 2.
Marginal likelihood estimates, Bayes factors relative to model 4, and false alarm probabilities. The last two columns list themedian of extra noise parameter, s , and the RMS residual. The MAP value of s appears below in parentheses.Model Periods Marginal Bayes factor False Alarm s RMS residual(d) Likelihood nominal Probability (m s − ) (m s − ) M . × − . × − M (7 .
2) (8 . +0 . − . ) × − . × − . × − M (7 . , . +0 . − . ) × − . × − . × − M (7 . , . , . +0 . − . ) × − . × − . × − M (7 . , . , , . × . × . ) × − . . × − M (7 . , . , , . ,
91) (3 . × . × . ) × − .
64 0 .
61 0.79 1.6(0.70) M (7 . , . , . , . , ,
91) (7 . × . × . ) × −
137 0 .
012 0.54 1.5(0.08) M (3 . , . , . , . , . , ,
91) (5 . × . × . ) × − .
96 0 .
993 0.23 1.4(0.03) à à à à à à à à à à àà àæ æ æ æ æ æ ææ æ æ ææ æì ì ì ì ì ì ìì ì ì ìì ìà à à à à à à à à à àà àæ æ æ æ æ æ ææ æ æ ææ æ % % % % % % % % % % % % % - - - - - - - - - - - @ Restricted Monte Carlo parameter volume D Log @ D M a r g i n a l L i ke li hood D à à à à à à à à à à àà àæ æ æ æ æ æ ææ æ æ ææ æì ì ì ì ì ì ìì ì ì ìì ìà à à à à à à à à à àà àæ æ æ æ æ æ ææ æ æ ææ æ % % % % % % % % % % % % % - - - - - - - - - - - @ Restricted Monte Carlo parameter volume D Log @ M a r g i n a l L i ke li hood D Figure 19.
Left panel shows the contribution of the individual nested intervals to the NRMC marginal likelihood for the 6 planet model.The right panel shows the integral of these contributions versus the parameter volume of the credible region. (12)For the 6 planet model we obtain a low FAP ≈ − .Table 2 gives the NMRC Marginal likelihood estimates,Bayes factors and false alarm probabilities for 0, 1, 2, 3, 4,5, 6, and 7 planet models which are designated M , · · · , M .The last two columns list the median estimate of the extranoise parameter, s (MAP value in parentheses below), andthe RMS residual. For each model the NRMC calculationwas repeated 5 times and the quoted errors give the spreadin the results, not the standard deviation. The Bayes factorsthat appear in the third column are all calculated relativeto model 4.A summary of the 6 planet model parameters and theiruncertainties are given in Table 3. The quoted value is themedian of the marginal probability distribution for the pa- rameter in question and the error bars identify the bound-aries of the 68.3% credible region . The value immediatelybelow in parenthesis is the MAP estimate, the value at themaximum of the joint posterior probability distribution. Forthe eccentricity parameter we also include the mode withindouble parentheses. It is not uncommon for the MAP esti-mate to fall close to the borders of the credible region. In one In practice, the probability density for any parameter X is rep-resented by a finite list of values p i representing the probability indiscrete intervals δX . A simple way to compute the 68.3% credi-ble region, in the case of a marginal with a single peak, is to sortthe p i values in descending order and then sum the values untilthey approximate 68.3%, keeping track of the upper and lowerboundaries of this region as the summation proceeds.c (cid:13) , 1– ?? Bayesian Re-analysis of the Gliese 667C HARPS Data % credible region
90 91 92 930.00.51.01.52.02.53.03.5 P H d L P D F Figure 20.
The marginal distribution of the 91.26 d signal (solid)is compared to the marginal distribution of the 30.82 d signal(dashed) after multiplying its period scale by a factor of 3. Thedashed and solid vertical lines indicate the 99% credible regionsof the two signals. case, the eccentricity of the fourth planet, the MAP estimatefalls just outside the 68.3% credible region which is one rea-son why we prefer to quote median or mode values as well.The semi-major axis and M sin i values are derived from themodel parameters assuming a stellar mass of 0 . ± . ⊙ (Anglada-Escud´e et al. 2012a). The quoted errors onthe semi-major axis and M sin i include the uncertainty inthe stellar mass. We first consider whether any of the 6 signals detected isthe result of stellar activity. Delfosse et al. (2012) analyzedseveral stellar activity diagnostics and found a high peak inone diagnostic (the full-width at half-maximum (FWHM)of the crosscorrelation function) with a period of ∼
105 d,which they interpret as the rotation period of the star. Inour Kepler periodgrams a 106 d period was detected startingat the 2 planet model which transitioned to a 53 d periodcommencing with the 5 planet model. Clearly these two pe-riods are harmonically related. If we assume that the 106 dperiod is the star rotation period then depending on the con-figuration of surface activity (eg., spots) it would not be toosurprising to detect RV variations at the second harmonicof the rotation frequency. This is a possible interpretationof the 53 d period feature in our analysis.In Sec. 4.5, we established that the 30.82 d period is notan alias of the 28.14 d period. Fig. 20 demonstrates that the30.82 is not a harmonic of the 91.26 d period, both of whichshow up in the 6 plan and 7 planet Kepler periodograms.As another test of the stability of the 6 planet modelresults, we subtracted one of the signals at a time (for the91, 53, and 28 d signals) from the data and carried out a5 planet Kepler periodogram of the modified data startingfrom an initial period set of 5, 10, 15, 20, and 25 d. In eachcase, we recovered the other 5 periods. Fig. 21 shows the5 planet periodogram result of the data with the 53 d sig-nal subtracted. The upper panel is a plot of the Log [Prior × Likelihood] versus iteration for a 5 planet Kepler peri-odogram. The middle panel shows the values of the five pe-riod parameters versus iteration number. The five startingperiods of 5, 10, 15, 20, 25 d are shown on the left hand side - - - Iterations H ‰ L Log @ P r i o r ‰ L i ke D Iterations H ‰ L P e r i od s P = . d Period H d L E cce n t r i c i t y Figure 21.
The upper panel is a plot of the Log [Prior × Like-lihood] versus iteration for a 5 planet Kepler periodogram of theHARPS data with 53 d signal subtracted. The middle shows thevalues of the five period parameters versus iteration number. Thefive starting periods of 5, 10, 15, 20, 25 d are shown on the lefthand side of the plot at a negative iteration number. The bottompanel is a plot of the eccentricity parameter versus period for postburn-in iterations. of the plot at a negative iteration number. The bottom panelis a plot of the eccentricity versus period for post burn-initerations. The marginal distributions for the 7.2, 28.1, 30.8,38.8, & 91 d signals are indistinguishable from those foundfor the 6 planet Kepler periodogram of the original data.This example also demonstrates the ability of the FMCMCbased Kepler periodogram to carry out a blind search in28 parameters, including 5 period parameters which have ahuge prior range. In fact, detection of the six signals foundin Sec. 4.5 would be a difficult feat for any analysis methodthat relies on a conventional periodogram of the residualsof an n signal fit to estimate the starting period for the n + 1 signal search. The ability of the FMCMC algorithmto explore a multi-modal environment has proven to be veryimportant in this analysis.The prospect exists that digging deeper with a morepowerful statistical algorithm might also uncover newKeplerian-like spectral artifacts of the radial velocity mea-surement system or of the star’s surface activity. There- c (cid:13) , 1– ?? P. C. Gregory
Table 3.
Six planet model parameter estimates. The value immediately below in parenthesis is the MAPestimate. For the eccentricity parameter, the value within double parentheses is the mode.Parameter planet 1 planet 2 planet 3 planet 4 planet 5 planet 6 P (d) 7 . + . − . . + . − . . + . − . . +0 . − . . +0 . − . . +0 . − . (7.1972) (28.120) (30.82) (38.87) (53.23) (91.37) K (m s − ) 3 . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . (3.95) (2.32) (1.45) (1.33) (1.55) (1.76) e . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . (0.068) (0.123) (0.22) (0.58) (0.48) (0.34)((0.068)) ((0.028)) ((0.115)) ((0.48)) ((0.45)) ((0.35)) ω (deg) 342 +35 − +68 − +48 − +31 − +20 − +20 − (338) (266) (220) (297) (296) (217) a (AU) 0 . + . − . . + . − . . + . − . . + . − . . +0 . − . . +0 . − . (0.0494) (0.123) (0.130) (0.152) (0.187) (0.269) M sin i (M ⊕ ) 5 . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . (5.44) (5.01) (3.17) (2.62) (3.67) (5.34)Periastron 4495 . +0 . − . +5 − +4 − +3 − +2 − +5 − passage (4503.0) (4499) (4503) (4480) (4499) (4477)(JD - 2,450,000) fore, N-body simulations will be required to determine whichcombination of the signals are consistent with a stable plan-etary system.In the rest of this discusion we speculate on the pos-sibility that some or all of the remaining 5 signals de-tected are planetary in origin. One surprise is the close-ness of the 28 d and 30.8 d orbits whose semi-major axesdiffer by only 0.007 AU. Closest approach would occur ev-ery 1 / (1 / . − / . M = 4 . . ⊕ , Carter et al.2012) with a minimum separation of 0.014 AU, Kepler 42b& d ( M . . ⊕ , Muirhead et al. 2012) witha minimum separation of 0.0038 AU, and KOI-55b & c( M ∼ .
44 & ∼ .
66 M ⊕ , Charpinet et al. 2012) with aminimum separation of 0.0016 AU. The masses of the par-ent star for Kepler 36, Kepler 42 and KOI-55 are 1.07, 0.13and 0.496 M ⊙ , respectively. Kepler 36b & c both show tran-sit timing variations consistent with a gravitational interac-tion between the two planets while no such timing variationswere detected for the three Kepler 42 planets. By definition, the habitable zone around a host star is theregion where liquid water can be stable on the surface of arocky planet (Huang 1959, Kasting et al. 1993). In the cur-rent absence of observational constraints, Selsis et al. (2007)chose to assess the habitable zone for planets with as few as-sumptions as possible on their physical and chemical nature.Only two conditions were assumed although they may derive from complex geophysical properties. Below we summarizethese assumptions and some of the key concepts.1) The amount of superficial water must be large enough sothat the surface can host liquid water for any temperaturebetween the temperature at the triple point of water, 273 K,and the critical temperature of water, 647 K. This conditionimplies that the water reservoir produces a surface pressurehigher than 220 bars when fully vaporized. With an Earthgravity, this corresponds to a 2.2 km layer of water, slightlylower than the mean depth of Earth oceans of 2.7 km. For agravity twice that of Earth, this pressure corresponds to halfthis depth. Planets with less water may still be habitable,but their HZ may be somewhat narrower because liquid wa-ter would disappear at a lower surface temperature.2) Atmospheric CO must accumulate in a planet’s atmo-sphere whenever the mean surface temperature falls belowthe freezing point of water. Walker et al. (1981) proposedthat Earth’s long-term climate was buffered by a negativefeedback mechanism involving atmospheric CO levels andthe dependence of silicate weathering rates on climate. Inthe carbonate-silicate cycle, the CO emitted by volcanoesis dissolved in rain water forming carbonic acid (H CO )which weathers silicate rocks. The products of silicate weath-ering, calcium (Ca ++ ), bicarbonate (HCO − ) ions and dis-solved silica (SiO ), are transported by streams and riversto the ocean. There, organisms, such as foraminifera, di-atoms and radiolarians, use the products to make shells ofcalcium carbonate (CaCO ). Most of the shells redissolve,but a fraction of them survive and are buried in sedimentson the seafloor. The combination of silicate weathering pluscarbonate precipitation can be represented chemically byCO + CaSiO → CaCO + SiO . The seabed is eventu-ally subducted where the heat pushes the carbonate-silicate c (cid:13) , 1– ?? Bayesian Re-analysis of the Gliese 667C HARPS Data cycle reversible reaction in the opposite direction, eventuallyreleasing CO through volcanism.This cycle replenishes all the CO in the combinedatmosphere-ocean system on a timescale of approximatelyhalf a million years. It is thus too slow to counteract human-induced global warming, but fast enough to have a dominat-ing effect on the billion-year timescale of planetary evolution(Kasting & Catling 2003). Weathering rates increase bothbecause of the direct effect of temperature on chemical re-action rates and because evaporation (and, hence, precipi-tation) rates increase as the surface temperature increases.If we increase the orbital radius, the surface temperaturefalls and the weathering decreases allowing atmospheric CO to build up leading to a stabilization of the temperaturethrough the greenhouse effect. This negative feedback in thecarbonate-silicate cycle stabilizes the long-term surface tem-perature and the amount of CO in the atmosphere of theEarth. Such an assumption implies that the planet is ge-ologically active and continuously outgassing CO . It alsoimplies that carbonates form in the presence of surface liq-uid water, which may require continental weathering.At a sufficiently large orbital radius the planet gets coldenough that water vapor disappears followed by the freezingof CO which results in the permanent loss of both green-house gases (positive feedback). With no atmospheric CO ,or with a fixed CO level as in Hart (1979), the HZ couldbe ∼
10 times narrower. In the absence of a greenhouse gaslike CO and H O, the present Earth would be frozen.With decreasing orbital radius, a positive feedbackeventually takes over in which too much water accumulatein the atmosphere causing the temperature to increase. At asufficiently high temperature water stops condensing. Withno rain the weathering ceases and CO builds up, furtherincreasing the temperature. Eventually the oceans evapo-rate completely. Water in the high atmosphere is dissociatedinto hydrogen and oxygen, the hydrogen escapes to space,the oxygen combines with minerals and all the water disap-pears.According to Selsis et al. (2007), planets with massesoutside the 0.5 - 10 M ⊕ range cannot host liquid surfacewater. Planets under the lower end of this range have tooweak a gravity to retain a sufficiently dense atmosphere, andthose above the upper end accrete a massive He-H envelope.In either case, the pressure at the surface is incompatiblewith liquid water. The 10 M ⊕ upper limit is somewhat fuzzy,since planets in the 3 - 10 M ⊕ range can have very differentdensities.Based on the above two major assumptions,(Selsis et al. 2007) derives boundaries for the HZ asshow in Fig. 22 for a range of stellar masses. Assuminga planetary origin for the Keplerian-like signals reportedhere, the 28.1, 30.8, and 38.8 d signals (labeled c, d, and e)translate to orbits in the centre of the HZ. The semi-majoraxis of a 91.3 d planet (f) would lie just within the outer-most extreme edge of the HZ, although its eccentric orbitwill result in it spending the majority of its time outsidethe HZ. The white object in the diagram corresponds to theorbital parameters implied by a 53 d Keplerian signal. Asexplained above, it is likely that the star’s surface activityis responsible for this signal as it is the second harmonic ofwhat is believed to be the star’s rotation period. Figure 22.
The darker area is the orbital region that remainscontinuously habitable during at least 5 Gyr as a function of thestellar mass (Selsis et al. 2007). The light grey region gives thetheoretical inner (runaway greenhouse) and outer limits with 50%cloudiness, with H O and CO clouds, respectively. The dottedboundaries correspond to the extreme theoretical limits, foundwith a 100% cloud cover. The dashed line indicates the distanceat which a 1 M ⊕ planet on a circular orbit becomes tidally lockedin less than 1 Gyr. A Bayesian re-analysis of published HARPS precision radialvelocity data Delfosse et al. (2012) for Gl 667C was carriedout with a multi-planet Kepler periodogram (from 0 to 7planets) based on our fusion Markov chain Monte Carlo al-gorithm. In all cases the analysis included an unknown pa-rameterized stellar jitter noise term and an unknown slope parameter to deal with the linear drift resulting from theorbital motion of Gl 667C relative to the center of mass ofthe Gl 667ABC system. The most probable number of sig-nals detected is 6 with a Bayesian false alarm probability of0.012. The residuals are shown to be consistent with whitenoise. The Keplerian signals detected include the 7.2 and28.1 d (planets b and c) periods reported previously plusperiods of 30.8 (d), 38.8 (e), 53.2, and 91.3 (f) d. The 53.2d period appears to correspond to the second harmonic ofthe stellar rotation period and is likely the result of surfaceactivity (eg., spots). The same set of periods were also de-tected in a subset of the data consisting of the first 143 datapoints.The existence of the additonal Keplerian-like signalssuggest the possibilty of further planets, two of which (dand e) could join Gl 667Cc in the central region of the hab-itable zone. M sin i values corresponding to signals b, c, d,e, and f are ∼ ⊕ , respectively.It is also possible that digging deeper with a more powerfulstatistical algorithm might have uncovered new spectral ar-tifacts of the radial velocity measurement system or of thestar’s surface activity. N-body simulations are required todetermine which of these signals are consistent with a sta-ble planetary system. c (cid:13) , 1– ?? P. C. Gregory
ACKNOWLEDGMENTS
The author would like to thank Wolfram Research for pro-viding a complementary license to run gridMathematica tofacilitate parallel processing.
REFERENCES
Anglada-Escud´e, G., Arriagada, P., Vogt, S. S., et al. 2012a,ApJ, 751, L16Anglada-Escud´e, G. & Butler, R. P. 2012b, ApJS, 200, 15Bonfils, X., Delfosse, X., Udry, S., et al., 2012,arXiv:1111.5019v2Bretthorst, G. L., 1988, Bayesian Spectrum Analysis andParameter Estimation, New York: Springer-VerlagButler, R. P., Wright, J. T., Marcy, G. W., Fischer, D. A.,Vogt, S. S., Tinney, C. G., Jones, H. R. A., Carter, B.D., Johnson, J. A., McCarthy, C., and Penny, A. J., 2006,ApJ, 646, 505Charpinet, S., Fontaine, G., Brassard, P., Green, E. M.,Van Grootel, V., Randall, S. K., Silvotti, R., Baran, A.S., stensen, R. H., Kawaler, S. D., Telting, J. H., 2011,Nature, 480, 496Carter, J. A.,Algol, E., Chaplin, W. J., Basu, S., Bedding,T. R., Buchhave, L. A., Christensen-Dalsgaard, J., Deck,K. M., Elsworth, Y., Fabrycky, D. C., Ford, E. B., Fortney,J. J., Hale, S. J., Handberg, R., Hekker, S., Holman, M.J., Huber, D., Karoff, C., Kawaler, S. D., Kjeldsen, H.,Lissauer, J. J., Lopez, E. D., Lund, M. N., Lundkvist,M., Metcalfe, T. S., Miglio, A., Rogers, L. A., Stello, D.,Borucki, W. J., Bryson, S., Christiansen, J. L., Cochran,W. D., Geary, J. C., Gilliland, R. L., Haas, M. R., Hall, J.,Howard, A. W., Jenkins, J. M., Klaus, T., Koch, D. G.,Latham, D. W., MacQueen, P. J., Sasselov, D., Steffen, J.H., Twicken, J. D., Winn, J. N., 2012, Science, 337, 556Cumming, A., 2004, MNRAS, 354, 1165Cumming, A., Dragomir, D., 2010, MNRAS, 401, 1029Clyde, M. A., Berger, J. O., Bullard, F., Ford, E. B., Jef-freys, W. H., Luo, R., Paulo, R., Loredo, T., 2007, in ‘Sta-tistical Challenges in Modern Astronomy IV,’ G. J. Babuand E. D. Feigelson (eds.), ASP Conf. Ser., 371, 224Dawson, R. I, and Fabrycky, D. C., 2010, ApJ, 722, 937Delfosse, X., Bonfils, X., Forveille, T., Udry, S., Mayor, M.,Bouchy, F., Gillon, F. M., Lovis, C., V. Neves, V., Pepe,F., Perrier, C., Queloz, D., Santos, N. C., and S´egransan,D. 2012, arXiv:1202.2467v1Fischer, D. A. & Marcy, G. W. 1992, ApJ, 396, 178Fischer,D. A., Marcy, G. W., Butler, R. P., Laughlin, G.L., and Vogt, S. S., 2002, ApJ, 564, 1028Ford, E. B., 2005, AJ, 129, 1706Ford, E. B., 2006, ApJ, 620, 481Ford, E. B., & Gregory, P. C., 2007, in ‘Statistical Chal-lenges in Modern Astronomy IV,’ G. J. Babu and E. D.Feigelson (eds.), ASP Conf. Ser., 371, 189Geyer, C. J., 1991, in E. M. Keramidas (ed.),‘ComputingScience and Statistics: 23rd symposium on the interface,Interface Foundation, Fairfax Station, pp. 156-163Gregory, P. C., 2005a, ‘Bayesian Logical Data Analysisfor the Physical Sciences: A Comparative approach with
Mathematica
Support’, Cambridge University PressGregory, P. C., 2005b, ApJ, 631, 1198 Gregory, P. C., 2007, MNRAS, 374, 1321Gregory, P. C., and Fischer, D. A., 2010, MNRAS, 403, 731Gregory, P. C., 2011, MNRAS, 410, 94Gregory, P. C., 2011b, MNRAS, 415, 2523Gregory, P. C., Chapter 7 in Astrostatistical Challengesfor the New Astronomy, Springer Series in Astrostatistics,Hilbe, J.M (ed), 2012, New York:SpringerHart, M. H. 1979, Icarus, 37, 351Huang, S. S. 1959, PASP, 71, 421Hukushima, K., and Nemoto, K., 1996, Journal of the Phys-ical Society of Japan, 65(4), 1604Jaynes, E.T., 1987, ‘Bayesian Spectrum & Chirp Analysis,’in
Maximum Entropy and Bayesian Spectral Analysis andEstimation Problems , ed. C.R. Smith and G.L. Erickson,D. Reidel, Dordrecht, p. 1Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993,Icarus, 101, 108Kasting, J. F., & Catling, D. (2003), Ann. Rev. of Astron-omy and Astrophysics, 41, 429Loredo, T., 2004, in ‘Bayesian Inference And MaximumEntropy Methods in Science and Engineering: 23rd Inter-national Workshop’, G.J. Erickson & Y. Zhai, eds, AIPConf. Proc. 707, 330 (astro-ph/0409386)Loredo, T. L. and Chernoff, D., 2003, in ‘Statistical Chal-lenges in Modern Astronomy III’, E. D. Feigelson and G.J. Babu (eds) , Springer, New York, p. 57Muirhead, P. S., Johnson, J. A., Apps, K., Carter, J. A.,Morton, T. D., Fabrycky, D. C., Pineda, J. S., Bottom,M., Rojas-Ayala, B., Schlawin, E., Hamren, K., Covey, K.R., Crepp, J. R., Stassun, K. G., Pepper, J., Hebb, L.,Kirby, E. N., Howard, A. W., Isaacson, H. T., Marcy, G.W., Levitan, D., Diaz-Santos, T., Armus, L., Lloyd, J. P.,2012, ApJ, 747, 144Selsis, F., Kasting, J. F., Levrard, B., et al. 2007, A&A,476, 1373Schneider, J., Le Sidaner, P., Savalle, R., Zolotukhin, I,2012, Astron. Soc. of the Pac. Conference Series, 461, 447Shen,Y., and Turner, E. L., ApJ, 685, 553Soderhjelm, S. 1999, A&A, 341, 121Tuomi, M. & Kotiranta, S. 2009, A&A, 496, L13Walker, J. C. G., Hays, P. B., & Kasting, J. F. 1981, J.Geophys. Res., 86, 9776Zakamska, N. L., Pan, M., Ford, E. B., 2011, MNRAS, 410,1895Zechmeister, M. & K¨urster, M., 2009, A&A, 496, 577 c (cid:13) , 1–, 1–