Free access
 Issue A&A Volume 494, Number 2, February I 2009 769 - 774 Astronomical instrumentation http://dx.doi.org/10.1051/0004-6361:200810288 11 December 2008

## Determining exoplanet mass with astrometric snapshots

M. Tuomi1,2 - S. Kotiranta2 - M. Kaasalainen1

1 - Department of Mathematics and Statistics, PO Box 68, 00014 University of Helsinki, Finland
2 - Tuorla Observatory, Väisäläntie 20, 21500, Piikkiö, University of Turku, Finland

Received 30 May 2008 / Accepted 3 November 2008

Abstract
Aims. It is commonly assumed that the two indirect exoplanet detection methods, the radial velocity method and astrometric method, require observational periods exceeding the orbital period to produce positive results. Here we test this assumption in detail. We also investigate the smallest ratio of observational timeline and orbital period required for positive detections.
Methods. We obtain full information on the orbital parameters by combining radial-velocity and astrometric measurements by means of Bayesian inference, and sample the parameter probability densities of orbital and other model parameters with a Markov Chain Monte Carlo (MCMC) method in simulated observational scenarios to test the detectability of planets with orbital periods longer than the observational timelines.
Results. We show that, when fitting model parameters simultaneously to measurements from both sources, it is possible to extract much more information from the measurements than when using either source alone. Currently available high-precision measurements of radial velocity (with 1 m s-1 precision) and astrometric measurements achievable with the SIM space telescope (with a precision of 1 as) can be used together to detect a Jupiter analog 30 pc away with an observational timeline of only three years, approximately one fourth of the orbital period. Such measurements are sufficient for determining all its orbital parameters, including inclination and the true mass. Also, with accurate radial velocity measurements covering a timeline of 20 years, the true mass could be determined by astrometric observations within a single year. These case studies demonstrate the potential power of the Bayesian inference of multiple data sources in exoplanet observations. As an example, we show that using the currently available radial velocity measurements, the inclination of HD 154345b could be determined with SIM in a year.

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

## 1 Introduction

Since the discovery of 51 Peg b (Mayor & Queloz 1995), high-precision radial velocity (RV) measurements have been used successfully to detect planetary companions of nearby stars (Butler et al. 2006). The finest instruments to date are capable of achieving an RV precision of 1 m s-1 (Santos et al. 2004; Moorwood & Masanori 2004). However, the exact nature of the companions remains unknown when only RV measurements are available, as these only yield the product of mass and sine of the orbital inclination, giving the lower limit for the mass. Recently, it has been claimed that an RV precision of 1 cm s-1 could be technically possible in the future (Li et al. 2008). However, the RV variations of stars of approximately 5 m s-1 to 50 m s-1 for K to F stars, caused by star spots and irregular convective zones, will likely prevent the detection of signals of the cm scale (Saar & Donahue 1997; Saar et al. 1998).

Astrometric (A) measurements of the position of the target star in the sky as a function of time can be used to detect the inclination, and consequently the true mass of the companion. But despite several trials (van de Kamp 1969; Neuhaeuser et al. 2008), the precision of these measurements has not been high enough to verify the planetary nature of RV companions.

With the aid of telescopes and instruments capable of high-precision astrometry (SIM, GAIA, PRIMA, etc.); however, this situation is about to change. With an estimated astrometric precision of future telescopes of 1.0 as (Unwin et al. 2008) or 8-10 as (Casertano et al. 2008; Derie et al. 2002), it will be possible to determine the masses of the already detected RV companions. It is commonly assumed that these detections require observational periods longer than the orbital period of the target system to be able to detect the periodic signal. Here we test this assumption in detail.

In this article we simulate astrometric and RV measurements to study the possibility of detecting planetary companions of nearby stars with various observational timelines. The goal is to find the minimum timeline required for detecting a planetary companion using high-precision RV and astrometric measurements. These two sources of data are combined by means of Bayesian inference, and the probability densities of orbital and inertial reference frame parameters are sampled using Markov Chain Monte Carlo (MCMC) (Metropolis et al. 1953; Hastings 1970) to find the full global solution to this multidata inverse problem (see Ford 20052006; Maness et al. 2007). Also, the probability densities are sampled to calculate the realistic error bars for parameters and to determine whether a positive detection has been made or not. We also analyse the RV measurements of HD 154345 (Wright et al. 2008) alone and together with simulated astrometric measurements to estimate the minimum observational timeline of astrometry for the detection of the true mass of HD 154345 b.

## 2 Modelling the data

The motion of a planet around a star was treated as a simple two-body system, with masses and for the star and the planetary companion, respectively. In Cartesian coordinates, when the gravitational forces between the possible other planets in the system are assumed negligible, the column position vector of the star with respect to the barycentre of the system can be expressed as a function of time (t) as (e.g. Green 1985)

 (1)

where E is the eccentric anomaly, satisfying the Kepler equation , and n = P-1 is the orbital frequency, P the orbital period and t0 the time of periastron. The velocity and position and w.r.t. the observer are some constant vectors defining the inertial reference frame. Vectors and are constant vectors, defined as

 (2)

The parameters , e, , , and are the orbital parameters of the system: semimajor axis of the star, eccentricity, longitude of ascending node, cosine of the inclination and the longitude of pericentre, respectively. The semimajor axis of the star can be expressed as a function of the masses of the gravitationally interacting bodies,

 (3)

The models for astrometric and RV data are now simply for astrometry and for RV. Here D is the distance of the system from the observer. There are now 12 independent parameters describing the system. The parameter vector of parameter space U can be written as , where , , , and . To fully describe the system, the masses and should be treated as independent unknown parameters, and the distance D should also be included in vector . Here we assume that and D are known with sufficient accuracy by some astrophysical techniques.

It is essential to include the parameters defining the inertial reference frame to be able to fully investigate the probability densities of the parameters. It is also assumed that the model parameters do not change as a function of time during the observational timeline.

It is intuitively clear that the availability of both RV and A data sets, instead of only one of them, should increase the amount and the quality of information on the observed system. This should be true even if one of the data sets is significantly more inaccurate than the other. In what follows, we state this principle more rigorously and present a practical method for combining the two data sources. In fact, since RV and A data are radial and tangential projection samples of orbital motion, they are separate sources only in the instrumental sense, so we could just as easily talk about full velocity data.

## 3 Simulations and fitting procedure

We simulated data sets to study the complementarity of astrometric and RV time series. Assuming that the errors for each data point i = 1,...,N are independent and identically distributed and that their probability density functions (PDF's) are Gaussian ( ), the observed RV data are modelled as

 (4)

Equally, the two-dimensional astrometric model can be written as

 (5)

To avoid inversion crimes, i.e. to avoid using exactly the same model to generate the measurements and to find the inverse solution (see e.g. Kaipio & Somersalo 2005), an additional planet of low mass was included in the model when generating the data. This choice was made because small systematic errors make the simulated measurements more realistic. Otherwise the simulated measurements and the corresponding solutions would only help for studying the model in Eq. (1), not necessarily situations encountered in reality.

Several data sets were generated, each with a different value for the data parameter T, representing the length of the observational timeline. It was further assumed that the observation time ti was evenly distributed in T.

When simulating the measurements, the timelines were set to , 10.0, 5.0, 3.0, 1.5, 1.2, 1.0 and 0.8 years. The values of all the other data and orbital parameters and masses were fixed in all these scenarios. These values were set to (1.0 m/s, 1.0 as, 100, 100) and (5.0 AU, 0.1, 1.0, 1.0, 1.0, 1000.0d, , , 30 pc), where is the mass of Jupiter. This simulated system is called S1.

### 3.1 What is a positive detection?

In the simplest possible case, when e = 0 and , the detection threshold of full velocity data can be calculated analytically. Let N be the number of data points, T the length of the timeline of observations, and the standard deviation of observations. Following the approach in Eisner & Kulkarni (2001a, b), the detected signal of the velocity variation amplitude is a false one produced by the uncertainties in the data with a probability <1% if

 (6)

where
 (7) (8) (9) (10)

and is in radians. This approach excludes the uncertainties in the orbital period and can therefore only yield the lower limit for the detection threshold. Hence, if Eq. (6) does not hold, it will be impossible to detect the signal. However, if it holds, the detectability of such a companion needs to be examined more closely by numerical simulations and by analysing the simulated data using methods such as MCMC and Bayesian model selection criterion.

To fully investigate the ability to detect planetary companions, we must define when a positive detection has been made. This question can be approached through Bayesian probabilities. Let be the model in Eq. (1) with one planetary companion (corresponding 12 parameters in the RV and astrometry models), and a model without a planetary companion (5 parameters). In general, let be a model with k planets.

Using the Bayes theorem, it can be seen that the conditional probability of model representing the data (m) best, out of the p+1 alternatives to be tested, can be written as

 (11)

where the Bayes factor Bk,j is defined as (e.g. Kass & Raftery 1995)

 (12)

and is the prior probability of the kth model, here set equal for all k, because it is assumed that there is no prior information available. Here the likelihood , with parameters for the kth model, is

 (13)

where is the parameter likelihood function and the prior density.

Since the model probability, defined in this way, automatically takes the Occamian principle of parsimony into account, the model with the smallest number of parameters out of those having almost equal probabilities will be selected. Hence, it can be said that a detection has been made if (Jeffreys 1961)

 (14)

This criterion is used throughout this article when deciding whether a statistically significant detection has been made or not.

### 3.2 Fitting method

The fitting was performed by requiring that the values of all the least-squares cost-functions Sx (astrometric x), Sy (astrometric y), (RV), and their sum be minimized simultaneously. This method, called multidata inversion, has been used successfully with astrometric and RV measurements when detecting stellar binaries (e.g. Torres 2007). See the discussion in Kaasalainen & Lamberg (2006), where the multidata inversion was applied to asteroid observations.

The models for astrometric position and RV of the two-body system of interest are non-linear, so an iterative method of fitting the model parameters is needed. The MCMC with Metropolis-Hastings (M-H) algorithm was chosen because it is a global method (Metropolis et al. 1953; Hastings 1970), it offers a direct estimate of the posterior probability density, and because it can be used to verify the existence and uniqueness of the solution. Since the probability densities given the measurements are available, MCMC can be used to calculate realistical error estimates for the model parameters. These estimates are typically much larger than those calculated using traditional methods (e.g. Ford 2006), implying that MCMC should be preferred when assessing the parameter errors. Assuming Gaussian errors with zero mean, the likelihood function of the parameters with respect to RV measurements can be written as

 (15)

When applying MCMC, a parameter value ( ) is selected for the first member of the chain. The next value is found by randomly selecting a proposal in the vicinity of . This is then accepted by comparing the likelihoods of the two parameter values. Proposed parameter values  with a greater likelihood than that of are always selected as the next chain member, but values with a smaller likelihood can also be selected according to the criterion of Hastings (1970). Samples of at least 105 points were generated when sampling the parameter space. For practical details on MCMC with astronomical data, see e.g. Gregory (2005).

The parameter space U in this Keplerian two-body model has a comparatively small dimension ( ), but in some cases it already makes the sampling computationally expensive. Especially when covariances between the parameters are large and of non-linear nature, the space of reasonable probability to be sampled can be very narrow and, as a result, the next proposed value of parameter vector in the Markov chain is likely to be outside this subspace and thus rejected, considerably increasing the time needed to generate a statistically representative chain. For this reason, when using a multivariate Gaussian density as a proposal, the acceptance rates were low, approximately 0.1 in the MCMC samplings.

## 4 Taking advantage of Bayesian inference

With more than one source of measurements available, it is possible to get more information from the system of interest than when relying on any single observation method alone. This is a consequence of Bayesian inference.

Denoting the astrometric measurements by and the RV measurements by , the conditional probability of having parameter vector , is a product of the impacts the two sets of measurements have on this hypothesis:

 (16)

Thus, there is always more information available on the system - either in a narrower parameter density or in the possibility of including more parameters in the model - when using multiple data sources. This is due to complementary rather than just additional information: the separate probability densities for complementary sources are quite different from each other, so their product (joint probability) is much more tightly bound than either factor alone.

Equation (16) is in fact just another way of stating that we simply minimized the sum , as can be seen by applying Eq. (15), while checking that each Si was still close to their minima. However, when calculating model probabilities and parameter densities, the formulation in Eq. (16) has to be used.

### 4.1 Correlations and complementarity

By having measurements made using two different observational techniques has an effect on the parameter PDF's. This happens because the two measurements are modelled using a different model with differing parameters for the inertial reference frame. Therefore, it is expected that the two measurements contain complementary information on different aspects of the system.

Generally, correlations between model parameters occur if, for some small displacement of the parameter vector, , the model used to describe the system satisfies

 (17)

Now the parameter PDF's are broadened or correlated until Eq. (17) no longer holds.

In a stellar system with a single planetary companion, the most obvious possible coupling, a positive correlation between t0 and , is a natural byproduct of the two-body Keplerian model. When assuming and using Eq. (2), Eq. (1) becomes

 (18)

Setting implies , resulting in a positive linear correlation between parameters and t0.

If it is also assumed that the observational timeline is much shorter than the orbital period, , more correlations take place in this long-period system. As for all i = 1, ..., N, when , each as well. Thus, this assumption justifies and and Eqs. (4) and (5) become

 (19)

where and are just the functions presented in Eq. (2) with the angle replaced with , and matrix .

Clearly, Eq. (17) holds if and . This means that it is possible to change the values of the components of vector by any amount and a corresponding negative change in the components of vector cancels this change exactly. As a result, the components of these vectors can correlate negatively. Also, the components of can correlate similarly with the components of . For RV, can correlate with , but the product has no corresponding parameter to correlate with. These are just the correlations described by Eisner & Kulkarni (2001a,b, 2002).

From Eq. (19) it is also clear that the orbital frequency can correlate with and making the detection of planetary signal harder and broadening the densities of the corresponding parameters.

Despite the existence of and due to the partially complementary nature of these correlations, it is possible to detect the periodic signal of a long-period planetary companion in the Bayesian model selection sense. The equiprobability contours of parameter combinations and with the simulated system S1 are shown in Fig. 1 for  years, which is approximately one fourth of the orbital period. The Bayesian model probabilities were found to satisfy the condition of Eq. (14), and the contours in Fig. 1 demonstrate that it is indeed possible to detect extrasolar planetary companions even if the observational timeline is shorter than the orbital period. Also, since clearly I < 1, the planetary nature of these companions can be verified in this scenario.

 Figure 1: Equiprobability contours containing 99%, 95%, 90%, and 50% of parameter PDF's showing the densities of and correlations between parameters and n and parameters I and . The simulated system has a long-period planet with (  years). Open with DEXTER

When using the two data sources, the planetary signal could not be detected for  years. For astrometric or RV measurements alone, this signal was found to be undetectable for  years, which clearly demonstrates the advantages of the Bayesian inference of multiple datasets.

### 4.2 Astrometric snapshots and detection thresholds

Astrometric observations with the property are called astrometric snapshots. This definition is made because the astrometric observations are now made in a fraction of the time interval of the RV observations incapable of separating m and . We modified the simulated system S1 by fixing  years and denoting this by S2.

Using Bayesian inference between RV and astrometric measurements made it possible to fit all the parameters in the model, including I and , to the measurements, even when the signal of the planetary companion could not be detected using astrometric observations alone. For the simulated system S2, the condition in Eq. (14) was satisfied for values of as low as 1.0 years, which is less than one tenth of the orbital period. The corresponding densities of I and in this scenario are shown in Fig. 2. It can be seen that the maximum a posteriori estimates of these densities are very close to the values selected for the simulated system. This implies that the true mass of the stellar companion is obtainable. Also, the density of I shows that certainly I < 1, implying that a companion of planetary mass has indeed been detected as claimed. The simulated RV and astrometric measurements and the maximum a posteriori orbit are shown in Fig. 3 for this snapshot scenario.

 Figure 2: PDF's of parameters I and in a snapshot scenario with ( years). The mode, mean (), standard deviation (), skewness (), and kurtosis () of the densities are shown. The solid curve is a Gaussian function . Open with DEXTER

 Figure 3: Simulated RV and astrometric measurements and the maximum a posteriori orbit. Open with DEXTER

The reason the parameters I and can be fitted is that the parameters in the RV model are now well-constrained by the RV measurements. The only possibility for Eq. (17) to be true is that , the amplitude of RV variations, remains unaltered even if and i do not. This implies that , which is the correlation observed in the parameter densities of and I (Fig. 4). Because of this correlation, the reference frame parameters of astrometry can correlate freely with . This is also demonstrated in Fig. 4.

 Figure 4: Equiprobability contours containing 99%, 95%, 90%, and 50% of parameter probability densities. The correlations between and I and of and in a snapshot scenario with (  years). Open with DEXTER

## 5 Case study: HD 154345b

Recently, Wright et al. (2008, hereafter W08) reported a detection of a Jupiter analog orbiting a G8 dwarf HD 154345. They claim that astrometric measurements over its 9-year period would determine the orientation of the orbital plane and as a consequence the true mass. The data published in W08 was re-examined and the orbital solution found using MCMC. These data have 55 measurements over a period of 10.4 years. The largest gap between two subsequent observations within these measurements is 352 days. The orbital parameters were calculated assuming the same jitter level as in W08, 2.5 m s-1, and are listed in Table 1. This table shows the MAP estimates of the parameters and their 99% confidence sets (CS). Missing confidence sets indicate that the posterior density of the corresponding parameter has significant values everywhere in its parameter space. The results of Wright et al. (2008) with 99% confidence limits are shown for comparison. The corresponding fit is shown in Fig. 5. The large uncertainties of parameters and t0 (their 99% Bayesian confidence sets are equal to their parameter spaces) stem from these parameters being meaningless for circular orbits.

 Figure 5: RV measurements of HD 154345 and the maximum a posteriori orbit of the planetary companion. Open with DEXTER

Regardless of the fact that in W08 the RV movement of the host star had been subtracted from the data, the constant movement parameter was fitted as well to be able to take its uncertainty and its effect on the orbital parameters into account. The orbital parameters were found consistent with but their confidence intervals smaller than those reported in W08. This difference in the confidence intervals takes place likely because the Bootstrap method was used to assess the parameter errors in W08 instead of direct sampling of the posterior density.

Table 1:   The solution and error estimates of the HD 154345 system.

Astrometric measurements with as, , and , 5.0, 3.0, 2.0, 1.5, 1.2, 1.0, and 0.8 years were generated to study the detectability of parameters I and . These measurements were generated assuming that there is a planetary companion with orbital parameters in Table 1, and  rad). This simulated data was used together with the real RV measurements published in W08 to find the limiting for which the true mass of the planet could still be measured with the SIM telescope. The selection I=0 was made because changing this value would result in a higher planetary mass and hence in a stronger astrometric signal, making the detection of orbital plane parameters even easier. Regardless of large error bars, we found that it is possible to detect the orbital plane parameters with years. With this short timeline, the 99% error bars of the parameters were [-0.30, 0.22] and [0.39, 1.54] for I and , respectively, demonstrating that it was indeed possible to determine their values.

## 6 Conclusions and discussion

The time needed to make a positive detection of an extrasolar planetary conpanion candidate depends essentially on its orbital period. It is commonly assumed that, to be able to detect the signature of such companion, an observational timeline longer than the orbital period is required. Also, since most of the exoplanet candidates have been detected using the RV method, only the lower limit of their mass is available. With the aid of future space telescopes and accurate astrometric measurements, it will be possible to detect the inclination and thus the true mass of planetary candidates.

We have shown that when high-precision RV and accurate astrometric measurements are both available, it is possible to detect the true mass of stellar companions with observational timelines considerably shorter than their orbital periods. Also, when the RV measurements have a long time span, astrometric measurements can reveal the true mass of a stellar companion in less time than one tenth of the orbital period of the system. This ability is also demonstrated using the RV measurements of HD 154345 as an example. We find that, having these measurements with  years in hand, astrometric observations with SIM telescope are sufficient for obtaining the true mass, within a single year.

Bayesian inference plays an important role when extracting information from several sources of measurements. The ability to use RV and astrometric measurements simultaneously makes it possible to employ observational timelines below the orbital ones and still be able to make positive exoplanet detections, thus helping to extract the maximum amount of information from measurements and increasing the time efficiency of observations.

In a forthcoming study, we plan to study the inclusion of additional transit-photometry measurements to further tighten the parameter probability densities in transiting scenarios. Also, the approach used here should be extended to systems with two or more planetary companions.

Acknowledgements

S. Kotiranta was supported by the Jenny & Antti Wihuri foundation, and M. Kaasalainen by the Academy of Finland (project New mathematical methods in planetary and galactic research''). We would like to thank Dr. J. Wright and the two anonymous referees for pointing out errors and inaccuracies in the manuscript.

## All Tables

Table 1:   The solution and error estimates of the HD 154345 system.

## All Figures

 Figure 1: Equiprobability contours containing 99%, 95%, 90%, and 50% of parameter PDF's showing the densities of and correlations between parameters and n and parameters I and . The simulated system has a long-period planet with (  years). Open with DEXTER In the text

 Figure 2: PDF's of parameters I and in a snapshot scenario with ( years). The mode, mean (), standard deviation (), skewness (), and kurtosis () of the densities are shown. The solid curve is a Gaussian function . Open with DEXTER In the text

 Figure 3: Simulated RV and astrometric measurements and the maximum a posteriori orbit. Open with DEXTER In the text

 Figure 4: Equiprobability contours containing 99%, 95%, 90%, and 50% of parameter probability densities. The correlations between and I and of and in a snapshot scenario with (  years). Open with DEXTER In the text

 Figure 5: RV measurements of HD 154345 and the maximum a posteriori orbit of the planetary companion. Open with DEXTER In the text

Copyright ESO 2009