Free Access
Volume 543, July 2012
Article Number A52
Number of page(s) 12
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201118518
Published online 25 June 2012

© ESO, 2012

1. Introduction

Over the recent years, radial velocity surveys of nearby stars have provided detections of several exoplanet systems with multiple low-mass planets, even few Earth-masses, in their orbits (e.g. Lovis et al. 2011; Mayor et al. 2009, 2011). These systems include a four-planet system around the M-dwarf GJ 581 (Mayor et al. 2009), which has been proposed to possibly have a habitable planet in its orbit (von Paris et al. 2011; Wordsworth et al. 2010), a system of likely as many as seven planets orbiting HD 10180 (Lovis et al. 2011), and several systems with 3 − 4 low-mass planets, e.g. HD 20792 with minimum planetary masses of 2.7, 2.4, and 4.8 M (Pepe et al. 2011) and HD 69830 with three Neptune-mass planets in its orbit (Lovis et al. 2006).

Currently, one of the most accurate spectrographs used in these surveys is the High Accuracy Radial Velocity Planet Searcher (HARPS) mounted on the ESO 3.6 m telescope at La Silla, Chile (Mayor et al. 2003). In this article, we re-analyse the HARPS radial velocities of HD 10180 published in Lovis et al. (2011). These measurements were reported to contain six strong signatures of low-mass exoplanets in orbits ranging from 5 days to roughly 2000 days and a possible seventh signal at 1.18 days. These planets include five 12 to 25 M planets classified in the category of Neptune-like planets, a more massive outer planet with a minimum mass of 65 M, and a likely super-Earth with a minimum mass of 1.35 M orbiting the star in close proximity (Lovis et al. 2011). While the confidence in the existence of the six more massive companions in this system is fairly high, it is less so for the innermost super-Earth (Feroz et al. 2011). Yet, even if the radial velocity signal corresponding to this low-mass companion were an artefact caused by noise and data sampling or periodic phenomena of the stellar surface, the HD 10180 system would be second only to the solar system with respect to the number of planets in its orbits, together with the transiting Kepler-11 six-planet system (Lissauer et al. 2011).

In this article, we re-analyse the radial velocity data of HD 10180 using posterior samplings and model probabilities. We perform these analyses to verify the results of Lovis et al. (2011) with another data analysis method, to calculate accurate uncertainty estimates for the planetary parameters, and to see if this data set contains additional statistically significant periodic signals that could be interpreted as being of planetary origin.

2. Observations of the HD 10180 planetary system

The G1 V star HD 10180 is a relatively nearby and bright target for radial velocity surveys with a Hipparcos parallax of 25.39  ±  0.62 mas and V = 7.33 (Lovis et al. 2011). It is a very inactive (log RHK =  −5.00) solar-type star with similar mass and metallicity (m = 1.06 ± 0.05, [Fe/H]  = 0.08 ± 0.01) and does not appear to show any well-defined activity cycles based on the HARPS observations (Lovis et al. 2011). When announcing the discovery of the planetary system around HD 10180, Lovis et al. (2011) estimated the excess variations in the HARPS radial velocities, usually referred to as stellar jitter, to be very low, approximately 1.0 m s-1. These properties make this star a suitable target for radial velocity surveys and enable the detection of very low-mass planets in its orbit.

Lovis et al. (2011) announced in 2010 that HD 10180 is host to six Neptune-mass planets in its orbit with orbital periods of 5.76, 16.36, 49.7, 123, 601, and 2200 days, respectively. In addition, they reported a 1.18 days power in the periodogram of the HARPS radial velocities that, if caused by a planet orbiting the star, would correspond to a minimum mass of only 35% more than that of the Earth. These claims were based on 190 HARPS measurements of the variations in the stellar radial velocity between November 2003 and June 2009.

The HARPS radial velocities have a baseline of more than 2400 days, which enabled Lovis et al. (2011) to constrain the orbital parameters of the outer companion in the system with almost similar orbital period. In addition, these HARPS velocities have an estimated average instrument uncertainty of 0.57 m s-1 and a relatively good phase coverage with only seven gaps of more than 100 days, corresponding to the annual visibility cycle of the star in Chile.

3. Statistical analyses

We analysed the HD 10180 radial velocities using a simple model that contains k Keplerian signals that are assumed to be caused by non-interacting planets orbiting the star. We also assumed that any post-Newtonian effects are negligible in the timescale of the observations. Our statistical models are then those described in e.g. Tuomi & Kotiranta (2009) and Tuomi (2011), where each radial velocity measurement was assumed to be caused by the Keplerian signals, some unknown reference velocity about the data mean, and two Gaussian random variables with zero means representing the instrument noise with a known variance as reported for the HD 10180 data by Lovis et al. (2011) together with the radial velocities, and an additional random variable with unknown variance that we treat as a free parameter of our model. This additional random variable describes the unknown excess noise in the data caused by the instrument and the telescope, atmospheric effects, and the stellar surface phenomena.

Clearly, the assumption that measurement noise has a Gaussian distribution might be limiting if it were actually more centrally concentrated, had longer tails, was skewed, or was dependent on time and other possible variables, such as stellar activity levels. However, with “only” 190 radial velocities it is unlikely that we could spot non-Gaussian features in the data reliably. For this reason, and because as far as we know the Gaussian one is the only noise model used when analysing radial velocity data, we restrict our statistical models to Gaussian ones.

3.1. Posterior samplings

We analysed the radial velocities of HD 10180 using the adaptive Metropolis posterior sampling algorithm of Haario et al. (2001) because it appears to be efficient in drawing a statistically representative sample from the parameter posterior density in practice (Tuomi 2011; Tuomi et al. 2011). This algorithm is essentialy an adaptive version of the famous Metropolis-Hastings algorithm (Metropolis et al. 1953; Hastings 1970), which adapts the proposal density to the information gathered up to the ith member of the chain when proposing the i + 1th member. It uses a Gaussian multivariate proposal density for the parameter vector θ and updates its covariance matrix Ci + 1 using (1)where is the mean of the parameter vector, (·)T is used to denote the transpose, I is the identity matrix of suitable dimension, ϵ is some very small number that enables the correct ergodicity properties of the resulting chain (Haario et al. 2001), and parameter s is a scaling parameter that can be chosen as 2.42K-1, where K is the dimension of θ, to optimise the mixing properties of the chain (Gelman et al. 1996).

We calculated the marginal integrals needed in model selection using the samples from posterior probability densities with the one-block Metropolis-Hastings (OBMH) method of Chib & Jeliazkov (2001), also discussed in Clyde et al. (2007). However, since the adaptive Metropolis algorithm is not exactly a Markovian process, only asymptotically so (Haario et al. 2001), the method of Chib & Jeliazkov (2001) does not necessarily yield reliable results. Therefore, after a suitable burn-in period used to find the global maximum of the posterior density, during which the proposal density also converges to a multivariate Gaussian that approximates the posterior, we fixed the covariance matrix to its present value, and continue the sampling with the Metropolis-Hastings algorithm, which enables the applicability of the OBMH method.

3.2. Prior probability densities

The prior probability densities of Keplerian models describing radial velocity data have received little attention in the literature. Ford & Gregory (2007) proposed choosing the Jeffreys’ prior for the period (P) of the planetary orbit, the radial velocity amplitude (K), and the amplitude of “jitter”, i.e. the excess noise in the measurements (σJ). This choice was justified because Jeffreys’ prior makes the logarithms of these parameters evenly distributed (Ford & Gregory 2007; Gregory 2007a). We used this functional form of prior densities for the orbital period and choose the cutoff periods such that they correspond to the one-day period, below which we do not expect to find any signals in this work, and a period of 10Tobs, where Tobs is the time span of the observations. We chose this upper limit because it enables the detection of linear trends in the data corresponding to long-period companions whose orbital period cannot be constrained, but is not much greater than necessary in practice (see e.g. Tuomi et al. 2009, for the detectability limits as a function of Tobs), which could slow down the posterior samplings by increasing the hypervolume of the parameter space with reasonably high likelihood values. Anyhow, if the period of the outermost companion cannot be constrained from above, it would violate our detection criterion of the previous subsection.

Unlike in Ford & Gregory (2007), we did not use Jeffreys’ prior for the radial velocity amplitudes nor the excess noise parameter. Instead, because the HARPS data of HD 10180 deviate about their mean less than 20 m s-1, we used uniform priors for these parameters as π(Ki) = π(σJ) = U(0,aRV), for all i, and used a similar uniform prior for the reference velocity (γ) about the mean as π(γ) = U( − aRV,aRV), where we chose the parameter of these priors as aRV = 20 m s-1. While the radial velocity amplitudes could in principle have values greater than 20  m s-1 while the data would still not deviate more than that about the mean, we do not consider that possibility a feasible one.

Following Ford & Gregory (2007), we chose uniform priors for the two angular parameters in the Keplerian model, the longitude of pericentre (ω) and the mean anomaly (M0). However, we did not set the prior of orbital eccentricity (e) equal to a uniform one between 0 and 1. Instead, we expect high eccentricities to be less likely in this case because there are at least six, likely as many as seven, known planets orbiting HD 10180. Therefore, high eccentricities would result in instability and therefore we did not consider their prior probabilities to be equal to the low eccentricity orbits. Our choice was then a Gaussian prior for the eccentricity, defined as (with the corresponding scaling in the unit interval), where the parameter of this prior model is set as σe = 0.3. This choice penalises the high-eccentricity orbits in practice but still allows them if the data claim it. In practice, with respect to the weight this prior puts on zero eccentricity, it gives the eccentricities of 0.2, 0.4, and 0.6 relative weights of 0.80, 0.41, and 0.14, indicating that this prior can only have a relatively minor effect on the posterior densities.

Finally, we required that the planetary systems corresponding to out Keplerian solution to the data did not have orbital crossings between any of the companions. We used this condition as additional prior information by estimating that the likelihood of having any two planets in orbits that cross one another is zero. We could have used a more restrictive criteria, such as the requirement that the planets do not enter each other’s Hill spheres at any given time, but decided to keep the situation simpler because we wanted to see whether the orbital periods of the proposed companions are constrained by data instead of stability criteria, as described in the next subsection.

This choice of restricting the solutions in such a way that the corresponding planetary system does not suffer from destabilising orbital crossings also helps reducing the computational requirements by making the posterior samplings simpler. After finding k Keplerian signals in the data, we simply searched for additional signals between them by limiting the period space of the additional signals between these k periods. We set the initial periods of the k planets close to the solution of the k-Keplerian model and performed k + 1 samplings where each begins with the period of the k + 1th signal in different “gaps” around the previously found k signals, i.e. corresponding to a planet inside the orbit of the innermost one, between the two innermost ones, and so forth. If a significant k + 1th signal is not found in one of these “gaps”, the corresponding solution can simply be neglected. However, if there are signals in two or more gaps, it is straightforward to determine the most significant one because they can be treated as different models containing the same exact number of parameters. We then chose the most significant periodicity as the k + 1th one and continued testing whether there are additional signals in the data. In this way, the problem of being able to rearrange the signals in any order, that would cause the posterior density to be actually highly multimodal (Feroz et al. 2011), actually disappears because in a given solution the orbital crossings are forbidden and the ordering of the companions remains fixed.

3.3. Detection threshold

While the Bayesian model probabilities can be used reliably in assessing the relative posterior probabilities of models with differing numbers of Keplerian signals (e.g. Ford & Gregory 2007; Gregory 2011; Loredo et al. 2011; Tuomi & Kotiranta 2009; Tuomi 2011), we introduced additional criteria to make sure that the signals we detect can be interpreted as being of planetary origin and not arising from unmodelled features in the measurement noise or as spurious signals caused by measurement sampling. Our basic criterion is that the posterior probability of a model with k + 1 Keplerian signals has to exceed 150 times that of a model with only k Keplerian signals to claim that there are k + 1 planets orbiting the target star. We chose this threshold probability based on the considerations of Kass & Raftery (1995).

We require that the signals we detect in the measurements have radial velocity amplitudes, Ki for all i, statistically distinguishable from zero. In practice this means that not only the maximum a posteriori estimate is clearly greater than zero, but that the corresponding Bayesian δ credibility sets, as defined in e.g. Tuomi & Kotiranta (2009), do not allow the amplitude to be negligible with δ = 0.99, i.e. with a probability of 99%. The second criterion is that the periods of all signals (Pi) are well-defined by the posterior samples in the sense that they can be constrained from above and below and are not constrained purely by the condition that orbital crossings corresponding to the planetary orbits are not allowed.

To further increase the confidence of our solutions, we did not set the prior probabilities of the different models equal in our analyses. We suspect a priori that detecting k + 1 planets would be less likely than detecting k planets in any given system. In other words, we estimated that any set of radial velocity data would be more likely to contain k Keplerian signals than k + 1. Therefore, we set the a priori probabilities of models with k and k + 1 Keplerian signals such that P(k) = 2P(k + 1) for all k, i.e. we penalise the model with one additional planet by a factor of two. Because of this subjective choice, if the model with k + 1 Keplerian signals receives a posterior probability that exceeds our detection threshold of being 150 times greater than that of the corresponding model with k Keplerian signals, we are likely underestimating the confidence level of the model with k + 1 Keplerian signals relative to a uniform prior.

Physically, this choice of prior probabilities for different models corresponds to the fact that the more planets there are orbiting a star, the less stable orbits there will be left. Therefore, we estimate that if k planets are being detected, there is naturally “less room” for an additional k + 1th companion. However, this might be true statistically, not in an individual case, which leaves room for discussion. Yet, this and the benefit that we underestimate the significance of any detected signals encourages us to use this prior.

3.4. Frequentist and Bayesian detection thresholds

In addition to the Bayesian analyses described in the previous subsections, we analysed the residuals of each model using the Lomb-Scargle periodograms (Lomb 1976; Scargle 1982). As in Lovis et al. (2011), we plotted the 10%, 1%, and 0.1% false-alarm probabilities (FAPs) to the periodograms to see the significance of the powers they contain. However, because Lovis et al. (2011) calculated the FAPs in a frequentist manner by performing random permutations to the residuals while keeping the exact epochs of the data fixed and by seeing how often this random permutation produces the observed powers, we first put this periodogram approach into its philosophical context.

Generating N random permutations of the residuals for each model aims at simulating a situation where it would be possible to receive N independent sets of measurements from the system of interest and seeing how often the measurement noise generates the signals corresponding to the highest powers in the periodogram. While they would not have been independent even in the case this process tries to simulate, because the exact epochs are fixed making the measurements actually dependent through the dimension of time (the measurement is actually a vector of two numbers, radial velocity, and time), this approach suffers from another more significant flaw. The uncertainties of the signals removed from the data cannot be taken into account, which means that the method assumes the removed signals were known correctly. Obviously this is not the case even with the strongest signals, and even less so for the weaker ones, making the process prone to biases. Therefore, while likely producing reliable results when the signals are clear and their periods can be constrained accurately, this method cannot be expected to provide reliable results for extremely weak signals with large (and unknown) uncertainties. As noted by Lovis et al. (2011), when testing the significance of the 600-day signal, they could not take into account the uncertainties of the parameters of each Keplerian signal, that of the reference velocity, or the uncertainty in estimating the excess noise in the data correctly.

The above “frequentist” way of performing the analyses and intepreting the consequent results is different from the Bayesian one. Because we only received one set of data, we have to base all our results on that and not some hypothetical data that would have corresponded to a repetition of the original measurements. The Bayesian philosophy is to infer all the information from the data by combining them with our prior beliefs on what might be producing them. For instance, as described above when discussing our choice of priors, we can expect tightly packed multiplanet systems to be more likely to contain planets on close-circular orbits than on very eccentric ones. Also, with the powerful posterior sampling algorithms available, it is possible to take the uncertainties in every parameter into account simultaneously, which enables the detection of weak signals in the data (e.g. Gregory 2005, 2007a,b; Tuomi & Kotiranta 2009) and prevents the detection of false positives, as happened in the case of Gliese 581 (Vogt et al. 2010; Gregory 2011; Tuomi 2011).

Yet, despite the above problems in the traditional periodogram analyses, we take advantage of the power spectra of the residuals in our posterior samplings. The highest peaks in the periodogram can be used very efficiently together with Bayesian methods by using the corresponding periodicities as initial states of the Markov chains in the adaptive Metropolis algorithm. Because of this choice, the initial parameter vector of the Markov chain starts very close to the likely maximum a posteriori (MAP) solution, which makes its convergence to the posterior density reasonably rapid and helps reducing computational requirements.

4. Results

When drawing a sample from the parameter posterior density and using it to calculate the corresponding model probabilities, it became crucial that this sample was a statistically representative one. While posterior samplings generally provide a global solution, it is always possible that the chain converges to a local maximum and stays in its vicinity within a sample of finite size. To make sure that we indeed received the global solution, we calculated several Markov chains starting from the vicinity of the apparent MAP solution and compared them to see that they indeed corresponded to the same posterior probability density. In practice, sampling the parameter spaces was computationally demanding because the probability that the parameter vectors drawn from the Gaussian proposal density are close to the posterior maximum decreases rapidly when the number of parameters with non-Gaussian probability densities increases. Therefore, while models with 0 − 6 planets were reasonably easy to sample and we received acceptance rates of 0.1 − 0.3, these rates decreased when increasing the number of signals in the model further. As a result, for models with 8 − 9 Keplerian signals, the acceptance rates decreased below 0.1 and forced us to increase the chain lenghts by two orders of magnitude from a typical 106 to as high as 108.

In the following subsections, the results are based on several samplings that yielded the same posterior densities, and also consistent model probabilities.

4.1. The number of significant periodicities

The posterior probabilities of the different models provide information on the number of Keplerian signals (k) in the data. Lovis et al. (2011) were confident that there are six planets orbiting HD 10180 based on their periodogram-based analyses of model residuals and the corresponding random permutations of them when calculating the significance levels of the periodogram powers. They also concluded that the six planets in the system with increasing periods of 5.75962  ±  0.00028, 16.3567  ±  0.0043, 49.747  ±  0.024, 122.72  ±  0.20, 602  ±  11, and 2248 days comprise a stable system given that the masses of the planets are within a factor of  ~ 3 from the minimum masses of 13.70  ±  0.63, 11.94  ±  0.75, 25.4  ±  1.4, 23.6  ±  1.7, 21.4  ±  3.0, and 65.3  ±  4.6 M, respectively.

Because Lovis et al. (2011) pointed out that there are actually two peaks in the periodogram corresponding to the seventh signal, i.e. 1.18 and 6.51 day periodicities that are the one-day aliases of each other, but noted that the 6.51 signal, if corresponding to a planet, would cause the system to be unstable on short timescales, we adopted the 1.18 periodicity as the seventh signal in the data. The relative probabilities of the models with k = 0,...,9 are shown in Table 1 together with the period (Ps) of the next Keplerian signal added to the model.

Table 1

The relative posterior probabilities of models with k = 0,...,9 Keplerian signals (ℳk) given radial velocities of HD 10180 (or data, d) together with the periods (Ps) of the signals added to the model when increasing the number of signals in the model by one.

According to the model probabilities in Table 1, the eight- and nine-Keplerian models are the most probable descriptions of the processes producing the data out of those considered. Improving the statistical model by adding the seventh signal, with a period of 1.18 days, increases the model probability by a factor of more than 106, which makes the credibility of this seventh signal high. Adding two more signals corresponding to 67.6 and 9.66 day periodicities increases the model probabilities even more. As a result, the nine-Keplerian model receives the greatest posterior probability of slightly more than 150 times more than the next best model, the eight-Keplerian model. This enables us to conclude that there is strong evidence in favour of the hypothesis that the 67.6 and 9.66 days periodicities are not produced by random processes, i.e. measurement noise.

4.2. Periodograms of residuals

We subtracted the models with six to nine periodic signals from the data and calculated the Lomb-Scargle periodograms (Lomb 1976; Scargle 1982) of these residuals (Fig. 1) together with the standard analytica FAPs. As already seen in Lovis et al. (2011), the two strong powers corresponding to a 1.18 day periodicity and its 1-day alias at 6.51 days are strong in the residuals of the six-Keplerian model (top panel in Fig. 1) and exceed the 1% FAP. These peaks are also removed from the residuals of the seven-Keplerian model (second panel). However, it can be seen that the 9.66 and 67.6 day periods have strong powers in these residuals, yet neither of them exceeds even the 10% FAP level. Modelling these periodicities as Keplerian signals and plotting the periodograms of the corresponding residuals of the nine-Keplerian model (bottom panel) shows that there are no strong powers left in the residuals. While this does not indicate that these two periodicities are significant, it shows that they are clearly present in the periodograms of the residuals and support the findings in the previous subsection based on the model probabilities.

thumbnail Fig. 1

Lomb-Scargle periodograms of the HD 10180 radial velocities for the residuals of the models with six (top) to nine (bottom) periodic signals. The dotted, dashed, and dot-dashed lines indicate the 10%, 1%, and 0.1% FAPs, respectively.

We note that there appear to be two almost equally strong peaks in the eight-Keplerian model residuals (Fig. 1, panel 3). However, these powers corresponding to periods of 9.66 and 1.11 days are one-day aliases of one another. This aliasing is clear because the 1.11 day power is absent in the periodogram of the nine-Keplerian model (Fig. 1, bottom panel).

4.3. Planetary interpretation: orbital parameters

Because of the samplings of parameter posterior densities of each statistical model, we were able to calculate the estimated shapes of the parameter distributions for each model and use these to estimate the features in the corresponding densities. We describe these densities using three numbers, the MAP estimates and the corresponding 99% Bayesian credibility sets (BCSs), as defined in e.g. Tuomi & Kotiranta (2009). The simple MAP point estimates and the corresponding 99% BCSs do not represent these dentities very well because some of the parameters are highly skewed and have tails on one or both sides. However, in this subsection we use these estimates when listing the parameters and interpret that all the signals we observe are of planetary origin.

When calculating the semimajor axes and minimum masses of the planets, we took the uncertainty in the stellar mass into account by treating it as an independent random variable. We assumed that this random variable had a Gaussian density with mean equal to the estimate given by Lovis et al. (2011) of 1.06 M and a standard deviation of 5% of this estimate. As a consequence, the densities of these parameters are broader than they would be if using a fixed value for the stellar mass, which indicates greater uncertainty in their values.

thumbnail Fig. 2

Phase-folded Keplerian signals of the nine-planet solution with the other eight signals removed.

4.3.1. The six planet solution

The parameter estimates of our 6-Keplerian model are listed in Table 21. This solution is consistent with the solution reported by Lovis et al. (2011) but the uncertainties are slightly greater, likely because we took into account the uncertainty in the jitter parameter σj and because Lovis et al. (2011) used more conservative uncertainty estimates from the covariance matrix of the parameters that does not take the nonlinear correlations between the parameters into account. The greatest difference is therefore in the uncertainty of the orbital period of the outermost companion, whose probability density has a long tail towards longer periods and periods as high as 2670 days cannot be ruled out with 99% confidence (the supremum of the 99% BCS).

Table 2

The six-planet solution of HD 10180 radial velocities.

4.3.2. The nine planet solution

Assuming that all periodic signals in the data are indeed caused by planetary companions orbiting the star, the parameters of out nine-Keplerian solution are listed in Table 3 and the phase-folded orbits of the nine Keplerian signals are plotted in Fig. 2.

Table 3

The nine-planet solution of HD 10180 radial velocities.

Because the simple estimates in Table 3 can be very misleading in practice, especially if there are nonlinear correlations between the parameters, we also plotted the projected distributions of some of the parameters in the appendix. The distributions of Pi, Ki, and ei, for each planet i = 1,...,9 (with indice 1 (9) referring to the shortest (longest) period) are shown in Fig. A.1 and show that the periods of all Keplerian signals are indeed well constrained, the radial velocity amplitudes differ significantly from zero, and the orbital eccentricities peak at or close to zero, indicating likely circular orbits. We also plotted the parameter describing the magnitude of the excess noise, or jitter, in Fig. A.1 to demonstrate that the MAP estimate of this parameter of 1.15 m s-1, with a BCS of [0.92, 1.42] m s-1, is consistent with the estimate of Lovis et al. (2011) of 1.0 m s-1, while this is not the case for the six-Keplerian model for which the σj receives a MAP estimate of 1.40 m s-1 with a BCS of [1.18, 1.70] m s-1 (Table 2).

The periods of all companions are constrained well but that of the 9.66 day signal remains bimodal with two peaks at 9.66 and 9.59 days (Fig. A.1). The 9.59 day peak was found to have lower probability based on our posterior samplings because the Markov chain had good mixing properties in the sense that it visited both maxima several times during all the samplings, the posterior density had always a global maximum at 9.66 days.

4.4. Dynamically allowed orbits

We performed tests of dynamical stability within the context of Lagrange stability to see whether the two additional signals in the HD 10180 radial velocities could correspond to low-mass planets orbiting the star. We used the analytically derived approximated Lagrange stability criterion of Barnes & Greenberg (2006) to test the stability of each subsequent pair of planets in the system. While this analytical criterion is only a rough approximation and only applicable for two-planet systems, it can nevertheless provide useful information on the likely stability or instability of the system. According to Barnes & Greenberg (2006), the orbits of two planets (denoted using subindices 1 and 2, respectively) satisfy approximately the Lagrange stability criterion if (2)where μi = miM-1, α = μ1 + μ2, , , M = m + m1 + m2, ei is the eccentricity, ai is the semimajor axis, mi is the planetary mass, and m is stellar mass.

Using the above relation, we calculated the threshold curves for the ith planet with both the next planet inside its orbit (i − 1th) and the next planet outside its orbit (i + 1th). We used the MAP parameter estimates for the ith planet and calculated the allowed eccentricities of the i − 1th and i + 1th planets as a function of their semimajor axes by using the MAP estimates for their masses.

Because Lovis et al. (2011) found their six-planet solution stable, we used it as a test case when calculating the Lagrange stability threshold curves. We used the six-planet solution in Table 2, and plotted the threshold curves together with the orbital parameters in Fig. 3. In this figure, the shaded areas indicate the likely unstable parameter space and the red circles indicate the positions of the modelled planets in the system.

thumbnail Fig. 3

Approximated Lagrange stability thresholds between each of the two planets and the MAP orbital parameters of the six-planet solution (Table 2).

It can be seen in Fig. 3 that the ith planet has orbital parameters that keeps it inside the Lagrange stability region of the neighbouring planetary companions for all i = 1,...,6. This result then agrees with the numerical integrations of Lovis et al. (2011) and, while only a rough approximation of the reality, encourages us to use the criterion in Eq. (2) for our seven-, eight-, and nine-companion solutions as well. Figure 3 also suggests that there might be stable regions between the orbits of these six planets for additional low-mass companions.

When interpreting all the signals in our nine-planet solution as being of planetary origin, the stability thresholds show some interesting features (Fig. 4). The periodicities at 9.66 and 67.6 days would correspond to planets that satisfy the condition in Eq. (2) if the orbital eccentricities of all the companions were close to or below the MAP estimate, which, according to the probability densities in Fig. A.1, appears to be likely based on the data alone. This means that the planetary origin of these periodicities cannot be ruled out by this analysis.

thumbnail Fig. 4

Approximated Lagrange stability thresholds between each of the two planets and the MAP orbital parameters of the nine-planet solution (Table 3).

In reality, the stability constraints are necessarily more limiting than those described by the simple Eq. (2) because of the gravitational interactions between all the planets, not only the nearby ones. Also, it does not take the stabilising or destabilising effects of mean motion resonances into account. However, the numerical integrations of the orbits of the seven planets performed by Lovis et al. (2011) do not rule out the 0.09 and 0.33 AU orbits (Table 3) out as unstable but show that there are regions of at least “reasonable stability” in the vicinity of these orbits given that they are close-circular and that the planetary masses in these orbits are small. When interpreted as being of planetary origin, the periodic signals at roughly 9.66 and 67.6 days satisfy these requirements.

We note that the selected prior density for the orbital eccentricities, namely for all i, in fact helps slightly in removing a priori unstable solutions from the parameter posterior density. However, this effect is not very significant in this case, because the prior requirement that the signals in the data do not correspond to planets with crossing orbits constrains the eccentricities much more strongly and the parameter posteriors we would receive with uniform eccentricity prior would therefore not differ significantly from those reported in Table 3 and Figs. 2 and A.1.

4.5. Avoiding unconstrained solutions

To further emphasise our confidence in the nine periodic signals in the data, we tried finding additional signals in the gaps between the nine Keplerians, especially, between the 123 and 600 day orbits. This part of the period-space is interesting because any habitable planet in the system would have its orbital period in this space and because the stability thresholds allow the existence of low-mass planets in close-circular orbits in this region (see Fig. 4).

As could already be suspected based on the periodograms of residuals in Fig. 1, we were unable to find any signals between the periods of 123 and 600 days. The sampling of the parameter space of this ten-Keplerian model was much more difficult than that of the models with fewer signals because the orbital period of this hypothetical tenth signal was only constrained by the fact that a priori we did not allow orbital crossings (Fig. 5). As a result, the probability density of the orbital period did not have one clear maximum, but several small maxima, whose relative significance is not known because we cannot be sure whether the Markov chain converged to the posterior in this case (Fig. 5). Therefore, the distribution of the orbital period in Fig. 5 is only a rough estimate of what the density might look like.

thumbnail Fig. 5

Posterior densities of the minimum mass and orbital period of the planet that could exist in the habitable zone of HD 10180 without having been detected using the current data. The solid curve is a Gaussian density with the same mean (μ) and variance (σ2) as the parameter distribution.

Because different samplings yielded similar but not equal densities for the orbital period, we could not be sure whether the chain had indeed converged to the posterior or not. For this reason, we did not consider the corresponding posterior probability of this model trustworthy and do not show it in Table 1. The posterior probabilities we received were roughly 1 − 5% of that of the nine-Keplerian model. However, these samplings still provide some interesting information in the sense that we can put an upper limit to the planetary masses that could exist between the 123 and 600 day periods and still not be detected confidently by the current observations.

According to the samplings of the parameter space, the probability of there being Keplerian signals between the 123 and 600 day periods with radial velocity amplitudes in excess of 1.1 m s-1 is less than 1%. Our MAP solution for this amplitude is 0.12 m s-1 with a BCS of [0.00, 1.10] m s-1. This means that we can rule out the existence of planets more massive than approximately 12.1 M in this period space because such companions would likely have been detected by the current data. Therefore, even the fact that a signal was not detected can help constraining the properties of the system, as seen in the probability density of the minimum mass of this hypothetical planet in Fig. 5 – this signal is clearly indistinguishable from one with negligible amplitude because the density is peaking close to zero. This result means that an Earth-mass planet could exist in the habitable zone of HD 10180 and most likely, if there is a low-mass companion orbiting the star between orbital periods of 123 and 600 days, it has a minimum mass of less than 12.1 M. However, we cannot say much about the possible orbit of this hypothetical companion, because all the orbits between 123 and 600 days without orbital crossings are almost equally probable (Fig. 5). All we can say is that orbital crossings limit the allowed periodicities and yield a 99% BCS of [128, 534] days for this orbital period, though, additional dynamical constraints would narrow this interval even more, as shown in Fig. 4 and Fig. 12 of Lovis et al. (2011).

4.6. Detectability of Keplerian signals

To further emphasise the significance of the nine signals we detect in the HD 10180 radial velocities, we generated an artificial data set to see if known signals could be extracted from it confidently given the definition of our criteria for the detection threshold. This data set was generated by using the same 190 epochs as in the HD 10180 data of Lovis et al. (2011). We generated the radial velocities corresponding to these epochs by using a superposition of the nine signals with parameters roughly as in Table 3. Furthermore, we added three noise components, Gaussian noise with zero mean as described by the uncertainty estimates of each original radial velocity of Lovis et al. (2011), Gaussian noise with zero mean and σ = 1.1 m s-1, and uniform noise as a random number drawn from the interval [− 0.2, 0.2] m s-1. The lattter two produce together the observed excess noise in the data of roughly 1.15 m s-1 when modelled as pure Gaussian noise. We used this different noise model when generating the data to not commit an inverse crime, i.e. to not generate the data using the same model used to analyse it which would correspond to studying the properties of the model only (see e.g. Kaipio & Somersalo 2005).

Analysing these artificial data yielded results confirming the trustworthiness of our methods. Using the detection criteria defined above, we could extract all nine signals from the artificial data with well-constrained amplitudes and periods that were consistent with the added signals in the sense that the 99% BCSs of the parameters contained the values of the added signals. The MAP estimates did differ from the added signals but not significantly so, given the uncertainties as described by distributions corresponding to the parameter posterior density. This simply represents the statistical nature of the solutions based on posterior samplings. As an example, we show the periods and amplitudes of the three weakest signals with the lowest radial velocity amplitudes as probability distributions (Fig. 6). This figure indicates that if the radial velocity noise is indeed dominated by Gaussian noise, the low-amplitude signals we report can be detected confidently. This conclusion was also supported by the corresponding model probabilities we received for the artificial data set. These probabilities indicated that the nine-Keplerian model had the greatest posterior probability exceeding the threshold of 150 times greater than the probability of the eight-Keplerian model.

thumbnail Fig. 6

Distributions of the periods and amplitudes received for the three weakest signals in the artificially generated data. The added signals had Pi = 1.18, 9.66, and 67.55 days, and Ki = 0.78, 0.53, and 0.75 m s-1 for i = 1,3,6, respectively.

4.7. Comparison with earlier results

The analyses of Lovis et al. (2011) of the same radial velocities yielded differing results, i.e. the number of significant powers in periodogram was found to be seven instead of the nine significant periodic signals reported here. While this difference is likely due to the fact that the power spectrums are calculated by fixing the parameters of the previous signals to some point estimates when searching for additional peaks, the analyses of Feroz et al. (2011) suffered from similar sources of bias. The Bayesian approach of Feroz et al. (2011) was basically a search of k + 1th signal in the residuals of the model with k Keplerian signals. These authors derived the probability densities of the residuals by assuming they were uncorrelated, which is unlikely to be the case, especially if there are signals left in the residuals. Therefore, any significance test, in this case the comparison of Bayesian evidence of a null hypothesis and an alternative one with a model containing one more signal, is similarly biased because these correlations are not fully accounted for. While this source of bias might be relatively small, in the approach of Feroz et al. (2011) the effective number of parameters is artificially decreased by the very fact that residuals are being analysed. This decrease, in turn, might make any comparisons of Bayesian evidence estimates biased. Modifying the uncertainties corresponding to the posterior density of the model residuals does not necessarily account for this decrease in the dimension when the weakest signals among the k detected ones are at or below the residual uncertainty and their contribution to the total uncertainty of the residuals becomes negligible in the first place.

Indeed, we were able to replicate the results of Feroz et al. (2011) using the OBMH estimates of marginal integrals. Using their simple method, we received a result that the six-Keplerian model residuals could not be found to contain an additional signal. While the Bayesian evidence for this additional signal did not exceed that of having the six-Keplerian model residuals consist of purely Gaussian noise, the amplitude of this additional signal in the residuals was also found indistinguishable from zero, in accordance with our criteria of not detecting a signal. Therefore, the results we present here do not actually conflict with those of Feroz et al. (2011). As an example, the log-Bayesian evidences of six- and seven-Keplerian models (Table 1) were  − 358.36 and  − 343.73, respectively. Analysing the residuals of the model with k = 6 using 0- and 1-Keplerian models should yield similar numbers, if the method of Feroz et al. (2011) were trustworthy. Instead, we received values of roughly  − 338 for both models. This means that the probability of the null-hypothesis is exaggerated because of the fact that the corresponding model contains only two free parameters, the jitter magnitude and the reference velocity. The Occamian principle cannot therefore penalise this model as much as it should, and the results are biased in favour of the null-hypothesis, which effectively prevents the detection of low-amplitude signals.

We also tested the method of Feroz et al. (2011) in analysing the artificial test data described in the previous subsection. The results were almost similar: the log-Bayesian evidence given the residuals of a model with k Keplerian signals was found to favour the six-Keplerian interpretation, while the artificial data was known to contain nine signals. This further emphasises that the method of Feroz et al. (2011), while capable of detecting the strongest signals in the data, cannot be considered trustworthy if it fails to make a positive detection of a low-amplitude signal. Yet, it is likely trustworthy if it provides a positive detection.

5. Conclusions and discussion

We have re-analysed the 190 HARPS radial velocities of HD 10180 published in Lovis et al. (2011) and reported our findings in this article. First, we revised the orbital parameters of the proposed six planetary companions to this star and calculated realistic uncertainty estimates based on samplings of the parameter space. We also verified the significance of the 1.18 day signal reported by Lovis et al. (2011) and interpreted as arising from a planetary companion with this orbital period and a minimum mass of as low as 1.3 M. In addition to these seven signals, we reported two additional periodic signals that are, according to our model probabilities in Table 1, statistically significant and unlikely to be caused by noise or data sampling or poor phase-coverage of the observations. Their amplitudes are well constrained and differ statistically from zero, which would not be the case unless they corresponded to actual periodicities in the data. We can also constrain their periods from above and below reasonably accurately.

A related analysis of the same radial velocities was recently carried out by Feroz et al. (2011) but they received results differing in the number of significant periodicities. They claimed that only six signals can be reliably detected in the data, as opposed to nine detected in the current work. However, they first analysed the data using a model with k Keplerian signals and analysed the remaining residuals to see if they contained one more, k + 1th signal. While, as in the analyses of Lovis et al. (2011), this approach does not fully account for the uncertainties of the first k signals when they have low amplitudes and their contribution to the residual uncertainty is negligible, it also actually assumes that the k-Keplerian model is a correct one and then tests if it is not so, which is a clear contradiction and, while useful for strong signals, as demonstrated by Feroz et al. (2011), likely prevents the detection of weak signals in the data. This is underlined by the fact that Feroz et al. (2011) assumed the residual vector to have an uncorrelated multivariate Gaussian distribution – this clearly cannot be the case if there are signals left in the residuals. Our analyses are not prone to similar weaknesses.

Because planetary companions orbiting the star would produce the kind of periodicities we observe in the radial velocities, the interpretation of the two new signals as caused by two new low-mass planets seems reasonable. As noted by Lovis et al. (2011), the star is a very quiet one without clear activity-induced periodicities, which makes it unlikely that one or some of the periodic signals in the data were caused by stellar phenomena. Also, the periodicities we report, namely 9.66 and 67.6 days, do not coincide with any periodicities arising from the movement of the bodies in the solar system. Therefore, we consider the interpretation of these two new signals of being as planetary origin to be the most credible explanation. If this were the case, these two signals would correspond to planets on close-circular orbits with minimum masses of 1.9 and 5.1 M, respectively, enabling the classification of them as super-Earths.

Apart from the significance of the signals we observe, there is another reasonably strong argument in favour of the interpretation that all nine signals in the data are indeed of planetary origin. Assuming that they were not, which based on stability reasons is the case with the 6.51 day signal that is quite certainly an alias of the 1.18 day periodicity likely caused by a planet (Lovis et al. 2011), we would expect the weakest signals to be at random periods independent of the six strong periodicities in the data and the seventh 1.18 periodicity. Instead, this is not the case but the two additional signals reported in the current work appear at periods that fall in between the existing ones and, if interpreted as being of planetary origin, likely have orbits that enable long-term stability of the system (Fig. 4) if their orbital eccentricities are close to or below the estimates in Table 3. As stated by Lovis et al. (2011), there are “empty” places in the HD 10180 system that allow dynamical stability of low-mass planets in the orbits corresponding to these empty places in the orbital parameter space, especially in the a − e space. The two periodic signals we observe in the data correspond exactly to those empty places if interpreted as being of planetary origin.

Additional measurements are needed to verify the significance of the two new periodic signals in the radial velocities of HD 10180 and to set tighter constraints on the orbital parameters of the planets in the system. Also, the possibility that all nine signals in the data correspond to planets should be tested by full-scale numerical integrations of their orbits. If all the configurations allowed by our solution in Table 3 were found to correspond to unstable systems in any timescales of less than the estimated stellar age of 4.3  ±  0.5 Gyr (Lovis et al. 2011), it would be a strong argument against the planetary interpretation of one or both of the signals we report in this article or the third low-amplitude 1.18 day signal. However, the results we present here and those in Lovis et al. (2011) suggest that this planetary interpretation of all the signals cannot be ruled out by dynamical analyses of the system.

If the significance of these signals increases when additional high-precision radial velocities become available, and their interpretation as being of planetary origin is confirmed, the planetary system aroung HD 10180 will be the first one to top the solar system in terms of number of planets in its orbits. Furthermore, according to the rough dynamical considerations of the current work and the more extensive numerical integrations of Lovis et al. (2011), there are stable orbits for a low-mass companion in or around the habitable zone of the star. If such a companion exists, its minimum mass is unlikely to exceed 12.1 M according to our posterior samplings of the corresponding parameter space.


Lovis et al. (2011) used letter b for the 1.18 day periodicity not present in the six-Keplerian model. Therefore, we use letters c − h in Table 2 to have the letters denote the same signals. Yet, the solution of Feroz et al. (2011) denotes the 5.76 day signal with letter b.


M. Tuomi is supported by RoPACS (Rocky Planets Around Cools Stars), a Marie Curie Initial Training Network funded by the European Commission’s Seventh Framework Programme. The author would like to thank the two anonymous referees for constructive comments and suggestions that resulted in significant improvements in the article.


  1. Barnes, R., & Greenberg, R. 2006, ApJ, 647, L163 [CrossRef] [Google Scholar]
  2. Chib, S., & Jeliazkov, I. 2001, J. Am. Stat. Ass., 96, 270 [Google Scholar]
  3. Clyde, M. A., Berger, J. O., Bullard, F., et al. 2007, Statistical Challenges in Modern Astronomy IV, ed. G. J. Babu, & E. D. Feigelson, ASP Conf. Ser., 371, 224 [Google Scholar]
  4. Feroz, F., Balan, S. T., & Hobson, M. P. 2011, MNRAS, 415, 3462 [NASA ADS] [CrossRef] [Google Scholar]
  5. Ford, E. B., & Gregory, P. C. 2007, Statistical Challenges in Modern Astronomy IV, ed. G. J. Babu, & E. D. Feigelson, ASP Conf. Ser., 371, 189 [Google Scholar]
  6. Gelman, A. G., Roberts, G. O., & Gilks, W. R. 1996, Efficient Metropolis jumping rules, in Bayesian Statistics V, ed. J. M. Bernardo, J. O. Berger, A. F. David, & A. F. M. Smith, 599 [Google Scholar]
  7. Gregory, P. C. 2005, ApJ, 631, 1198 [NASA ADS] [CrossRef] [Google Scholar]
  8. Gregory, P. C. 2007a, MNRAS, 381, 1607 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  9. Gregory, P. C. 2007b, MNRAS, 374, 1321 [NASA ADS] [CrossRef] [Google Scholar]
  10. Gregory, P. C. 2011, MNRAS, 415, 2523 [NASA ADS] [CrossRef] [Google Scholar]
  11. Haario, H., Saksman, E., & Tamminen, J. 2001, Bernoulli, 7, 223 [CrossRef] [MathSciNet] [Google Scholar]
  12. Hastings, W. 1970, Biometrika 57, 97 [Google Scholar]
  13. Kaipio, J., & Somersalo, E. 2005, Statistical and Computational Inverse Problems, Applied Mathematical Sciences (Springer), 160 [Google Scholar]
  14. Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Ass., 430, 773 [Google Scholar]
  15. Lissauer, J. J., Fabrycky, D. C., Ford, E. B., et al. 2011, Nature, 470, 53 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  16. Lomb, N. R. 1976, Astrophys. Space Sci., 39, 447 [Google Scholar]
  17. Loredo, T. J., Berger, J. O., Chernoff, D. F., et al. 2011, Statistical Methodology, accepted, DOI: 10.1016/j.stamet.2011.07.005 [Google Scholar]
  18. Lovis, C., Mayor, M., Pepe, F., et al. 2006, Nature, 441, 305 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  19. Lovis, C., Ségransan, D., Mayor, M., et al. 2011, A&A, 528, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Mayor, M., Pepe, F., Queloz, D., et al. 2003, Messenger, 114, 20 [Google Scholar]
  21. Mayor, M., Bonfils, X., Forveille, T., et al. 2009, A&A, 507, 487 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Mayor, M., Marmier, M., Lovis, C., et al. 2011, A&A, submitted [arXiv:1109.2497] [Google Scholar]
  23. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., et al. 1953, J. Chem. Phys., 21, 1087 [NASA ADS] [CrossRef] [Google Scholar]
  24. Pepe, F., Lovis, C., Ségransan, D., et al. 2011, A&A, 534, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
  26. Tuomi, M. 2011, A&A, 528, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Tuomi, M., & Kotiranta, S. 2009, A&A, 496, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Tuomi, M., Kotiranta, S., & Kaasalainen, M. 2009, A&A, 494, 769 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Tuomi, M., Pinfield, D., & Jones, H. R. A. 2011, A&A, 532, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. von Paris, P., Gebauer, S., Godolt, M., et al. 2011, A&A, 532, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Vogt, S. S., Butler, R. P., Rivera, E. J., et al. 2010, ApJ, 723, 954 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Wordsworth, R., Forget, F., Selsis, F., et al. 2010, A&A, 522, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

Appendix A: Parameter distributions

thumbnail Fig. A.1

Posterior densities of orbital periods (Pi), eccentricities (ei), and radial velocity amplitudes (Ki) corresponding to the nine Keplerian signals in the data and of the magnitude of the excess jitter (σj). The solid curve is a Gaussian density with the same mean (μ) and variance (σ2) as the parameter distribution has.

All Tables

Table 1

The relative posterior probabilities of models with k = 0,...,9 Keplerian signals (ℳk) given radial velocities of HD 10180 (or data, d) together with the periods (Ps) of the signals added to the model when increasing the number of signals in the model by one.

Table 2

The six-planet solution of HD 10180 radial velocities.

Table 3

The nine-planet solution of HD 10180 radial velocities.

All Figures

thumbnail Fig. 1

Lomb-Scargle periodograms of the HD 10180 radial velocities for the residuals of the models with six (top) to nine (bottom) periodic signals. The dotted, dashed, and dot-dashed lines indicate the 10%, 1%, and 0.1% FAPs, respectively.

In the text
thumbnail Fig. 2

Phase-folded Keplerian signals of the nine-planet solution with the other eight signals removed.

In the text
thumbnail Fig. 3

Approximated Lagrange stability thresholds between each of the two planets and the MAP orbital parameters of the six-planet solution (Table 2).

In the text
thumbnail Fig. 4

Approximated Lagrange stability thresholds between each of the two planets and the MAP orbital parameters of the nine-planet solution (Table 3).

In the text
thumbnail Fig. 5

Posterior densities of the minimum mass and orbital period of the planet that could exist in the habitable zone of HD 10180 without having been detected using the current data. The solid curve is a Gaussian density with the same mean (μ) and variance (σ2) as the parameter distribution.

In the text
thumbnail Fig. 6

Distributions of the periods and amplitudes received for the three weakest signals in the artificially generated data. The added signals had Pi = 1.18, 9.66, and 67.55 days, and Ki = 0.78, 0.53, and 0.75 m s-1 for i = 1,3,6, respectively.

In the text
thumbnail Fig. A.1

Posterior densities of orbital periods (Pi), eccentricities (ei), and radial velocity amplitudes (Ki) corresponding to the nine Keplerian signals in the data and of the magnitude of the excess jitter (σj). The solid curve is a Gaussian density with the same mean (μ) and variance (σ2) as the parameter distribution has.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.