Issue 
A&A
Volume 528, April 2011



Article Number  L5  
Number of page(s)  5  
Section  Letters  
DOI  https://doi.org/10.1051/00046361/201015995  
Published online  24 February 2011 
Letters to the Editor
Bayesian reanalysis of the radial velocities of Gliese 581
Evidence in favour of only four planetary companions
^{1}
University of Hertfordshire, Centre for Astrophysics Research, Science and
Technology Research Institute, College Lane,
AL10 9AB
Hatfield
UK
^{2}
University of Turku, Tuorla Observatory, Department of Physics and
Astronomy, Väisäläntie
20, 21500
Piikkiö,
Finland
email: mikko.tuomi@utu.fi
Received: 25 October 2010
Accepted: 16 February 2011
Aims. The Gliese 581 planetary system has received attention because it has been proposed to host a lowmass planet in its habitable zone. We reanalyse the radial velocity measurements reported to contain six planetary signals to see whether these conclusions remain valid when the analyses are made using Bayesian tools instead of the common periodogram analyses.
Methods. We analyse the combined radial velocity data set obtained using the HARPS and HIRES spectrographs using posterior sampling techniques and computation of the posterior probabilities of models with differing numbers of Keplerian signals. We do not fix the orbital eccentricities and stellar jitter to certain values but treat these as free parameters of our statistical models. Hence, we can take the uncertainties of these parameters into account when assessing the number of planetary signals present in the data, the point estimates of all of the model parameters, and the uncertainties of these parameters.
Results. We conclude that based on the Bayesian model probabilities and the nature of the posterior densities of the different models, there is evidence in favour of four planets orbiting GJ 581. The HARPS and HIRES data do not imply the conclusion that there are two additional companions orbiting GJ 581. We also revise the orbital parameters of the four companions in the system. Especially, according to our results, the eccentricities of all the companions in the system are consistent with zero.
Key words: planets and satellites: detection / methods: statistical / techniques: radial velocities / stars: individual: Gliese 581
© ESO, 2011
1. Introduction
The nearby M dwarf Gliese 581 has received plenty of attention during recent years. In 2005 a Neptunemass planet candidate was found in its orbit using the HARPS spectrograph (Bonfils et al. 2005). Two years later it was reported to be a host to two additional superEarths (Udry et al. 2007). In 2009 an new Earthmass planet with a minimum mass of 1.9 M_{⊕} was found in its orbit, making it a system with four planetary companions of relatively low mass (Mayor et al. 2009).
Ever since the discovery of the first companion orbiting GJ 581, the star has been a target of intensive radial velocity (RV) surveys because few M dwarf stars are known to be hosts to planetary systems despite the fact that they are numerous even in the Solar neighbourhood. As a result, in 2010, GJ 581 was reported to have two more companions of planetary mass with GJ 581 g, a 3.1 M_{⊕} planet, in the habitable zone of the star (Vogt et al. 2010).
The discovery of this 6planet system was made by analysing two highprecision RV data sets made using the HARPS spectrograph (Mayor et al. 2009) and the HIRES spectrograph (Vogt et al. 2010). The combined data set consists of 241 RV measurements. In Vogt et al. (2010), the Keplerian signals of the six planets were discovered by studying the periodogram of the combined timeseries and by fitting the orbital parameters of the proposed companions to the data.
The purpose of this Letter is to see whether Bayesian data analysis gives consistent results with those reported by Vogt et al. (2010). Our major concern is, that by fixing eccentricities to zero the uncertainties of these eccentricities and their effect on the detectability of planetary signals were not taken into account by Vogt et al. (2010). They did let the eccentricities float freely and concluded that this did not produce any significan improvement to the fit. However, they did not discuss whether the uncertainty about the eccentricities could have an effect on the probabilities of finding the Keplerian signals in the data in the first place. Also, in Vogt et al. (2010), the stellar jitter was estimated to have a value of 1.4 m s^{1} – simply because it yielded a reduced χ^{2} value of unity. Our second concern is that the uncertainty of the jitter was not taken into account either in the analyses. If its value was under or overestimated, it could have a significant effect on the detectability of the planetary signals as demonstrated by e.g. Ford (2006); Gregory (2007a,b); Tuomi & Kotiranta (2009).
In this Letter, we reanalyse the combined RV data set of GJ 581 using Bayesian tools – posterior samplings and model probabilities. First, we sample the parameter spaces of Keplerian models with k planetary companions by letting k = 0,...,6. Second, we calculate the Bayesian probabilities for each of these models ż_{k}. Finally, we compare our results with those of Vogt et al. (2010) to see whether their periodogram analyses and fitting algorithms and Bayesian ones do yield similar results in this case.
2. Modelling and Bayesian model comparison
2.1. Statistical model and posterior samplings
Following Tuomi & Kotiranta (2009), we assume that the planets do not interact with one another in the timescale of the measurements and model the superposition of k Keplerian signals simply by summing their effect on the RV. Consequently, there are five parameters describing the signature of an individual planet: RV amplitude (K), orbital period (P), orbital eccentricity (e), longitude of pericentre (ω), and the mean anomaly (M_{0}). As suggested by Ford (2006), to improve the efficiency of the sampling of the parameter space, we use the logarithm of the period in the samplings because it is a scaleinvariant parameter.
Our statistical model consists of the sum of Keplerian signals and two sources of uncertainty. These sources are the instrument noise and the noise caused by the stellar surface – the stellar jitter. We model these as independent random variables with Gaussian density and zero mean. We assume that the standard deviation of the instrument noise is known and use the values reported in Mayor et al. (2009); Vogt et al. (2010). However, we adopt the standard deviation of the stellar jitter as a free parameter of our statistical model and denote it as σ_{J}. The statistical model can be written as (1)where r_{i,l} is the measurement at time t_{i} made using telescopeinstrument combination l, ż_{k} represents the k Keplerian signals, parameter γ_{l} is the reference velocity, and ϵ_{i} and ϵ_{J} are Gaussian random variables describing the instrument noise, as reported by the observers, and the stellar jitter, respectively. As a result, there are 5k + 3 free parameters in the parameter vectors θ_{k} of our models – five parameters for each planet, jitter magnitude, and the reference velocities of the HARPS and HIRES measurements.
We sample the parameter spaces of the different models using the adaptive Metropolis algorithm of Haario et al. (2001). This sampling algorithm is similar to the famous MetropolisHastings algorithm (Metropolis et al. 1953; Hastings 1970) but it adapts to the information gathered during the first n − 1 samples from the parameter posterior density by approximating this sample as a multivariate Gaussian density. The nth sample is then drawn by using this multivariate Gaussian as a proposal density with the n − 1th value as a mean. While only an approximation, this algorithm converges to the posterior relatively rapidly – apparently even in the case of multimodal posterior sampled in this Letter.
When calculating the posterior densities for the parameters representing semimajor axes and RV masses of the planets, we took into account the uncertainty of the stellar mass. This mass is estimated to be 0.31 ± 0.02 M_{⊙} for GJ 581 (Delfosse et al. 2000). We sampled the densities of the semimajor axes and RV masses by drawing random numbers from the estimated density of the stellar mass that had a mean of 0.31 M_{⊙} and a standard deviation of 0.02 M_{⊙}. This enabled us to calculate more reliable estimates for the uncertainties of the semimajor axes and RV masses.
2.2. Model comparisons
We compare the models with differing number of planetary companions using the Bayesian model probabilities. The probability of the kth model is defined as (2)where the marginal integral is (3)f(r θ_{k},ż_{k}) is the likelihood function, and p(θ_{k} ż_{k}) the prior density of the model parameters. In addition, P(ż_{k}) is the prior probability of the kth model and p is the greatest number of planetary signals in our analyses.
We interpret the probabilities of Eq. (2) as proper probabilities and require that the probability of confidently finding a kth planetary signal requires that P(ż_{k}) ≫ P(ż_{k − 1}). In practice, according to Jeffreys (1961), we require that the probability of finding k signals is at least 150 times greater than that of finding k − 1 signals to be able to claim confidently that there are k planets orbiting the target star. In addition, we require that the probability density of the kth planet has a unique maximum that can be interpreted as a Keplerian signal and is not likely to be caused by gaps in the data or by pure noise.
Since the integral in Eq. (3) cannot be calculated analytically, we calculate its value from the sample from the posterior density using the technique discussed in Chib & Jeliazkov (2001).
3. Results
We modelled the system by making two different assumptions about the orbital eccentricities of the planets. First, we treated the eccentricities as free parameters of the model, letting them vary freely in the samplings of the parameter spaces. Second, following the analysis of Vogt et al. (2010), we tested a case where the eccentricities were only allowed to have low values – values below 0.2 – for the sake of dynamical stability of the system.
3.1. Eccentricities as free parameters
When adopting the orbital eccentricities of the k planets in the system as free parameters of the model, we receive the model probabilities (P_{A}) in Table 1. The probabilities of both model sets P_{A} and P_{B}, with different assumptions on the eccentricity, are scaled to unity according to Eq. (2).
Bayesian model probabilities of k planet models with free eccentricity (P_{A}) and eccentricity limited to values below 0.2 (P_{B}).
The fourplanet solution of GJ 581 RV’s. Maximum a posteriori estimates of the parameters and their sets.
In Table 1, we present the probabilities of the different models with k = 0,...,6 planetary companions. We divide the model with k = 5 into two different solutions. In these solutions, 5a and 5b are used to denote the 5companion models including the four planets reported by Mayor et al. (2009) and a fifth Keplerian signal corresponding to one or other of the two signals proposed by Vogt et al. (2010) with periods of roughly 37 and 433 days, respectively. The 6 planet model includes all six planets reported by Vogt et al. (2010). Clearly, according to the probabilities, it cannot be concluded that there are the signals of six planetary companions in the data set. Instead, allowing the values of the eccentricities to be determined freely by the data leads to a conclusion that there are clear signals of at least four companions in the data but the 4companion model has such high probability that it cannot be ruled out confidently enough to claim that there are more than four Keplerian signals in the data. Hence, taking into account the uncertainties of the orbital and other parameters of the models leads to results that contradict with those of Vogt et al. (2010).
The 5planet solution with eccentricities as free parameters consists broadly of two clear probability maxima for the 5th planet (the roughly 37 and 433 day periodicities) but we cannot conclude that the corresponding signals are real as opposed to artefacts of the data. Interestingly, the 433 day periodicity proposed by Vogt et al. (2010) has a greater probability than the 37 day periodicity but the former appears to consist of two closely spaced maxima in the periodicity space. There is one additional maxima in the vicinity of this period (at roughly 465 days), as can be seen in Fig. 2.
Fig. 1 The distribution of the orbital eccentricity of GJ 581 d from the 4companion solution. A Gaussian curve with the same mean and variance is shown for comparison together with the mode, mean (μ), standard deviation (σ), skewness (μ^{3}), and kurtosis (μ^{4}) of the distribution. 

Open with DEXTER 
If the orbital stability of the system is not considered, there is evidence in the data in favour of only four planetary companions. Curiously, we receive an interesting probability distribution for the eccentricity of planet GJ 581 d. This distribution is shown in Fig. 1 and it supports the conclusion of Mayor et al. (2009), who claimed that this companion has an eccentricity of 0.38 ± 0.09. However, there is another maximum in this distribution close to zero, which makes the eccentricity of GJ 581 d consistent with zero. The eccentricities of the other companions were also consistent with zero, yet that of GJ 581 c peaked at 0.1 – a value consistent with the results of Mayor et al. (2009).
According to our results, the point estimates of the orbital parameters and especially their uncertainty estimates need to be revised. The revised parameters and their 99% Bayesian credibility sets , as defined in e.g. Tuomi & Kotiranta (2009), are shown in Table 2.
Fig. 2 Same as Fig. 1 but for the period of GJ 581 f in the 5companion solution. 

Open with DEXTER 
3.2. Low eccentricities
We further assumed that the orbital eccentricities had values lower than 0.2 and repeated the analyses in the previous subsection. The results of these analyses show that the five and sixplanet models do not have great enough posterior probabilities to be able to claim that there are more than four planetary companions orbiting GJ 581 (Table 1). These results show that the sixcompanion solution proposed by Vogt et al. (2010) cannot be considered to imply the existence of six planets in the system because the fourcompanion model cannot be shown to be an insufficient description of the data confidently enough. Instead, the solution of Mayor et al. (2009) remains the most convincing solution even with the combined HARPS and HIRES dataset. Also, the updated parameter and uncertainty estimates of this solution are those presented in Table 2.
Regarding the existence of the proposed companion in the habitable zone of GJ 581 with an orbital period of roughly 37 days, the posterior probabilities in Table 1 imply that there is no evidence in favour of the existence of this companion. It is more probable that there is a companion corresponding to the periodicity of 433 days but even this periodicity appears to not be probable enough and also rather poorly constrained, as shown in Fig. 2.
We also note, that limiting the eccentricities of the companions to values lower than 0.2 favours the sixcompanion interpretation presented in Vogt et al. (2010). The reason is that the eccentricity of GJ 581 d likely differs from zero and the resulting solution where eccentricities are not limited describes the data much better than the limited case.
Based on the data alone, the freely varying eccentricity is a much more likely scenario than the limited one with roughly 25 times greater probability for the fourcompanion model and more than 100 times greater probability for the fivecompanion model. However, they are almost equal for the sixcompanion model, which means that the hypothetical sixcompanion model favours negligible eccentricities.
3.3. Stellar jitter
For the 4companion solution, our estimate of stellar jitter of 1.89 [1.51, 2.36] m s^{1} appears to differ significantly from the value reported by Mayor et al. (2009) of 1.2 m s^{1}. First, we note that the jitter does not only correspond to the noise caused by the stellar surface, but it contains all the excess variations in the data not explained by the Keplerian model or the instrument noise. With measurements from two telescopeinstrument combinations, any small differences in the calibration of the telescopes and instruments may cause systematic differences to the measurements and lead to increased values for the jitter. Also, undetected planets may increase the jitter between our results and those of Mayor et al. (2009) because the observational timeline is longer in the combined dataset analysed here.
Second, our lower limit of the jitter, or more accurately excess noise, in the combined dataset is 1.51 m s^{1} – a value reasonably close to the estimate of Mayor et al. (2009). The fact that we used posterior sampling technique, and took into account the uncertainties of all the parameters, can result in a greater estimate for the jitter magnitude because the common χ^{2} fitting, where jitter is not a free parameter, yields the lowest possible value for the jitter – fitting only the orbital parameters and reference velocities essentially corresponds to minimising the excess noise, not finding the most probable solution.
3.4. Orbital stability
We studied the orbital stability of out fourplanet solution (Table 2) briefly by using numerical integrations of the planetary orbits. We used the BulirschStoer algorithm (Bulirsch & Stoer 1966) because it is a relatively fast and reliable algorithm for studying the dynamics of planetary systems. To assess the instability of the system, we used the concept of Lagrange stability, in which the orbital parameters remain inside some bounded subspace of the parameter space for stable systems. Therefore, we considered the system not stable if there was orbital crossing, collision between the planets, accretion of a planet by the star, or if a planet escaped from the system with a exceeding 100 AU.
We drew a sample of 50 parameter vectors from the parameter posterior density by weighting the sample towards the higheccentricity orbits – the orbital configurations most likely to be unstable – allowed by our solution. We integrated the orbits for 10 000 years for each parameter vector because during that period of time even the outer planet would have completed more than 50 000 orbits. We checked whether the semimajor axes or orbital ecentricities evolved significantly during these integrations.
In our the numerical integrations of the planetary orbits, the semimajor axes and eccentricities remained bounded in all but ten integrations. The semimajor axes remained almost constant and the orbital ecentricities of the three inner planets librated roughly between 0 and 0.2, whereas the eccentricity of the outer planet librated only slightly around its initial value. Also, despite the possible moderate eccentricity of the outer planet of even more than 0.5, it did not appear to have a significant effect on the orbits of the three inner companions that would have resulted in orbital crossings. In ten cases, the initial eccentricities of planet c and e were greater than 0.2, where their probability densities have already become very low with respect to the MAP values. These eccentricities resulted in orbital crossings and it can be therefore concluded that these planets are likely to have eccentricities lower than 0.2. This is a consequence of the tight packing of the three innermost planets in the system.
With these constraints, our sample from the posterior appeared to correspond to stable orbital configurations. Therefore, though having analysed the stability only briefly, we conclude that our revised solution likely corresponds to a stable system.
4. Conclusions and discussion
According to our analyses of the combined RV data for Gliese 581 from HARPS and HIRES spectrographs, it cannot be concluded that there are six planetary companions orbiting Gliese 581. We find that there are confidently four keplerian signals in the combined data set and that there is no evidence in favour of the existence of a lowmass planet in the habitable zone on Gliese 581 with an orbital period of 37 days. Therefore, we conclude that our fourcompanion solution is consistent with the results of Mayor et al. (2009), though the uncertainties of the orbital parameters and RV masses require revision. We also conclude that the interpretation of Vogt et al. (2010), who claimed that there are as many as six planets orbiting Gliese 581, appears to have been not supported by the data strongly enough. There is a periodicity of roughly 433 days in the combined data but it is not probable enough to claim that it is a real signal as opposed to an artefact of the data spacing or pure noise.
We do realise that while the combined dataset provides evidence in favour of the existence of only four companions orbiting GJ 581, it cannot be claimed that the proposed solution of Vogt et al. (2010) is not real – it is fairly possible that the proposed companions f and g do exist in the system. However, the current data do not imply their existence in a statitically significant way, and the solution of the combined dataset remains that in Table 2.
The most likely reasons for the fact that we could not verify the results of Vogt et al. (2010) are that we i) treated the orbital eccentricities of the planets as free parameters throughout the analyses; ii) adopted the magnitude of the stellar RV jitter as a free parameter of our statistical model instead of fixing its value to some a priori determined value or minimising it; and iii) calculated the posterior probabilities of the different models that take into account the Occamian principle of parsimony in a consistent manner – i.e. penalise the models more strongly the more free parameters they contain. We note that our estimate for the stellar jitter of 1.89 m s^{1} contradicts the value obtained by Vogt et al. (2010) of 1.4 m s^{1} but is very close to that estimated from the data of Wright (2005) of 1.9 m s^{1}.
According to the dynamical studies of our solution (Table 2), our revised solution is likely stable – a result consistent with the analyses of Mayor et al. (2009). However, the availability of a sample from the posterior density of the orbital parameters allows a more detailed dynamical study of the Gliese 581 system. Such a study could help ruling out some subsets of the parameter space that appear to have a high probability based on the data alone and could not be shown to correspond to unstable configurations in this study because of the relatively low integration times.
Additional highprecision measurements of the lowmass target Gliese 581 are needed to assess the number of superEarths in its orbit. Also, the first lowmass planet confidently orbiting its hoststar within the limits of the stellar habitable zone remains to be confirmed with confidence.
Acknowledgments
The author was supported by RoPACS (ROcky Planets Around Cool Stars), a Marie Curie Initial Training Network funded by the European Commission’s Seventh Framework Programme.
References
 Bonfils, X., Forveille, T., Delfosse, X., et al. 2005, A&A, 443, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bulirsch, R., & Stoer, J. 1966, Numer. Math., 8, 1 [CrossRef] [Google Scholar]
 Chib, S., & Jeliazkov, I. 2001, J. Am. Stat. Ass., 96, 270 [CrossRef] [Google Scholar]
 Delfosse, X., Forveille, T., Segransan, D., et al. 2000, A&A, 364, 217 [NASA ADS] [Google Scholar]
 Ford, E. B. 2006, ApJ, 642, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Gregory, P. C. 2007a, MNRAS, 374, 1321 [NASA ADS] [CrossRef] [Google Scholar]
 Gregory, P. C. 2007b, MNRAS, 381, 1607 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Haario, H., Saksman, E., & Tamminen, J. 2001, Bernoulli, 7, 223 [CrossRef] [MathSciNet] [Google Scholar]
 Hastings, W. 1970, Biometrika, 57, 97 [CrossRef] [MathSciNet] [Google Scholar]
 Jeffreys, H. 1961, The Theory of Probability (The Oxford Univ. Press) [Google Scholar]
 Mayor, M., Bonfils, X., Forveille, T., et al. 2009, A&A, 507, 487 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., et al. 1953, J. Chem. Phys., 21, 1087 [NASA ADS] [CrossRef] [Google Scholar]
 Tuomi, M., & Kotiranta, S. 2009, A&A, 496, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Udry, S., Bonfils, X., Delfosse, X., et al. 2007, A&A, 469, L43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vogt, S., Butler, P., Rivera, E., et al. 2010, ApJ, 723, 954 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, J. T. 2005, PASP, 117, 657 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Bayesian model probabilities of k planet models with free eccentricity (P_{A}) and eccentricity limited to values below 0.2 (P_{B}).
The fourplanet solution of GJ 581 RV’s. Maximum a posteriori estimates of the parameters and their sets.
All Figures
Fig. 1 The distribution of the orbital eccentricity of GJ 581 d from the 4companion solution. A Gaussian curve with the same mean and variance is shown for comparison together with the mode, mean (μ), standard deviation (σ), skewness (μ^{3}), and kurtosis (μ^{4}) of the distribution. 

Open with DEXTER  
In the text 
Fig. 2 Same as Fig. 1 but for the period of GJ 581 f in the 5companion solution. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.