Bayesian analysis of the radial velocities of HD 11506 reveals another planetary companion

M. Tuomi; S. Kotiranta

doi:10.1051/0004-6361/200811531

Home

All issues

Volume 496 / No 2 (March III 2009)

A&A, 496 2 (2009) L13-L16

Full HTML

Free Access

Issue		A&A Volume 496, Number 2, March III 2009


Page(s)		L13 - L16
Section		Letters
DOI		https://doi.org/10.1051/0004-6361/200811531
Published online		24 February 2009

LETTER TO THE EDITOR

Bayesian analysis of the radial velocities of HD 11506 reveals another planetary companion

M. Tuomi^1,2 - S. Kotiranta²

1 - Department of Mathematics and Statistics, PO Box 68, 00014 University of Helsinki, Finland
2 - Tuorla Observatory, Department of Physics and Astronomy, University of Turku, 21500 Piikkiö, Finland

Received 16 December 2008 / Accepted 14 February 2009

Abstract
Aims. We aim to demonstrate the efficiency of a Bayesian approach in analysing radial velocity data by reanalysing a set of radial velocity measurements.
Methods. We present Bayesian analysis of a recently published set of radial velocity measurements known to contain the signal of one extrasolar planetary candidate, namely, HD 11506. The analysis is conducted using the Markov chain Monte Carlo method and the resulting distributions of orbital parameters are tested by performing direct integration of randomly selected samples with the Bulirsch-Stoer method. The magnitude of the stellar radial velocity variability, known as jitter, is treated as a free parameter with no assumptions about its magnitude.
Results. We show that the orbital parameters of the planet known to be present in the data correspond to a different solution when the jitter is allowed to be a free parameter. We also show evidence of an additional candidate, a 0.8 $M_{\rm Jup}$ planet with period of about 0.5 yr in orbit around HD 11506. This second planet is inferred to be present with a high level of confidence.

Key words: planetary systems - methods: statistical - techniques: radial velocities - stars: individual: HD 11506

1 Introduction

The question of whether an extrasolar planet is detectable or not depends on a delicate mixture of observational technology, effective tools for data analysis, and theoretical understanding on the related phenomena. Traditionally, the instrumentation has been the most celebrated part of this trinity (Santos et al. 2004; Li et al. 2008), but the fainter the signals that one is able to detect, the more the attention should be directed towards the two other factors involved in successful discoveries. It is already known that of poor understanding of the noise-generating physics may produce misleading results (Queloz et al. 2001), although this is also true for the means the data analysis. In many cases the statistical methods involved in the data reduction and analysis of the observations have a tendency to identify local solutions instead of global ones. This is particularly true especially when information is extracted using gradient-based algorithms.

In this letter, we reanalyze the radial velocity (RV) data of a detected extrasolar planet host, namely HD 11506 (Fischer et al. 2007). Our analysis is based on Bayesian model probabilities and full inverse solutions, i.e., the full a posteriori probability densities of model parameters. The Bayesian model probabilities provide a strict mathematical criterion for deciding how many planetary signals are present in the data. We calculate the inverse solutions using the Markov chain Monte Carlo (MCMC) method with the Metropolis-Hastings algorithm (Metropolis et al. 1953; Hastings 1970). These inverse solutions are used to present reliable error estimates for all the model parameters in the form of Bayesian confidence sets. As a result, we present strong evidence for a new and previously unknown planetary companion orbiting HD 11506. We also demonstrate the importance of treating the magnitude of stellar RV noise, the RV jitter, as a free parameter in the model.

2 The model and Bayesian model comparison

When assuming the gravitational interactions between the planetary companions to be negligible, the RV measurements with k such companions can be modelled as (e.g. Green 1985)

$\displaystyle \dot{z}(t)$ = $\displaystyle \sum_{j=1}^{k} K_{j} \big[ \cos ( \nu_{j}(t) + \omega_{j} ) + e_{j} \cos \omega_{j} \big] ,$

(1)

where $\nu_{j}$ is the true anomaly, K_j is the RV semi-amplitude, e_j is the eccentricity, and $\omega_{j}$ is the longitude of pericentre. Index j refers to the jth planetary companion. Hence, the RV signal of the jth companion is fully described using five parameters, K_j, $\omega_{j}$ , e_j, mean anomaly M_0,j, and the orbital period P_j.

Following Gregory (2007b,a,2005a), we calculate the Bayesian model probabilities for different statistical models. These models include the first and second companion Keplerian signals and an additive RV jitter (models M1 and M2, respectively). This jitter is assumed to be Gaussian noise with a zero mean and its deviation, $\sigma_{J}$ , is another free parameter in the model. The jitter cannot be expected to provide an accurate description of the RV noise caused by the stellar surface. It also includes the possible signatures of additional planetary companions if their signals cannot be extracted from the measurements. Hence, this parameter represents the upper limit to the true RV variations at the stellar surface. We also test a model without companions (M0) for comparison.

The instrument error, whose magnitude is usually known reasonably well, is also included as a Gaussian random variable in our analysis. Hence, for a model with k planetary companions, the RV measurements r_i at t_i are described by

$\begin{displaymath} r_{i} = \dot{z}_{k}(t_{i}) + \gamma + \epsilon_{I} + \epsilon_{J} , \end{displaymath}$

(2)

where $\epsilon_{I} \sim N(0, \sigma_{I}^{2})$ is the instrument error, $\epsilon_{J} \sim N(0, \sigma_{J}^{2})$ describes the remaining uncertainty in the measurements, called the RV jitter, and parameter $\gamma$ is a reference velocity parameter.

When deciding whether a signal of a planetary companion has been detected or not, we adopt the Bayesian model selection criterion (Jeffreys 1961). This criterion states that a model with k+1 planetary companions should be used instead of one with k companions if

$\begin{displaymath} P(\dot{z}_{k+1} \vert \vec{r}) \gg P(\dot{z}_{k} \vert \vec{r}) , \end{displaymath}$

(3)

where $\vec{r}$ is a vector consisting of the RV measurements. For more information about Bayesian model comparison, see e.g. Kass & Raftery (1995); Gregory (2005a).

After identifying the most appropriate model in the context of Eq. (3), we must find the Bayesian credibility sets that can assess accurately the uncertainties in the orbital parameters and the RV masses of the planetary companions. The Bayesian credibility sets are robust uncertainty estimates because they show the uncertainty of the model parameters given the measurements. A Bayesian credibility set $\mathcal{D}_{\delta}$ containing a fraction $\delta \in [0,1]$ of the probability density of parameter vector $\vec{u}$ can be defined to a subset of the parameter space U that satisfies the criteria (e.g. Kaipio & Somersalo 2005)

$\begin{displaymath} \left\{ \begin{array}{l} \int_{\vec{u} \in \mathcal{D}_{\d... ... \in \partial \mathcal{D}_{\delta}} = c, \end{array} \right. \end{displaymath}$

(4)

where $p(\vec{u} \vert m)$ is the probability density of the parameters, c is a constant, m is the measurements, and $\partial \mathcal{D}_{\delta}$ is the edge of the set $\mathcal{D}_{\delta}$ . We use $\delta = 0.99$ throughout this article when discussing the parameter errors.

The question of model comparison is not only statistical but also physical. Because the inverse solution is based on non-interacting planets - a simplification that is fairly adequate in most cases - it might be physically impossible in reality. This is important to understand because many extrasolar systems contain large eccentricities that may cause close encounters. In these cases, a detailed description of the dynamics is required. We drew 50 random values of $\vec{u}$ from its posterior probability density and integrated these directly to investigate the orbital evolution. The integration was conducted using the Bulirsch-Stoer integrator (Bulirsch & Stoer 1966), which is famous for its high reliability even when the dynamics include several close encounters. The stability analysis that we used in this analysis is not exact but estimates reliably the expected long-term behaviour from a relatively short ensemble or orbit. In this simple analysis, we identify the variation in the semi-major axis and the eccentricity from randomly selected initial conditions to the end of integration. If the variations are not significant the system is then assumed to be long-term stable (Haghighipour, priv. comm.). The advantage this formulation is the short integration time required when compared to more adequate, for example, Lyapunov exponent based, integration methods.

Table 1: Bayesian model probabilities of 1 and 2 planet models.

The most important error source in this type of study is a phenomenon known as stellar jitter. It is caused by a combination of convection, rotation, and magnetic activity on the stellar surface (see e.g., Wright 2005, and references therein). Although the role of stellar activity as an error source is well known, the magnitude of this error is almost always assumed to be constant, based on certain studies (Wright 2005). Here, we use a more conservative approach and consider the magnitude of the jitter to be free parameter.

Before this work, four exoplanets were found using Bayesian approach instead of a more traditional periodogram, the first of them being HD 73526 c (Gregory 2005a). The candidate HD 208487 c was proposed by Gregory (2005b) and later confirmed by additional observations (Butler et al. 2006; Gregory 2007a). For the system HD 11964, the Bayesian analysis revealed evidence of three planets, despite only one being previously known (Gregory 2007b). None of these works include a discussion of the dynamics, although it is clear that including dynamics in the work enhance the quality of results. This point also is interesting because if the dynamical analysis excludes a part of the parameter space as physically impossible, this restriction can be inserted into Bayesian model as additional a priori constraint that will, with the data, provide tighter confidence limits.

3 Orbital solutions

The star HD 11506 is a quiescent main sequence star of spectral class G0 V. It is a relatively nearby star with a Hipparcos parallax 18.58 mas, which corresponds to a distance of 53.8 pc. It has $T_{\rm eff} = 6060$ K and [Fe/H] = 0.31 (Valenti & Fischer 2005) and its mass is estimated to be 1.19 $M_{\odot}$ (Fischer et al. 2007). We use this mass estimate throughout the paper when calculating planetary masses.

The planetary companion HD 11506 b was first announced by Fischer et al. (2007). They speculated an additional companion could be present because the $\chi^{2}$ value of their single-companion model fit was large (10.3). However, their $\chi^{2}$ value was calculated by assuming a fixed jitter level. Fischer et al. (2007) assumed that $\sigma_{J} = 2.0$ ms^-1. Our solution for this parameter is consistent with this estimate (Table 2).

The Bayesian model probabilities for k=1, 2 are listed in Table 1. These probabilities imply that M2 provides tha most accurate description of the RV data of HD 11506. These probabilities favour two planetary companions, implying that, according to the measurements of Fischer et al. (2007), this star does host two companions. It also verifies that the large $\chi^{2}$ value was consistent with there being a signal of an additional companion, as suspected by Fischer et al. (2007). The corresponding $\chi^{2}$ value of our two-companion solution is 3.1. They also mentioned that a 170 day period could be present in the data but were unable to verify the existence of a second companion. This period corresponds to our two-companion solution (Table 2).

When assuming the fixed jitter level adopted by Fischer et al. (2007), the probability of the one-companion model was again found to be considerably lower than in Table 1 (less than 10^-36). This result emphasizes the fact that the jitter cannot be fixed but must be considered as a free parameter. Furthermore, by fixing the jitter to some a priori estimated value, the probability densities of the orbital parameters become far narrower, which underestimates the uncertainty in the solution.

Table 2: The RV two-planet solution of HD 11506. MAP estimates of the parameters and their $\mathcal{D}_{0.99}$ sets.

In Fig. 1 we show the maximum a posteriori (MAP) orbital solution of model M2 and the velocity curve of HD 11506 c after the signal of b companion has been subtracted. Interestingly, our solution for the companion HD 11506 b differs from that of Fischer et al. (2007). For instance, we found its RV mass to be 3.4 $M_{\rm Jup}$ , whereas they reported a mass of 4.7 $M_{\rm Jup}$ , which is outside the margins of $\mathcal{D}_{0.99}$ set in Table 2. The jitter parameter also appears to have a higher value than estimated. This could be indicative of the difficulties in estimating the jitter magnitude but the jitter may also contain the signal of a third companion.

$\begin{figure} \par\includegraphics[angle=270,width=6.5cm,clip]{1531f1a.ps}\par\vspace*{2mm} \includegraphics[angle=270,width=6.5cm,clip]{1531f1b.ps} \end{figure}$	Figure 1: Radial velocity measurements of HD11506 (Fischer et al. 2008). Two companion MAP orbital solution ( top), and signal of companion HD 11506 c ( bottom).
Open with DEXTER

We were unable to determine any strong correlations between the model parameters, and all of the parameter densities, apart from P₂ (Fig. 2), were reasonably close to Gaussian.

$\begin{figure} \par\includegraphics[angle=270,width=6.5cm,clip]{1531f2.ps} \end{figure}$	Figure 2: The posterior probability density of parameter P₂ and its mode, mean, deviation, and two higher moments. The solid curve is a Gaussian curve with the mean and variance of the density.
Open with DEXTER

We subtracted the MAP one-planet solution from the RV measurements and calculated the Scargle-Lomb periodogram (Lomb 1976; Scargle 1982) for these residuals. The highest peak corresponded to a period of 0.45 yr, which is close to the period of the second companion (0.47 yr). However, the FAP of this peak was as high as 0.54, which ensures that it was impossible to detect the signal of this companion with periodogram. In contrast, this solution was easily found using the MCMC method. In agreement with the MAP solution, we were unable to find other probability maxima in the parameter space of the two-companion model.

4 Orbital stability

To present additional evidence of HD 11506 c, we selected 50 random combinations of parameters taken from the parameter probability densities. These values were used to simulate the dynamical behaviour of the planetary system by direct integration. A random example of these simulations is shown in Fig. 3 (top), where 100 000 yr excerpts are presented for two randomly selected parameter combinations. The orbital ellipses precess slowly but both the semimajor axis and the eccentricities remain almost constant during the evolution. This feature does not prove but suggests strongly that the two planetary companions orbiting HD 11506 provide a physically stable system.

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{1531f3a.eps}\vspace*{2mm} \includegraphics[width=6.5cm,clip]{1531f3b.eps} \end{figure}$	Figure 3: Two random examples of two-planet solution for HD 11506. Green and blue lines indicate orbits of HD 11506 b and c, respectively. Red cross denotes the star. Both figures present the orbit with 100 000 year separation, and the slow precession of the apsid line is clearly visible. Unit of both axes is AU.
Open with DEXTER

We tested the stability further by selecting 50 values of the parameter vector from the region of parameter space most likely prone to instability, i.e., where the quantity Q = a₁(1-e₁)/(a₂(1+e₂)) has its smallest value. We drew the values from a probability density, whose maximum is at $\min Q$ and that decreases linearly to zero at $\max Q$ . This results in a sample from the posterior density that has values mainly close to $\min Q$ . This test was designed to obtain further information about the orbital parameters by excluding parameter values that excluded an unstable system. However, all 50 parameter values resulted in bound orbits even though they precessed strongly. Therefore, it was impossible to extract additional information about the orbital parameters from these simulations. See random example in Fig. 3 (bottom).

5 Conclusions

We have presented complete Bayesian reanalysis of radial velocities of HD 11506, and identified the orbital parameters of a previously unknown exoplanet candidate.

These analyses demonstrated the importance of taking stellar jitter, the most important error source in RV measurements, into account as a free parameter in the analysis. As an unknown noise parameter, jitter cannot be predetermined to some estimated value because of its strong effect on the Bayesian model probabilities of models with different numbers of companions. If overestimated, it may prevent the detection of a additional companions. Its uncertainty must also be taken into account when calculating the error bars for the orbital parameters to prevent the underestimation of their errors.

In the case of HD 11506 the two-planet model is by far more probable for the given data (Fischer et al. 2007) than the original one-planet solution, a result that remained unchanged regardless of whether we considered the jitter as fixed or free. The fit of the two-planet model is similar to any known two-planet signal and a dynamical analysis has demonstrated the solution to be physically possible. We therefore claim that the RV data, which was presented by Fischer et al. (2007), contains the signals of two planetary companions.

Another interesting feature of the full inverse solution of the two-companion model is that the RV mass of HD 11506 b differs significantly from that obtained by Fischer et al. (2007). This difference implies that the Bayesian model selection criterion should be used when assessing the number of planetary signals in RV data. Without accurate knowledge on the best-fit model, the orbital solution can be biased and the resulting statistical conclusions about the properties of extrasolar planetary systems can be misleading.

More observations are required to tighten further the constraints on the parameter space of HD 11506 system.

Acknowledgements

S.K. has obtained funding from Jenny and Antti Wihuri Foundation. The authors would like to thank the referee for valuable suggestions and comments and to acknowledge Proffan Kellari for inspiring environment during the group meetings.

References

Bulirsch, R., & Stoer, J. 1966, Numer. Math., 8, 1 [CrossRef] (In the text)
Butler, R. P., Wright, J. T., Marcy, G. W., et al. 2006, ApJ, 646, 505 [NASA ADS] [CrossRef]
Fischer, D., Vogt, S. S., Marcy, G. W., et al. 2007, ApJ, 669, 1336 [NASA ADS] [CrossRef] (In the text)
Green, R. M. 1985, Spherical Astronomy (Cambridge: Cambridge Univ. Press) (In the text)
Gregory, P. C. 2005a, ApJ, 631, 1198 [NASA ADS] [CrossRef]
Gregory, P. C. 2005b, AIP Conf. Proc., 803 [arXiv:astro-ph/0509412v1] (In the text)
Gregory, P. C. 2007a, MNRAS, 374, 1321 [NASA ADS] [CrossRef]
Gregory, P. C. 2007b, MNRAS, 381, 1607 [NASA ADS] [CrossRef]
Hastings, W. 1970, Biometrika 57, 97 [CrossRef]
Jeffreys, H. 1961, The Theory of Probability (The Oxford Univ. Press) (In the text)
Kaipio, J., & Somersalo E. 2005, Statistical and Computational Inverse Problems, Applied Mathematical Sciences, 160 (In the text)
Kass, R. E., & Raftery, A. E. 1995. J. Am. Stat. Ass., 430, 773
Li, C.-H., Benedick, A. J., Fendel, P., et al. 2008, Nature, 452, 610 [NASA ADS] [CrossRef]
Lomb, N. R. 1976, Ast. phys. Space Sci., 39, 447 [NASA ADS] [CrossRef]
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., et al. 1953, J. Chem. Phys., 21, 1087 [NASA ADS] [CrossRef]
Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
Santos, N. C., Bouchy, F., Mayor, M., et al. 2004, A&A, 426, L19 [NASA ADS] [CrossRef] [EDP Sciences]
Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef]
Valenti, J. A., & Fischer, D. A. 2005. ApJS, 159, 141 (In the text)
Wright, J. T. 2005, PASP, 117, 657 [NASA ADS] [CrossRef] (In the text)

All Tables

Table 1: Bayesian model probabilities of 1 and 2 planet models.

Table 2: The RV two-planet solution of HD 11506. MAP estimates of the parameters and their $\mathcal{D}_{0.99}$ sets.

All Figures

$\begin{figure} \par\includegraphics[angle=270,width=6.5cm,clip]{1531f1a.ps}\par\vspace*{2mm} \includegraphics[angle=270,width=6.5cm,clip]{1531f1b.ps} \end{figure}$	Figure 1: Radial velocity measurements of HD11506 (Fischer et al. 2008). Two companion MAP orbital solution ( top), and signal of companion HD 11506 c ( bottom).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[angle=270,width=6.5cm,clip]{1531f2.ps} \end{figure}$	Figure 2: The posterior probability density of parameter P₂ and its mode, mean, deviation, and two higher moments. The solid curve is a Gaussian curve with the mean and variance of the density.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{1531f3a.eps}\vspace*{2mm} \includegraphics[width=6.5cm,clip]{1531f3b.eps} \end{figure}$	Figure 3: Two random examples of two-planet solution for HD 11506. Green and blue lines indicate orbits of HD 11506 b and c, respectively. Red cross denotes the star. Both figures present the orbit with 100 000 year separation, and the slow precession of the apsid line is clearly visible. Unit of both axes is AU.
Open with DEXTER
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.