Analytical determination of orbital elements using Fourier analysis
I. The radial velocity case
^{1} Observatoire de l’Université de Genève, 51 chemin des Maillettes, 1290 Sauverny, Switzerland
email: jeanbaptiste.delisle@unige.ch
^{2} ASD, IMCCECNRS UMR 8028, Observatoire de Paris, UPMC, 77 Av. DenfertRochereau, 75014 Paris, France
Received: 10 December 2015
Accepted: 23 March 2016
We describe an analytical method for computing the orbital parameters of a planet from the periodogram of a radial velocity signal. The method is very efficient and provides a good approximation of the orbital parameters. The accuracy is mainly limited by the accuracy of the computation of the Fourier decomposition of the signal which is sensitive to sampling and noise. Our method is complementary with more accurate (and more expensive in computer time) numerical algorithms (e.g. LevenbergMarquardt, Markov chain Monte Carlo, genetic algorithms). Indeed, the analytical approximation can be used as an initial condition to accelerate the convergence of these numerical methods. Our method can be applied iteratively to search for multiple planets in the same system.
Key words: celestial mechanics / planets and satellites: general / methods: analytical / techniques: radial velocities
© ESO, 2016
1. Introduction
In this article we are interested in retrieving planetary orbital parameters from a radial velocity time series. These parameters are usually obtained using numerical leastsquares minimization methods (LevenbergMarquardt, Markov chain Monte Carlo, genetic algorithms, etc.). These methods enable one to explore the parameter space and find the best fitting solution as well as estimates of the error made in the parameters. However, numerical methods need to be initialized with starting values for the parameters. When initial values are far from the solution, these algorithms can become very expensive in terms of computation time, or even unable to converge (nonlinear fit can have several local χ_{2} minima, etc.). On the contrary, using a good guess of the parameters as initial conditions significantly improves the efficiency of numerical methods.
Here we describe an analytical method that allows a very efficient determination of the orbital parameters. This method is based on the Fourier decomposition of the radial velocity data, which is obtained using linear leastsquares spectral analysis. This idea of recovering the orbital parameters of the planet from the Fourier decomposition of the radial velocity signal has already been proposed by Correia (2008) and used in the analysis of different planetary systems (see Correia et al. 2005, 2008, 2009, 2010). A similar idea was also developed in the case of interferometric astrometry by Konacki et al. (2002). In this article, we complete the sketch proposed by Correia (2008) and provide a fully analytical method to retrieve the orbital parameters. In addition to this Fourier method, we propose an alternative algorithm that is based on information contained in the extrema of the radial velocity curve and that can be used when the Fourier decomposition is unreliable (e.g. illsampled, very eccentric planets). These algorithms are implemented in the DACE web platform (see Buchschacher et al. 2015) ^{1}.
In Sect. 2 we derive an analytical decomposition of a Keplerian radial velocity signal in Fourier series. In Sect. 3, we describe an analytical method to compute orbital parameters from the Fourier coefficients of the fundamental and first harmonics of a Keplerian signal. In Sect. 4, we show how to find the fundamental frequency and compute the Fourier coefficients. In Sect. 5, we illustrate the performances and limitations of our methods on observed planetary systems. In Sect. 6, we summarize and discuss our results.
2. Fourier decomposition of a Keplerian radial velocity signal
We assume here that the observed system is composed of a star and a planet (no perturbation). The motion of the star with respect to the center of mass is Keplerian. The radial velocity of the star in the reference frame of the center of mass reads (1)with (2)and a is the semimajor axis of the planet with respect to the star, e the eccentricity, i the inclination, ν the true anomaly, ω the argument of periastron, P the orbital period, m_{p} the planet mass, m_{s} the star mass. We denote by n the meanmotion of the planet: (3)with (4)and is the gravitational constant. The radial velocity signal is Pperiodic and can be decomposed in discrete Fourier series: (5)with (6)where M is the mean anomaly and M_{0} is the mean anomaly at the reference time t = 0. The coefficients V_{k} are complex numbers. Since V(t) is real, we have (where is the complex conjugate of V_{k}). We denote by the Hansen coefficients (7)These coefficients are real numbers that only depend on the eccentricity (e), and we have . The Fourier expansion of the radial velocity can be rewritten using Hansen coefficients We observe that the radial velocity only depends on Hansen coefficients of the form (k ∈ Z^{∗}). In the following, we drop the exponents and use the notation .
We note that V_{0} = 0 is only valid if the observer is fixed in the reference frame of the center of mass. If we account for a constant motion of the observed system with respect to the solar system, we have V_{0} ≠ 0 but all other coefficients are unaffected.
3. From Fourier coefficients to orbital parameters
Let us assume that we are able to determine the period of the planet as well as the values of V_{1} (fundamental) and V_{2} (first harmonics) from an observed radial velocity signal. We thus have access to five observables (because V_{k} are complex numbers). Therefore, this information should be sufficient to determine the five parameters: P, K, e, ω, M_{0}.
3.1. Expanding Hansen coefficients
Hansen coefficients are functions of the eccentricity alone and can be expanded in power series of the eccentricity. Since we restrict our study to the fundamental and the first harmonics, we only need the expansion of X_{± 1} and X_{± 2}. We have (see Appendix A)
3.2. Eccentricity and phase
Since we assume that the complex coefficients V_{1}, V_{2} are known, the ratio ρ ≡ V_{2}/V_{1} is also known. From Eq. (9), we have (14)From Eqs. (10)−(13) we obtain We thus have (18)and (19)We introduce the complex coefficient (20)and use, in the following, the approximation (21)where order 5 (and more) terms are neglected. We first suppose that C(ω) is known and show how to determine the eccentricity and the angle M_{0}. Using approximation (21), we only have to solve a third order polynomial equation to obtain the eccentricity (22)where ℜ(C) denotes the real part of C, (23)Since we always have ℜ(C) < 1/3, there is at most one solution for e in the interval [0,1], and only one, provided that  ρ  ≤ 1−ℜ(C). One can easily verify that this solution is given by (24)where the hat denotes an estimation of the considered parameter (here the eccentricity e). Using the same approximation (see Eq. (21)), the angle M_{0} is estimated by (25)These estimates (Eqs. (24) and (25)) are based on the assumption that C(ω) is known. As a first (very crude) approximation, we can use C(ω) ≈ 1/4, since C varies in a small disk in the complex plane around 1/4 (see Eq. (20)). However a better approximation is given by first computing an estimate of ω and using Eq. (20) to obtain C. We introduce (26)which enables us to eliminate M_{0}. We know that (see Eqs. (16), (17)), thus, at leading order in eccentricity (27)and the phase of η provides a very crude estimate of ω(28)(where the double hat stands for a very crude estimate) which can be used to compute the complex coefficient C (Eq. (20)), and the estimates ê, .
3.3. Other parameters
The remaining parameters (K and ω) are easily obtained now that we have an estimation for the eccentricity ê and the mean anomaly at the reference time (). We only consider the coefficient of the fundamental (29)and rewrite it has a function of Kcosω, Ksinω(30)Therefore We obtain estimates for K and ω by introducing the estimates of e and M_{0} in Eqs. (31), (32): The Hansen coefficients X_{± 1} can be evaluated at ê using the power series expansions of Eqs. (10), (12) or by using numerical estimations of the integral of Eq. (7) (see Appendix A). In the case of vanishing eccentricity, M_{0} and ω are meaningless quantities, while the mean longitude at reference time λ_{0} = M_{0} + ω is the only important angle. Our analytical method provides an accurate estimate for λ_{0} even in the low eccentricity case. Indeed, at low eccentricity, the estimate (see Eq. (25)) may be incorrect because the coefficient V_{2} is small and more sensitive to noise. However, for small e, we have X_{1}(e) ≈ 1 and X_{1}(e) ≪ X_{1}(e), and thus (35)The phase of V_{1} directly provides an estimate of λ_{0}. Still assuming small eccentricity, we have (see Eqs. (33), (34)) (36)Therefore, is also incorrect, but the estimate is not affected: (37)In principle, higher order harmonics could be used to check and/or improve the solution. In particular, the second harmonics could be used to discriminate between two planets in a 2:1 resonance and a single eccentric planet. Indeed, if the second harmonics coefficient is not compatible with the orbital parameters that were determined from the fundamental and first harmonics, it could be evidence of the presence of a second planet that is in resonance. However, the amplitude of the kth harmonics is of order Ke^{k}. Thus the Fourier coefficients of higher order harmonics are much more sensitive to the noise and less reliable, especially for low eccentricities.
3.4. Refinement of estimates
Our analytical determination of the parameters gives very good results for small to moderate eccentricities (see Sect. 3.5). However, at high eccentricities, the estimates exhibit significant divergences from the true parameters (while still being in their neighbourhood). To refine the estimates, one can use a simple NewtonRaphson method. We denote by x the vector of parameters and y the vector of Fourier coefficients: (38)We denote by J the Jacobian matrix of our system (39)The NewtonRaphson method provides corrections for the estimates of the parameters x(40)From our estimates of the parameters, we can compute estimates for V_{1}, V_{2} (see Eq. (9)), and thus determine the ŷ vector and compare it to the known Fourier coefficients (vector y). The estimation of V_{1}, V_{2} requires estimates for the Hansen coefficients X_{± 1}, X_{± 2} which can be obtained analytically or numerically (see Appendix A). The computation of the matrix J requires estimates of the partial derivatives of V_{1}, V_{2}. From Eq. (9) we deduce We thus need to determine estimates of the derivatives of Hansen coefficients (for k = −2,−1,1,2). We use the same method as for Hansen coefficients computation (see Appendix B). Then, we only need to invert the 4 × 4 matrix J to obtain the corrections for the orbital parameters. This process can be reiterated to increase the accuracy. This numerical method converges very rapidly since we already have good approximations for the parameters.
3.5. Accuracy of estimates
To test the accuracy of our algorithm, we computed the Fourier coefficients V_{1} and V_{2} according to Eq. (9) for different values of e and ω and using a numerical computation of the Hansen coefficients (see Appendix A). We arbitrarily set M_{0} = 0, and K = 1 because the value of these parameters should not affect the accuracy of the algorithm. We then applied our algorithm to obtain estimates of K, e, ω, M_{0} and compared them to the expected values. Figure 1 shows the results of this comparison, using the estimate C ≈ 1/4, or using (see Sect. 3.2). In the latter case, we also show the results of the NewtonRaphson refinement after one and two iterations. We observe that for e = 0.5, all methods approximate the orbital parameters very well (see Fig. 1, left). Above 0.8, the estimate using C ≈ 1/4 does not constrain the eccentricity and the argument of periastron well, while the use of improves the results significantly (see Fig. 1). Finally, we see that, with only two iterations of the NewtonRaphson correction, we obtain very accurate parameters, even for e = 0.95 (see Fig. 1, right). Therefore, for applications to real data (see Sect. 5) we systematically use these two steps of NewtonRaphson correction in our algorithm. These encouraging results are obtained assuming that the Fourier coefficients V_{1}, V_{2} have been determined very accurately. In real cases the errors might be much greater if the estimates of V_{1}, V_{2} are not precise enough.
Fig. 1 Comparison between our analytical estimates (see Eqs. (24), (25), (33), (34)) and the expected values of the eccentricity (e), the mean longitude at the reference time (λ = M_{0} + ω), the semiamplitude K, and the argument of periastron (ω) (from top to bottom), as a function of ω and for e = 0.5, 0.8, 0.9, 0.95 (from left to right). We show the errors obtained on the parameters using the approximations C = 1/4, and (see Sect. 3.2). For the approximation , we also show the errors after one and two iterations of the NewtonRaphson correcting algorithm (see Sect. 3.4). 

Open with DEXTER 
4. Period finding and Fourier coefficients
Our estimation of the planetary orbital parameters is based on the assumption that we are able to determine the planet period and the Fourier coefficients V_{1}, V_{2} (see Eq. (6)) from the radial velocity signal. The classical method for this is to use a leastsquare spectral analysis (such as the LombScargle periodogram, e.g. Zechmeister & Kürster 2009). We denote by t_{i}, v_{i}, and σ_{i} (i ∈ [1,N_{m}]) the time, the value of radial velocity, and the estimated (Gaussian) noise of each measurement. We assume that (45)with ϵ_{i} randomly drawn from a normal distribution . We aimed to decompose the radial velocity signal on a basis of functions f_{k}(t) (k ∈ [1,N_{f}]), such that (46)To find the coefficients C_{k}, we perform a weighted leastsquares regression. We introduce the matrices (47)(48)The leastsquares regression provides an estimator for C(49)The functions f_{k} that are used in the decomposition basis are usually the sine and cosine of a given period: In addition to these sinusoidal functions, one can introduce shifts for each instrument (which also account for the proper motion of the system with respect to the solar system), as well as power law drifts (t^{k}, k = 1,2,3...) to better model the signal.
We then compute the squared difference between the model and the data (χ^{2}(P) which depends on the chosen period P). The bestfitting period is found by minimizing the χ^{2}(P) function. One can construct the periodogram of the signal by computing the power for a range of periods (e.g. Zechmeister & Kürster 2009). The power is simply a renormalization of the χ^{2} function (see Zechmeister & Kürster 2009): (52)where is the squared difference between the data and the model without the periodic functions (only shifts and drifts).
We assume that the period of the planet corresponds to a peak in the periodogram (which is equivalent to a minimum of χ^{2}). We denote by the period that corresponds with this peak. We now need estimates for V_{1} and V_{2}. We obtain these by using the same leastsquares algorithm, but we introduce in the function basis f the sines and cosines corresponding to both the fundamental period and the first harmonics : Again, we can complete this basis with instrumental shifts and drifts and we obtain estimates for the coefficients C_{k} (see Eq. (49)). By definition, V_{1} and V_{2} are then given by and the method described in Sect. 3 can be used to retrieve the orbital parameters of the planet from these coefficients.
As shown in Sect. 3.5 (Fig. 1), our estimates of orbital parameters are very accurate if the Fourier coefficients are perfectly determined. In real cases, the accuracy of Fourier coefficients is affected by the noise and the sampling of the data. The determination of Fourier coefficients is thus the main source of error in our algorithm. In particular, for very eccentric planets, if the first harmonics amplitude is overestimated (compared to the fundamental amplitude), there may not be any meaningful solution for the orbital parameters (especially the eccentricity), and our algorithm cannot be used. In this case, we can use an alternative method to find estimates for the orbital elements. In Appendix C, we describe a method that makes use of the values and timings of the minimum and maximum of the (folded) radial velocity curve to provide rough estimates of the orbital parameters (to use when the Fourier method breaks).
5. Illustrations and performances
In this section, we apply our method to actual radial velocity data to illustrate its accuracy and limitations. All the algorithms described in this article have been implemented in the DACE platform (http://dace.unige.ch, see Buchschacher et al. 2015). The DACE platform also implements a LevenbergMarquardt algorithm that enables us to refine the solution initially found with our methods. We denote by Fourier fit (FF; see Sect. 3), min/max fit (MMF; see Appendix C), and LevenbergMarquardt (LM) the different algorithms we use to fit the data and find the orbital elements. All the following results are obtained using the DACE platform. We selected four representative cases: GJ 3021 (one eccentric planet), HD 156846 (one illsampled, very eccentric planet), HD 192310 (two planets with a low S/N), and HD 147018 (two planets with equal semiamplitudes). For all these analyses, we quadratically added an instrumental systematic error to the radial velocity uncertainty. For CORALIE98, this instrumental error is set to 6 m s^{1}, for CORALIE07 we use 5 m s^{1}, and for HARPS 0.75 m s^{1}.
We note that our FF method provides values for K, e, ω, M_{0}, but other sets of parameters can then be derived and used for the LM fit. In particular, for low eccentricities and/or low S/N, the set K, k = ecosω, h = esinω, λ_{0} = M_{0} + ω is better suited. Indeed, for low eccentricities, M_{0} and ω present strong anticorrelations since the periastron direction is illdefined. Using k, h, λ_{0} allows to avoid this issue. In the following illustrations, we used the same set of parameters as the published ones to be able to compare our results with published solutions.
5.1. GJ 3021
Fig. 2 Radialvelocity data of GJ 3021 (top) and corresponding periodogram (bottom). The data are taken from Naef et al. (2001; 61 CORALIE98 measurements). 

Open with DEXTER 
Orbital parameters of GJ 3021b obtained using our FF algorithm, using the LM refinement algorithm (starting from the FF solution), and compared with the published solution (PUB, Naef et al. 2001, on the same 61 CORALIE98 measurements).
Fig. 3 Radialvelocity data of GJ 3021 superimposed with the Keplerian model obtained with the FF algorithm (top), and the model obtained with the LM algorithm starting from the FF solution (bottom). The orbital elements corresponding to both solutions are presented in Table 1, and compared to the published solution (Naef et al. 2001). 

Open with DEXTER 
GJ 3021b is an eccentric (e ≈ 0.5) planet discovered by the CORALIE survey (see Naef et al. 2001). We applied our method (using DACE) to the exact same data as Naef et al. (2001), which consist of 61 CORALIE98 measurements. Figure 2 shows these data and the corresponding periodogram. From this periodogram we selected the period that corresponds to the highest peak (about 130 days) and applied our FF algorithm. We then used this FF solution as an initial step for an LM algorithm to refine our solution. The orbital parameters provided by both algorithms (FF and LM) are shown in Table 1, and the corresponding radial velocity curves are shown (superimposed with the data) in Fig. 3.
The LM solution is almost identical to the published solution and the differences are well below the error bars (see Table 1). The FF solution exhibits more significant differences but remains a very good approximation (see Table 1), especially considering the fact that only one period is well sampled (see Fig. 3).
The case of GJ 3021 shows that the FF solution can provide a very good approximation of the orbital parameters even for quite high eccentricities (e ≈ 0.5). This approximation can then be refined with numerical methods (such as the LM algorithm), which converge much more rapidly to the solution than without this first approximation.
5.2. HD 156846
Fig. 4 Same as Fig. 2 but for HD 156846. The data are taken from Tamuz et al. (2008). They consist of 53 CORALIE98 measurements (red), and 11 CORALIE07 measurements (blue). 

Open with DEXTER 
Fig. 5 Phasefolded radialvelocity curve of HD 156846b superimposed with the Keplerian model that was obtained with the MMF algorithm (top), and the model obtained with the LM algorithm starting from the MMF solution (bottom). The orbital elements corresponding to both solutions are presented in Table 2, and compared to the published solution (Tamuz et al. 2008). 

Open with DEXTER 
HD 156846b is a very eccentric planet (e ≈ 0.85) with a period very close to one year (P ≈ 360 d) discovered by the CORALIE survey (see Tamuz et al. 2008). We used the same data as Tamuz et al. (2008; 53 CORALIE98 measurements, and 11 CORALIE07 measurements). These data are shown in Fig. 4, together with the corresponding periodogram. We note that in Tamuz et al. (2008), CORALIE98 and CORALIE07 were treated as a single instrument because the offset between them is relatively small. In our study we release this constraint and consider both instruments individually (for the offsets).
The proximity of the period to one year induces sampling problems (i.e. a gap in the phasefolded curve). Together with the high eccentricity, the sampling issues prevent an accurate determination of the Fourier coefficients. In particular, the amplitude of the fundamental is underestimated (compared to the harmonics, see the periodogram in Fig. 4 bottom). The highest peak in the periodogram (see Fig. 4 bottom) corresponds to the third harmonics of the planet (period of P/ 4). The determination of the correct fundamental period is thus more difficult than in the case of GJ 3021 (see Sect. 5.1). However, the peak at the fundamental period is still visible, as well as the peaks that correspond to the first and second harmonics. One can thus observe that the four peaks are related (P, P/ 2, P/ 3, P/ 4), and select the correct fundamental period. In addition to this issue, even if the correct period is chosen, the determined Fourier coefficients are not compatible with a Keplerian signal (the ratio between the first harmonics and the fundamental amplitudes is too large). The FF method cannot be used in this case. As an alternative to the FF method, we used the MMF method, which is based on the extrema information (see Appendix C). This method provides crude estimates of the orbital parameters that we use as an initial step for the LM algorithm. Table 2 shows the orbital elements provided by the MMF and LM algorithms compared with the published ones (from Tamuz et al. 2008). Figure 5 shows the phasefolded data superimposed with both solutions (MMF and LM). In Table 2 and Fig. 5, we observe that the orbital parameters provided by the MMF algorithm are a good first approximation. The main issue is in the determination of instrumental offsets (see Fig. 5). Finally, the LM solution is consistent with the published solution (see Table 2).
The case of HD 156846 shows that sampling issues (associated with high eccentricity) may prevent the use of the FF algorithm. In such cases, the MMF solution seems to provide acceptable estimates for the orbital parameters. These estimates provide a very good initialization step for the LM algorithm.
5.3. HD 192310
Fig. 6 Radialvelocity data of HD 192310. The data are taken from Pepe et al. (2011; 139 HARPS measurements). 

Open with DEXTER 
Fig. 7 Periodogram of the radialvelocities of HD 192310 (A)), of the residuals of the oneKeplerian model (LM1 see Table 3, B)), and of the residuals of the twoKeplerians model (LM2 see Table 3, C)). The orbital elements corresponding to the successive solutions are presented in Table 3, and compared to the published solution (Pepe et al. 2011). 

Open with DEXTER 
Orbital parameters of HD 192310b, c obtained using our FF algorithm, using the LM refinement algorithm (starting from the FF solution), and compared with the published solution (PUB, Pepe et al. 2011, on the same 139 HARPS measurements).
Fig. 8 Phase folded radialvelocity curves of HD 192310b (top), and c (bottom), superimposed with our final twoKeplerians model (LM2). See Table 3 for the corresponding orbital elements. 

Open with DEXTER 
HD 192310 hosts two planets discovered with HARPS (see Pepe et al. 2011). We use the 139 HARPS measurements provided by Pepe et al. (2011) for our analysis. These data are shown in Fig. 6 and the corresponding periodogram is shown in Fig. 7 (panel A). The semiamplitudes of the radial velocity signal corresponding to both planets are small (about 3 m s^{1}, see Pepe et al. 2011). Thus, the radial velocity signal has a low S/N, which is challenging for our algorithm. We applied the FF method on the inner planet (Planet b) first, which provided us with a first solution FF1. We then refined this solution with the LM algorithm (LM1). Figure 7 (panel B) shows the periodogram of the residuals of this oneplanet solution (LM1). We then applied the FF algorithm again on the residuals to obtain the orbital parameters of Planet c (FF2). Finally, we refined the global solution (Planets b and c) using the LM method (LM2). Figure 7 (panel C) shows the periodogram of the residuals of this twoplanet solution (LM2). The four successive solutions (FF1, LM1, FF2, and LM2) are shown in Table 3 and compared with the published solution (PUB, Pepe et al. 2011). The phasefolded radial velocity data superimposed with the final model (LM2) are shown in Fig. 8.
The final LM2 solution agrees very well with the published solution (see Table 3). The intermediate solutions (FF1, LM1, and FF2) exhibit larger differences but remain close to the LM2 solution. The most significant differences is observed in the eccentricities. The eccentricity of Planet b is slightly overestimated in FF1 and LM1, while the eccentricity of Planet c is slightly underestimated in FF2 (see Table 3). The complete fit (LM2) is necessary to solve these issue and refine the final model.
The case of HD 192310 shows that the FF method can be used even for multiplanetary systems with low S/N, and still provides good estimates for the orbital parameters.
5.4. HD 147018
Fig. 9 Same as Fig. 6 but for HD 147018. The data are taken from Ségransan et al. (2010) and consist of 6 CORALIE98 measurements (red) and 95 CORALIE07 measurements (blue). 

Open with DEXTER 
Fig. 10 Same as Fig. 7 but for HD 147018. The orbital elements corresponding to the successive solutions are presented in Table 4, and compared to the published solution (Ségransan et al. 2010). 

Open with DEXTER 
Fig. 11 Same as Fig. 8 but for HD 147018. See Table 4 for the corresponding orbital elements. 

Open with DEXTER 
HD 147018 hosts two planets discovered by the CORALIE survey (see Ségransan et al. 2010). We analyzed the same data as Ségransan et al. (2010; 6 CORALIE98 and 95 CORALIE07 measurements), which are shown in Fig. 9. The periodogram is shown in Fig. 10 (panel A). Both planets have equivalent semiamplitudes (≈140 m s^{1}, see Ségransan et al. 2010). We applied the same steps as for HD 192310 (see Sect. 5.3). The four successive solutions (FF1, LM1, FF2, and LM2) are shown in Table 4. The periodogram of the residuals of solution LM1 and LM2 are shown in Fig. 10 (panels B and C). The phasefolded radial velocity data superimposed with the LM2 model is shown in Fig. 11.
The final LM2 solution is consistent with the published solution (see Table 4). The intermediate solutions (FF1, LM1, and FF2) exhibit significant differences (e.g. on the outer period or the angles ω, see Table 4). This is probably due to the iterative method we used: the planets are fitted one after the other by looking at the residuals of the previous model. The case of HD 147018 is representative of the worst situation for such an iterative method. Indeed, both planets have very similar semiamplitudes so they both contribute to the same amount in the radial velocity signal. When fitting the data for the first planet (FF1, LM1), the second planet signal acts as a noise (which, in the case of HD 147018, has the same amplitude as the first planet signal).
The case of HD 147018 shows that the FF method alone cannot be used to reliably fit multiplanetary systems with equivalent semiamplitudes. However it is still a good first approximation to use as an initial step for an overall (all planets at a time) numerical fit (e.g. LM).
6. Discussion
In this article we describe an analytical method to retrieve the orbital parameters of a planet from the periodogram of a radial velocity signal. The method can be outlined as follows:

1.
Find the period (P) of the planet from the periodogram.

2.
Compute the Fourier coefficients of the fundamental and the first harmonics (V_{1}, V_{2}, see Eqs. (57), (58)).

3.
Compute the orbital parameters (e, M_{0}, K, ω) from these Fourier coefficients (see Eqs. (24), (25), (33), (34)).

4.
Possibly refine the parameters with a numerical method.

5.
Compute the residuals and reiterate to find other planets.
This analytical method is very efficient for retrieving parameters compared to numerical methods. The determination of the orbital parameters from the Fourier coefficients V_{1}, V_{2} is very accurate, even at high eccentricities (see Fig. 1). The main limitation of the technique comes from the estimation of the Fourier coefficients. The accuracy of these coefficients (and thus of the final parameters) depends greatly on the sampling of the signal. In some illsampled cases, the inaccuracy of the Fourier coefficients prevents the use of our method (see Sect. 5.2). In such cases, an alternative method that is based on the extrema information can be used to obtain estimates for the orbital parameters (see Appendix C).
Our analytical method is complementary with numerical algorithms since it provides very good initialization values for them. We have illustrated this complementarity by using the LevenbergMarquardt numerical method in combination with our analytical algorithm (see Sect. 5). In particular, for multiplanetary cases, a global numerical fit can significantly improve the solution, especially when the semiamplitudes of the different planets are comparable (see Sect. 5.4). We emphasize the fact that the aim of this method is not to replace more complex numerical algorithms, but to complement them. Our method can be used in combination with any numerical method (LevenbergMarquardt, Markov chain Monte Carlo, genetic algorithms, etc.), and can significantly reduce the parameter space that needs to be explored, and improve the convergence speed. This is especially relevant for multiplanetary systems for which the parameter space has many dimensions and the convergence of numerical algorithm can be slow.
Our methods are all implemented in the DACE platform (Buchschacher et al. 2015), available at https://dace.unige.ch (under the project Observations).
The DACE platform is available at http://dace.unige.ch. The algorithms described in this article are part of the project “Observations” of DACE.
Acknowledgments
We thank the anonymous referee for his/her useful comments. This work has, in part, been carried out within the framework of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge financial support from the SNSF. This publication makes use of DACE, a Data Analysis Center for Exoplanets, a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, based at the University of Geneva (CH).
References
 Buchschacher, N., Ségransan, D., Udry, S., & Díaz, R. 2015, in Astronomical Data Analysis Software an Systems XXIV (ADASS XXIV), eds. A. R. Taylor, & E. Rosolowsky, ASP Conf. Ser., 495, 7 [Google Scholar]
 Correia, A. C. M. 2008, in Precision Spectroscopy in Astrophysics, eds. N. C. Santos, L. Pasquini, A. C. M. Correia, & M. Romaniello, 207 [Google Scholar]
 Correia, A. C. M., Udry, S., Mayor, M., et al. 2005, A&A, 440, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Udry, S., Mayor, M., et al. 2008, A&A, 479, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Udry, S., Mayor, M., et al. 2009, A&A, 496, 521 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Couetdic, J., Laskar, J., et al. 2010, A&A, 511, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Konacki, M., Maciejewski, A. J., & Wolszczan, A. 2002, ApJ, 567, 566 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 2005, Celest. Mech. Dyn. Astron., 91, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Naef, D., Mayor, M., Pepe, F., et al. 2001, A&A, 375, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pepe, F., Lovis, C., Ségransan, D., et al. 2011, A&A, 534, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ségransan, D., Udry, S., Mayor, M., et al. 2010, A&A, 511, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tamuz, O., Ségransan, D., Udry, S., et al. 2008, A&A, 480, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Computation of Hansen coefficients
In this Appendix, we show how to compute Hansen coefficients, both analytically (power series of eccentricity) and numerically. There is a very broad literature on this topic (e.g. Laskar 2005; Laskar & Boué 2010). Here we simply aim at illustrating how to compute a subfamily of Hansen coefficients: (k ∈Z). In this restricted case, the definition of Eq. (7) rewrites (A.1)Using
we obtain (A.5)From this last expression, one only has to expand and e^{ikesinE} in power series of e and integrate over E to obtain the expansion of X_{k}. Equations (10)−(13) provide this expansion for the four coefficients of interest in this study (X_{± 1}, X_{± 2}).
Equation (A.5) is also useful to numerically estimate X_{k}. One simply needs to sample E over [0,2π] and compute the integral (rectangle method) of Eq. (A.5): (A.6)with .
Fig. A.1 Convergence of the numerical estimation of X_{1}(0.99). We observe that the optimum value of N (see Eq. (A.6)) is about 100. For too small N (≲20), the rectangle method does not provide enough precision. For too high N (≳300), the machine errors accumulate and the precision decreases. 

Open with DEXTER 
Figure A.1 illustrates how the numerical estimates of Eq. (A.6) evolves as a function of N (for X_{1}(0,99)). We observe that when N is too small, the rectangle method is not precise enough and when N is too large, machine errors accumulate and the precision is lost. In this study, we use N = 100, which seems to be a good compromise (we tested different Hansen coefficients and different values of e).
Appendix B: Derivatives of Hansen coefficients
The computation of derivatives of Hansen coefficients is very similar to the computation of the coefficients themselves (see Appendix A). From Eq. (A.5) we deduce (B.1)As explained in Appendix A, this expression can be used both to expand the derivatives in power series of e, and to numerically estimate the integral (rectangle method).
Appendix C: Orbital elements from extrema
In some cases, the Fourier coefficients cannot be determined correctly (e.g. owing to a bad sampling). An alternative method, using extrema of the radial velocity curve, can be used to obtain a rough estimate of orbital parameters. We assume that the orbital period (as well as instrumental offsets and drifts) has been correctly determined but not the Fourier coefficients. By looking at the timings and values of the maximum and minimum of the folded radial velocity curve, we again obtain four observables that allow the four parameters to be determined: K, e, ω, M_{0}. We note V_{min}, V_{max}, t_{min}, t_{max} the values and timings (in [0,P]) of the minimum and maximum, assuming that the mean value of V is zero. From Eq. (1), we obtain We have (see Eqs. (C.1), (C.2)) To determine the remaining parameters (esinω, M_{0}), we need to use the information contained in the timing of extrema. We have (C.7)thus These mean anomalies (M) can be expressed as functions of the eccentricity and the true anomalies (v), which are given by Eqs. (C.3), (C.4) Developing Eqs. (C.10), (C.11) at second order in eccentricity, we obtain (C.12)which gives
thus (see Eqs. (C.8), (C.9)) (C.16)This provides an estimate for esinω, and thus for e and ω (since we already have an estimate for ecosω see Eq. (C.6)).
Finally, using the estimates for e and ω in Eqs. (C.10), (C.11), we obtain the values of M_{min}, M_{max}. The remaining parameter M_{0} is straightforwardly derived from these values (C.17)The estimates given in Eqs. (C.5), (C.6), (C.16), (C.17) rely on the determination of the period P and the values and timings of the extrema (V_{min}, V_{max}, t_{min}, and t_{max}). We use the period estimate given by the periodogram method (see Sect. 4). The extrema determination can be affected by offsets, drifts, and outliers. We thus substract from the data, the offsets and drifts determined by the linear fit used for the Fourier method (see Sect. 4). Moreover, to reduce the impact of outliers, we compute the minimum and maximum (values and timings) as a weighted mean over the lowest and highest N points (we chose N = 2 for HD 156846b, see Sect. 5.2). This method provides very crude estimates for the parameters and is mainly useful when the Fourier method fails.
All Tables
Orbital parameters of GJ 3021b obtained using our FF algorithm, using the LM refinement algorithm (starting from the FF solution), and compared with the published solution (PUB, Naef et al. 2001, on the same 61 CORALIE98 measurements).
Orbital parameters of HD 192310b, c obtained using our FF algorithm, using the LM refinement algorithm (starting from the FF solution), and compared with the published solution (PUB, Pepe et al. 2011, on the same 139 HARPS measurements).
All Figures
Fig. 1 Comparison between our analytical estimates (see Eqs. (24), (25), (33), (34)) and the expected values of the eccentricity (e), the mean longitude at the reference time (λ = M_{0} + ω), the semiamplitude K, and the argument of periastron (ω) (from top to bottom), as a function of ω and for e = 0.5, 0.8, 0.9, 0.95 (from left to right). We show the errors obtained on the parameters using the approximations C = 1/4, and (see Sect. 3.2). For the approximation , we also show the errors after one and two iterations of the NewtonRaphson correcting algorithm (see Sect. 3.4). 

Open with DEXTER  
In the text 
Fig. 2 Radialvelocity data of GJ 3021 (top) and corresponding periodogram (bottom). The data are taken from Naef et al. (2001; 61 CORALIE98 measurements). 

Open with DEXTER  
In the text 
Fig. 3 Radialvelocity data of GJ 3021 superimposed with the Keplerian model obtained with the FF algorithm (top), and the model obtained with the LM algorithm starting from the FF solution (bottom). The orbital elements corresponding to both solutions are presented in Table 1, and compared to the published solution (Naef et al. 2001). 

Open with DEXTER  
In the text 
Fig. 4 Same as Fig. 2 but for HD 156846. The data are taken from Tamuz et al. (2008). They consist of 53 CORALIE98 measurements (red), and 11 CORALIE07 measurements (blue). 

Open with DEXTER  
In the text 
Fig. 5 Phasefolded radialvelocity curve of HD 156846b superimposed with the Keplerian model that was obtained with the MMF algorithm (top), and the model obtained with the LM algorithm starting from the MMF solution (bottom). The orbital elements corresponding to both solutions are presented in Table 2, and compared to the published solution (Tamuz et al. 2008). 

Open with DEXTER  
In the text 
Fig. 6 Radialvelocity data of HD 192310. The data are taken from Pepe et al. (2011; 139 HARPS measurements). 

Open with DEXTER  
In the text 
Fig. 7 Periodogram of the radialvelocities of HD 192310 (A)), of the residuals of the oneKeplerian model (LM1 see Table 3, B)), and of the residuals of the twoKeplerians model (LM2 see Table 3, C)). The orbital elements corresponding to the successive solutions are presented in Table 3, and compared to the published solution (Pepe et al. 2011). 

Open with DEXTER  
In the text 
Fig. 8 Phase folded radialvelocity curves of HD 192310b (top), and c (bottom), superimposed with our final twoKeplerians model (LM2). See Table 3 for the corresponding orbital elements. 

Open with DEXTER  
In the text 
Fig. 9 Same as Fig. 6 but for HD 147018. The data are taken from Ségransan et al. (2010) and consist of 6 CORALIE98 measurements (red) and 95 CORALIE07 measurements (blue). 

Open with DEXTER  
In the text 
Fig. 10 Same as Fig. 7 but for HD 147018. The orbital elements corresponding to the successive solutions are presented in Table 4, and compared to the published solution (Ségransan et al. 2010). 

Open with DEXTER  
In the text 
Fig. 11 Same as Fig. 8 but for HD 147018. See Table 4 for the corresponding orbital elements. 

Open with DEXTER  
In the text 
Fig. A.1 Convergence of the numerical estimation of X_{1}(0.99). We observe that the optimum value of N (see Eq. (A.6)) is about 100. For too small N (≲20), the rectangle method does not provide enough precision. For too high N (≳300), the machine errors accumulate and the precision decreases. 

Open with DEXTER  
In the text 