© 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 least-squares minimization methods (Levenberg-Marquardt, 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 (non-linear fit can have several local χ₂ 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 least-squares 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. ill-sampled, very eccentric planets). These algorithms are implemented in the DACE web platform (see Buchschacher et al. 2015) ¹.

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 $\mathit{V}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{K}\left(\cos \mathrm{(}\mathit{\nu }\mathrm{+}\mathit{\omega }\mathrm{)}\mathrm{+}\mathit{e}\cos \mathrm{\left(}\mathit{\omega }\mathrm{\right)}\right)\mathit{,}$(1)with $\mathit{K}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{p}}}{{\mathit{m}}_{\mathrm{s}}\mathrm{+}{\mathit{m}}_{\mathrm{p}}}\frac{\mathrm{2}\mathit{\pi a}\sin \mathit{i}}{\mathit{P}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}}\mathit{,}$(2)and a is the semi-major 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 mean-motion of the planet: $\mathit{n}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathit{P}}\mathrm{=}\sqrt{\frac{\mathit{\mu }}{{\mathit{a}}^{\mathrm{3}}}}\mathit{,}$(3)with $\mathit{\mu }\mathrm{=}\mathrm{𝒢}\mathrm{(}{\mathit{m}}_{\mathrm{s}}\mathrm{+}{\mathit{m}}_{\mathrm{p}}\mathrm{)}\mathit{,}$(4)and is the gravitational constant. The radial velocity signal is P-periodic and can be decomposed in discrete Fourier series: $\mathit{V}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}\mathrm{\in }Z}{\mathit{V}}_{\mathit{k}}{\mathrm{e}}^{\mathrm{i}\mathit{knt}}\mathit{,}$(5)with $\begin{array}{ccc}{\mathit{V}}_{\mathit{k}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{P}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{P}}\mathit{V}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{knt}}\mathrm{d}\mathit{t}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathit{V}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k}\mathrm{(}\mathit{M}\mathrm{-}{\mathit{M}}_{\mathrm{0}}\mathrm{)}}\mathrm{d}\mathit{M,}\end{array}$(6)where M is the mean anomaly and M₀ is the mean anomaly at the reference time t = 0. The coefficients V_k are complex numbers. Since V(t) is real, we have ${\mathit{V}}_{\mathrm{-}\mathit{k}}\mathrm{=}{\overline{\mathit{V}}}_{\mathit{k}}$ (where ${\overline{\mathit{V}}}_{\mathit{k}}$ is the complex conjugate of V_k). We denote by ${\mathit{X}}_{\mathit{k}}^{\mathit{l,q}}$ the Hansen coefficients ${\mathit{X}}_{\mathit{k}}^{\mathit{l,q}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\left(\frac{\mathit{r}}{\mathit{a}}\right)}^{\mathit{l}}{\mathrm{e}}^{\mathrm{i}\mathit{q\nu }}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{kM}}\mathrm{d}\mathit{M}\mathit{.}$(7)These coefficients are real numbers that only depend on the eccentricity (e), and we have ${\mathit{X}}_{\mathrm{-}\mathit{k}}^{\mathit{l,}\mathrm{-}\mathit{q}}\mathrm{=}{\mathit{X}}_{\mathit{k}}^{\mathit{l,q}}$. The Fourier expansion of the radial velocity can be rewritten using Hansen coefficients $\begin{array}{ccc}{\mathit{V}}_{\mathrm{0}}& \mathrm{=}& \mathrm{0}\\ {\mathit{V}}_{\mathit{k}}& \mathrm{=}& \frac{\mathit{K}{\mathrm{e}}^{\mathrm{i}\mathit{k}{\mathit{M}}_{\mathrm{0}}}}{\mathrm{2}}\mathrm{(}{\mathit{X}}_{\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }}\mathrm{+}{\mathit{X}}_{\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{-}\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}\mathrm{)}\\ & \mathrm{=}& \frac{\mathit{K}{\mathrm{e}}^{\mathrm{i}\mathit{k}{\mathit{M}}_{\mathrm{0}}}}{\mathrm{2}}\mathrm{(}{\mathit{X}}_{\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }}\mathrm{+}{\mathit{X}}_{\mathrm{-}\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}\mathrm{)} \mathrm{(}\mathit{k}\mathrm{\in }{Z}^{\mathrm{\ast }}\mathrm{)}\mathit{.}\end{array}$We observe that the radial velocity only depends on Hansen coefficients of the form ${\mathit{X}}_{\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{1}}$ (k ∈ Z^∗). In the following, we drop the exponents and use the notation ${\mathit{X}}_{\mathit{k}}\mathrm{\equiv }{\mathit{X}}_{\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{1}}$.

We note that V₀ = 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 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₁ (fundamental) and V₂ (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₀.

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) $\begin{array}{ccc}{\mathit{X}}_{\mathrm{1}}& \mathrm{=}& \mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{4}}\mathrm{)}\\ {\mathit{X}}_{\mathrm{2}}& \mathrm{=}& \mathit{e}\mathrm{-}\frac{\mathrm{5}}{\mathrm{4}}{\mathit{e}}^{\mathrm{3}}\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{5}}\mathrm{)}\\ {\mathit{X}}_{-1}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{8}}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{48}}{\mathit{e}}^{\mathrm{4}}\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{6}}\mathrm{)}\\ {\mathit{X}}_{-2}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{12}}{\mathit{e}}^{\mathrm{3}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{48}}{\mathit{e}}^{\mathrm{5}}\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{7}}\mathrm{)}\mathit{.}\end{array}$

3.2. Eccentricity and phase

Since we assume that the complex coefficients V₁, V₂ are known, the ratio ρ ≡ V₂/V₁ is also known. From Eq. (9), we have $\mathit{\rho }\mathrm{=}\frac{{\mathit{V}}_{\mathrm{2}}}{{\mathit{V}}_{\mathrm{1}}}\mathrm{=}{\mathrm{e}}^{\mathrm{i}{\mathit{M}}_{\mathrm{0}}}\frac{{\mathit{X}}_{\mathrm{2}}}{{\mathit{X}}_{\mathrm{1}}}\left(\frac{\mathrm{1}\mathrm{+}{\mathit{X}}_{-2}\mathit{/}{\mathit{X}}_{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathrm{i}\mathit{\omega }}}{\mathrm{1}\mathrm{+}{\mathit{X}}_{-1}\mathit{/}{\mathit{X}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathrm{i}\mathit{\omega }}}\right)\mathrm{·}$(14)From Eqs. (10)−(13) we obtain $\begin{array}{ccc}\frac{{\mathit{X}}_{\mathrm{2}}}{{\mathit{X}}_{\mathrm{1}}}& \mathrm{=}& \mathit{e}\mathrm{-}\frac{{\mathit{e}}^{\mathrm{3}}}{\mathrm{4}}\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{5}}\mathrm{)}\\ \frac{{\mathit{X}}_{-1}}{{\mathit{X}}_{\mathrm{1}}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{e}}^{\mathrm{2}}}{\mathrm{8}}\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{4}}\mathrm{)}\\ \frac{{\mathit{X}}_{-2}}{{\mathit{X}}_{\mathrm{2}}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{e}}^{\mathrm{2}}}{\mathrm{12}}\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{4}}\mathrm{)}\mathit{.}\end{array}$We thus have $\frac{\mathrm{1}\mathrm{+}{\mathit{X}}_{-2}\mathit{/}{\mathit{X}}_{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathrm{i}\mathit{\omega }}}{\mathrm{1}\mathrm{+}{\mathit{X}}_{-1}\mathit{/}{\mathit{X}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathrm{i}\mathit{\omega }}}\mathrm{=}\mathrm{1}\mathrm{+}\frac{{\mathit{e}}^{\mathrm{2}}}{\mathrm{24}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathrm{i}\mathit{\omega }}\mathrm{+}\mathrm{O}\mathrm{\left(}{\mathit{e}}^{\mathrm{4}}\mathrm{\right)}\mathit{,}$(18)and $\mathit{\rho }\mathrm{=}{\mathrm{e}}^{\mathrm{i}{\mathit{M}}_{\mathrm{0}}}\left(\mathit{e}\mathrm{-}\frac{{\mathit{e}}^{\mathrm{3}}}{\mathrm{4}}\mathrm{+}\frac{{\mathit{e}}^{\mathrm{3}}}{\mathrm{24}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathrm{i}\mathit{\omega }}\mathrm{+}\mathrm{O}\mathrm{\left(}{\mathit{e}}^{\mathrm{5}}\mathrm{\right)}\right)\mathit{.}$(19)We introduce the complex coefficient $\mathit{C}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\left(\mathrm{1}\mathrm{-}\frac{{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathrm{i}\mathit{\omega }}}{\mathrm{6}}\right)\mathit{,}$(20)and use, in the following, the approximation $\mathit{\rho }\mathrm{\approx }{\mathrm{e}}^{\mathrm{i}{\mathit{M}}_{\mathrm{0}}}\mathrm{(}\mathit{e}\mathrm{-}\mathit{C}\mathrm{(}\mathit{\omega }\mathrm{\left)}{\mathit{e}}^{\mathrm{3}}\mathrm{\right)}\mathit{,}$(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₀. Using approximation (21), we only have to solve a third order polynomial equation to obtain the eccentricity $\mathrm{\left|}\mathit{\rho }\mathrm{\right|}\mathrm{\approx }\mathit{e}\mathrm{-}\mathrm{\Re }\mathrm{\left(}\mathit{C}\mathrm{\right)}{\mathit{e}}^{\mathrm{3}}\mathit{,}$(22)where ℜ(C) denotes the real part of C, $\mathrm{\Re }\mathrm{\left(}\mathit{C}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\left(\mathrm{1}\mathrm{-}\frac{\cos \mathrm{\left(}\mathrm{2}\mathit{\omega }\mathrm{\right)}}{\mathrm{6}}\right) \mathrm{\in }\left[\frac{\mathrm{5}}{\mathrm{24}}\mathit{,}\frac{\mathrm{7}}{\mathrm{24}}\right]\mathrm{·}$(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 $\mathit{ê}\mathrm{=}\frac{\mathrm{2}}{\sqrt{\mathrm{3}\mathrm{\Re }\mathrm{\left(}\mathit{C}\mathrm{\right)}}}\cos \left(\frac{\mathit{\pi }\mathrm{+}\arccos \left(\frac{\mathrm{3}\sqrt{\mathrm{3}\mathrm{\Re }\mathrm{\left(}\mathit{C}\mathrm{\right)}}}{\mathrm{2}}\mathrm{\left|}\mathit{\rho }\mathrm{\right|}\right)}{\mathrm{3}}\right)\mathit{,}$(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₀ is estimated by $\mathit{M̂}\mathrm{0}\mathrm{=}\arg \left(\frac{\mathit{\rho }}{\mathit{e}\mathrm{-}\mathit{C}{\mathit{e}}^{\mathrm{3}}}\right)\mathrm{·}$(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 $\mathit{\eta }\mathrm{=}\frac{{\mathit{V}}_{\mathrm{2}}}{{\mathit{V}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{=}\frac{\mathrm{2}}{\mathit{K}}\frac{{\mathit{X}}_{\mathrm{2}}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }}\mathrm{+}{\mathit{X}}_{-2}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}}{{\left({\mathit{X}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }}\mathrm{+}{\mathit{X}}_{-1}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}\right)}^{\mathrm{2}}}\mathit{,}$(26)which enables us to eliminate M₀. We know that ${\mathit{X}}_{\mathrm{-}\mathit{k}}\mathit{/}{\mathit{X}}_{\mathit{k}}\mathrm{=}\mathrm{O}{}^{\mathrm{(}}\mathit{e}^{\mathrm{2}}{}^{\mathrm{)}}$ (see Eqs. (16), (17)), thus, at leading order in eccentricity $\mathit{\eta }\mathrm{\approx }\frac{\mathrm{2}{\mathit{X}}_{\mathrm{2}}}{\mathit{K}{\mathit{X}}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}\mathit{,}$(27)and the phase of η provides a very crude estimate of ω$\mathit{\omega ̂}\mathrm{=}\mathrm{-}\arg \mathrm{\left(}\mathit{\eta }\mathrm{\right)}\mathit{,}$(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 ${\mathit{V}}_{\mathrm{1}}\mathrm{=}\frac{\mathit{K}{\mathrm{e}}^{\mathrm{i}{\mathit{M}}_{\mathrm{0}}}}{\mathrm{2}}\mathrm{\left(}{\mathit{X}}_{\mathrm{1}}\mathrm{\right(}\mathit{e}\mathrm{)}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }}\mathrm{+}{\mathit{X}}_{-1}\mathrm{(}\mathit{e}\mathrm{\left)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}\mathrm{\right)}\mathit{,}$(29)and rewrite it has a function of Kcosω, Ksinω$\begin{array}{ccc}{\mathit{V}}_{\mathrm{1}}\mathrm{=}\frac{{\mathrm{e}}^{\mathrm{i}{\mathit{M}}_{\mathrm{0}}}}{\mathrm{2}}& & (\mathit{K}\cos \mathit{\omega }\left({\mathit{X}}_{\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{+}{\mathit{X}}_{-1}\mathrm{\left(}\mathit{e}\mathrm{\right)}\right)\\ & & \mathrm{+}\mathrm{i}\mathit{K}\sin \mathit{\omega }\left({\mathit{X}}_{\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{-}{\mathit{X}}_{-1}\mathrm{\left(}\mathit{e}\mathrm{\right)}\right))\mathit{.}\end{array}$(30)Therefore $\begin{array}{ccc}\mathit{K}\cos \mathit{\omega }& \mathrm{=}& \frac{\mathrm{2}}{{\mathit{X}}_{\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{+}{\mathit{X}}_{-1}\mathrm{\left(}\mathit{e}\mathrm{\right)}}\mathrm{\Re }\mathrm{(}{\mathit{V}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{M}}_{\mathrm{0}}}\mathrm{)}\\ \mathit{K}\sin \mathit{\omega }& \mathrm{=}& \frac{\mathrm{2}}{{\mathit{X}}_{\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{-}{\mathit{X}}_{-1}\mathrm{\left(}\mathit{e}\mathrm{\right)}}\mathrm{\Im }\mathrm{(}{\mathit{V}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{M}}_{\mathrm{0}}}\mathrm{)}\mathit{.}\end{array}$We obtain estimates for K and ω by introducing the estimates of e and M₀ in Eqs. (31), (32): $\begin{array}{ccc}\mathit{K̂}\cos \mathit{\omega ̂}& \mathrm{=}& \frac{\mathrm{2}}{{\mathit{X}}_{\mathrm{1}}\mathrm{\left(}\mathit{ê}\mathrm{\right)}\mathrm{+}{\mathit{X}}_{-1}\mathrm{\left(}\mathit{ê}\mathrm{\right)}}\mathrm{\Re }\mathrm{(}{\mathit{V}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{M̂}\mathrm{0}}\mathrm{)}\\ \mathit{K̂}\sin \mathit{\omega ̂}& \mathrm{=}& \frac{\mathrm{2}}{{\mathit{X}}_{\mathrm{1}}\mathrm{\left(}\mathit{ê}\mathrm{\right)}\mathrm{-}{\mathit{X}}_{-1}\mathrm{\left(}\mathit{ê}\mathrm{\right)}}\mathrm{\Im }\mathrm{(}{\mathit{V}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{M̂}\mathrm{0}}\mathrm{)}\mathit{.}\end{array}$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₀ and ω are meaningless quantities, while the mean longitude at reference time λ₀ = M₀ + ω is the only important angle. Our analytical method provides an accurate estimate for λ₀ even in the low eccentricity case. Indeed, at low eccentricity, the estimate (see Eq. (25)) may be incorrect because the coefficient V₂ is small and more sensitive to noise. However, for small e, we have X₁(e) ≈ 1 and X_-1(e) ≪ X₁(e), and thus ${\mathit{V}}_{\mathrm{1}}\mathrm{\approx }\frac{\mathit{K}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{i}{\mathit{\lambda }}_{\mathrm{0}}}\mathit{.}$(35)The phase of V₁ directly provides an estimate of λ₀. Still assuming small eccentricity, we have (see Eqs. (33), (34)) $\mathit{K̂}{\mathrm{e}}^{\mathrm{i}\mathit{\omega ̂}}\mathrm{\approx }\mathrm{2}{\mathit{V}}_{\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{M̂}\mathrm{0}}\mathit{.}$(36)Therefore, is also incorrect, but the estimate is not affected: $\mathit{K̂}{\mathrm{e}}^{\mathrm{i}\mathit{\lambda ̂}\mathrm{0}}\mathrm{\approx }\mathit{K}{\mathrm{e}}^{\mathrm{i}{\mathit{\lambda }}_{\mathrm{0}}}\mathit{.}$(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 Newton-Raphson method. We denote by x the vector of parameters and y the vector of Fourier coefficients: $\mathit{x}\mathrm{=}\left(\begin{array}{c}\\ \mathit{K}\\ \mathit{e}\\ \mathit{\omega }\\ {\mathit{M}}_{\mathrm{0}}\end{array}\right)\mathit{,} \mathit{y}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{\Re }\mathrm{\left(}{\mathit{V}}_{\mathrm{1}}\mathrm{\right)}\\ \mathrm{\Im }\mathrm{\left(}{\mathit{V}}_{\mathrm{1}}\mathrm{\right)}\\ \mathrm{\Re }\mathrm{\left(}{\mathit{V}}_{\mathrm{2}}\mathrm{\right)}\\ \mathrm{\Im }\mathrm{\left(}{\mathit{V}}_{\mathrm{2}}\mathrm{\right)}\end{array}\right)\mathit{.}$(38)We denote by J the Jacobian matrix of our system $\mathit{J}\mathrm{=}\left(\begin{array}{c}\\ \frac{\mathit{\partial }{\mathit{y}}_{\mathrm{1}}}{\mathit{\partial }{\mathit{x}}_{\mathrm{1}}}& \mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}& \frac{\mathit{\partial }{\mathit{y}}_{\mathrm{1}}}{\mathit{\partial }{\mathit{x}}_{\mathrm{4}}}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& {\begin{array}{c}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \frac{\mathit{\partial }{\mathit{y}}_{\mathrm{4}}}{\mathit{\partial }{\mathit{x}}_{\mathrm{1}}}& \mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}& \frac{\mathit{\partial }{\mathit{y}}_{\mathrm{4}}}{\mathit{\partial }{\mathit{x}}_{\mathrm{4}}}\end{array}\right)\mathit{.}$(39)The Newton-Raphson method provides corrections for the estimates of the parameters x$\mathit{\delta }\mathit{x̂}\mathrm{=}{\mathit{J}}^{-1}\mathrm{(}\mathit{y}\mathrm{-}\mathit{ŷ}\mathrm{)}\mathit{.}$(40)From our estimates of the parameters, we can compute estimates for V₁, V₂ (see Eq. (9)), and thus determine the ŷ vector and compare it to the known Fourier coefficients (vector y). The estimation of V₁, V₂ 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₁, V₂. From Eq. (9) we deduce $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{V}}_{\mathit{k}}}{\mathit{\partial K}}& \mathrm{=}& \frac{{\mathit{V}}_{\mathit{k}}}{\mathit{K}}\mathrm{=}\frac{{\mathrm{e}}^{\mathrm{i}\mathit{k}{\mathit{M}}_{\mathrm{0}}}}{\mathrm{2}}\mathrm{(}{\mathit{X}}_{\mathit{k}}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }}\mathrm{+}{\mathit{X}}_{\mathrm{-}\mathit{k}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}\mathrm{)}\mathit{,}\\ \frac{\mathit{\partial }{\mathit{V}}_{\mathit{k}}}{\mathit{\partial e}}& \mathrm{=}& \frac{\mathit{K}{\mathrm{e}}^{\mathrm{i}\mathit{k}{\mathit{M}}_{\mathrm{0}}}}{\mathrm{2}}\left(\frac{\mathrm{d}{\mathit{X}}_{\mathit{k}}}{\mathrm{d}\mathit{e}}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }}\mathrm{+}\frac{\mathrm{d}{\mathit{X}}_{\mathrm{-}\mathit{k}}}{\mathrm{d}\mathit{e}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}\right)\mathit{,}\\ \frac{\mathit{\partial }{\mathit{V}}_{\mathit{k}}}{\mathit{\partial \omega }}& \mathrm{=}& \mathrm{i}\frac{\mathit{K}{\mathrm{e}}^{\mathrm{i}\mathit{k}{\mathit{M}}_{\mathrm{0}}}}{\mathrm{2}}\mathrm{(}{\mathit{X}}_{\mathit{k}}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }}\mathrm{-}{\mathit{X}}_{\mathrm{-}\mathit{k}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}\mathrm{)}\mathit{,}\\ \frac{\mathit{\partial }{\mathit{V}}_{\mathit{k}}}{\mathit{\partial }{\mathit{M}}_{\mathrm{0}}}& \mathrm{=}& \mathrm{i}\mathit{k}{\mathit{V}}_{\mathit{k}}\mathrm{=}\mathrm{i}\mathit{k}\frac{\mathit{K}{\mathrm{e}}^{\mathrm{i}\mathit{k}{\mathit{M}}_{\mathrm{0}}}}{\mathrm{2}}\mathrm{(}{\mathit{X}}_{\mathit{k}}{\mathrm{e}}^{\mathrm{i}\mathit{\omega }}\mathrm{+}{\mathit{X}}_{\mathrm{-}\mathit{k}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\omega }}\mathrm{)}\mathit{.}\end{array}$We thus need to determine estimates of the derivatives of Hansen coefficients $\frac{\mathrm{d}{\mathit{X}}_{\mathit{k}}}{\mathrm{d}\mathit{e}}$ (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₁ and V₂ 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, 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₀ 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 Newton-Raphson 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 Newton-Raphson 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 Newton-Raphson correction in our algorithm. These encouraging results are obtained assuming that the Fourier coefficients V₁, V₂ have been determined very accurately. In real cases the errors might be much greater if the estimates of V₁, V₂ 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 (λ = M0 + ω), the semi-amplitude 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 Newton-Raphson correcting algorithm (see Sect. 3.4).

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₁, V₂ (see Eq. (6)) from the radial velocity signal. The classical method for this is to use a least-square spectral analysis (such as the Lomb-Scargle 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 ${\mathit{v}}_{\mathit{i}}\mathrm{=}\mathit{V}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathrm{+}{\mathit{ϵ}}_{\mathit{i}}\mathit{,}$(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 $\mathit{V}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{f}}}{\mathit{C}}_{\mathit{k}}{\mathit{f}}_{\mathit{k}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}$(46)To find the coefficients C_k, we perform a weighted least-squares regression. We introduce the matrices $\mathit{v}\mathrm{=}\left(\begin{array}{c}\\ {\mathit{v}}_{\mathrm{1}}\mathit{/}{\mathit{\sigma }}_{\mathrm{1}}\\ {\mathit{v}}_{\mathrm{2}}\mathit{/}{\mathit{\sigma }}_{\mathrm{2}}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{v}}_{{\mathit{N}}_{\mathrm{m}}}\mathit{/}{\mathit{\sigma }}_{{\mathit{N}}_{\mathrm{m}}}\end{array}\right)\mathit{,} \mathit{C}\mathrm{=}\left(\begin{array}{c}\\ {\mathit{C}}_{\mathrm{1}}\\ {\mathit{C}}_{\mathrm{2}}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{C}}_{{\mathit{N}}_{\mathit{f}}}\end{array}\right)\mathit{,}$(47)$\mathit{f}\mathrm{=}\left(\begin{array}{c}\\ {\mathit{f}}_{\mathrm{1}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{\mathrm{1}}& {\mathit{f}}_{\mathrm{2}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{\mathrm{1}}& \mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}& {\mathit{f}}_{{\mathit{N}}_{\mathit{f}}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{\mathrm{1}}\\ {\mathit{f}}_{\mathrm{1}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{\mathrm{2}}& {\mathit{f}}_{\mathrm{2}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{\mathrm{2}}& \mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}& {\mathit{f}}_{{\mathit{N}}_{\mathit{f}}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{\mathrm{2}}\\ \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& {\begin{array}{c}\mathit{.}\end{array}}^{\mathit{.}}\mathit{.}& \begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{f}}_{\mathrm{1}}\mathrm{\left(}{\mathit{t}}_{{\mathit{N}}_{\mathrm{m}}}\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{{\mathit{N}}_{\mathrm{m}}}& {\mathit{f}}_{\mathrm{2}}\mathrm{\left(}{\mathit{t}}_{{\mathit{N}}_{\mathrm{m}}}\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{{\mathit{N}}_{\mathrm{m}}}& \mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}& {\mathit{f}}_{{\mathit{N}}_{\mathit{f}}}\mathrm{\left(}{\mathit{t}}_{{\mathit{N}}_{\mathrm{m}}}\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{{\mathit{N}}_{\mathrm{m}}}\end{array}\right)\mathit{.}$(48)The least-squares regression provides an estimator for C$\mathit{Ĉ}\mathrm{=}{\mathrm{(}{\mathit{f}}^{\mathrm{T}}{\mathit{f}}^{\mathrm{)}}}^{-1}{\mathit{f}}^{\mathrm{T}}\mathit{v}\mathit{.}$(49)The functions f_k that are used in the decomposition basis are usually the sine and cosine of a given period: $\begin{array}{ccc}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \cos \left(\mathrm{2}\mathit{\pi }\frac{\mathit{t}}{\mathit{P}}\right)\\ {\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \sin \left(\mathrm{2}\mathit{\pi }\frac{\mathit{t}}{\mathit{P}}\right)\mathrm{·}\end{array}$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 (χ²(P) which depends on the chosen period P). The best-fitting period is found by minimizing the χ²(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 χ² function (see Zechmeister & Kürster 2009): $\mathrm{𝒫}\mathrm{\left(}\mathit{P}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\chi }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{P}\mathrm{\right)}}{{\mathit{\chi }}_{\mathrm{0}}^{\mathrm{2}}}\mathit{,}$(52)where ${\mathit{\chi }}_{\mathrm{0}}^{\mathrm{2}}$ 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 χ²). We denote by the period that corresponds with this peak. We now need estimates for V₁ and V₂. We obtain these by using the same least-squares algorithm, but we introduce in the function basis f the sines and cosines corresponding to both the fundamental period and the first harmonics : $\begin{array}{ccc}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \cos \left(\mathrm{2}\mathit{\pi }\frac{\mathit{t}}{\mathit{P̂}}\right)\\ {\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \sin \left(\mathrm{2}\mathit{\pi }\frac{\mathit{t}}{\mathit{P̂}}\right)\\ {\mathit{f}}_{\mathrm{3}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \cos \left(\mathrm{4}\mathit{\pi }\frac{\mathit{t}}{\mathit{P̂}}\right)\\ {\mathit{f}}_{\mathrm{4}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \sin \left(\mathrm{4}\mathit{\pi }\frac{\mathit{t}}{\mathit{P̂}}\right)\mathrm{·}\end{array}$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₁ and V₂ are then given by $\begin{array}{ccc}\mathit{V̂}\mathrm{1}& \mathrm{=}& \frac{\mathit{Ĉ}\mathrm{1}\mathrm{-}\mathrm{i}\mathit{Ĉ}\mathrm{2}}{\mathrm{2}}\\ \mathit{V̂}\mathrm{2}& \mathrm{=}& \frac{\mathit{Ĉ}\mathrm{3}\mathrm{-}\mathrm{i}\mathit{Ĉ}\mathrm{4}}{\mathrm{2}}\mathit{,}\end{array}$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 Levenberg-Marquardt 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 Levenberg-Marquardt (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 ill-sampled, very eccentric planet), HD 192310 (two planets with a low S/N), and HD 147018 (two planets with equal semi-amplitudes). For all these analyses, we quadratically added an instrumental systematic error to the radial velocity uncertainty. For CORALIE-98, this instrumental error is set to 6 m s^-1, for CORALIE-07 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₀, 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ω, λ₀ = M₀ + ω is better suited. Indeed, for low eccentricities, M₀ and ω present strong anti-correlations since the periastron direction is ill-defined. Using k, h, λ₀ 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 Radial-velocity data of GJ 3021 (top) and corresponding periodogram (bottom). The data are taken from Naef et al. (2001; 61 CORALIE-98 measurements).

Fig. 3 Radial-velocity 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).

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 CORALIE-98 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 CORALIE-98 measurements (red), and 11 CORALIE-07 measurements (blue).

Fig. 5 Phase-folded radial-velocity 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).

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 CORALIE-98 measurements, and 11 CORALIE-07 measurements). These data are shown in Fig. 4, together with the corresponding periodogram. We note that in Tamuz et al. (2008), CORALIE-98 and CORALIE-07 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 phase-folded 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 phase-folded 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 Radial-velocity data of HD 192310. The data are taken from Pepe et al. (2011; 139 HARPS measurements).

Fig. 7 Periodogram of the radial-velocities of HD 192310 (A)), of the residuals of the one-Keplerian model (LM1 see Table 3, B)), and of the residuals of the two-Keplerians 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).

Fig. 8 Phase folded radial-velocity curves of HD 192310b (top), and c (bottom), superimposed with our final two-Keplerians model (LM2). See Table 3 for the corresponding orbital elements.

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 semi-amplitudes 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 one-planet 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 two-planet 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 phase-folded 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 multi-planetary 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 CORALIE-98 measurements (red) and 95 CORALIE-07 measurements (blue).

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).

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

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 CORALIE-98 and 95 CORALIE-07 measurements), which are shown in Fig. 9. The periodogram is shown in Fig. 10 (panel A). Both planets have equivalent semi-amplitudes (≈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 phase-folded 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 semi-amplitudes 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 multi-planetary systems with equivalent semi-amplitudes. 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₁, V₂, see Eqs. (57), (58)).

• 3.

  Compute the orbital parameters (e, M₀, 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₁, V₂ 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 ill-sampled 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 Levenberg-Marquardt numerical method in combination with our analytical algorithm (see Sect. 5). In particular, for multi-planetary cases, a global numerical fit can significantly improve the solution, especially when the semi-amplitudes 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 (Levenberg-Marquardt, 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 multi-planetary 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).

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

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 sub-family of Hansen coefficients: ${\mathit{X}}_{\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{1}}$ (k ∈Z). In this restricted case, the definition of Eq. (7) rewrites ${\mathit{X}}_{\mathit{k}}\mathrm{=}{\mathit{X}}_{\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{1}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\mathrm{e}}^{\mathrm{i}\mathit{\nu }}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{kM}}\mathrm{d}\mathit{M}\mathit{.}$(A.1)Using

$\begin{array}{ccc}& & \mathit{M}\mathrm{=}\mathit{E}\mathrm{-}\mathit{e}\sin \mathit{E,}\\ & & \mathrm{d}\mathit{M}\mathrm{=}\frac{\mathit{r}}{\mathit{a}}\mathrm{d}\mathit{E,}\\ & & {\mathrm{e}}^{\mathrm{i}\mathit{\nu }}\mathrm{=}\frac{\mathit{a}}{\mathit{r}}\left(\cos \mathit{E}\mathrm{-}\mathit{e}\mathrm{+}\mathrm{i}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\sin \mathit{E}\right)\mathit{,}\end{array}$we obtain ${\mathit{X}}_{\mathit{k}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\left(\cos \mathit{E}\mathrm{-}\mathit{e}\mathrm{+}\mathrm{i}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\sin \mathit{E}\right){\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k}\mathrm{(}\mathit{E}\mathrm{-}\mathit{e}\sin \mathit{E}\mathrm{)}}\mathrm{d}\mathit{E}\mathit{.}$(A.5)From this last expression, one only has to expand $\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}$ 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): ${\mathit{X}}_{\mathit{k}}\mathrm{\approx }\frac{\mathrm{1}}{\mathit{N}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\left(\cos {\mathit{E}}_{\mathit{k}}\mathrm{-}\mathit{e}\mathrm{+}\mathrm{i}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\sin {\mathit{E}}_{\mathit{k}}\right){\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k}\mathrm{(}{\mathit{E}}_{\mathit{k}}\mathrm{-}\mathit{e}\sin {\mathit{E}}_{\mathit{k}}\mathrm{)}}\mathit{,}$(A.6)with ${\mathit{E}}_{\mathit{k}}\mathrm{=}\mathrm{2}\mathit{\pi }\frac{\mathit{k}\mathrm{-}\mathrm{1}}{\mathit{N}}$.

Fig. A.1 Convergence of the numerical estimation of X1(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.

Figure A.1 illustrates how the numerical estimates of Eq. (A.6) evolves as a function of N (for X₁(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 $\begin{array}{c}\\ \frac{\mathrm{d}{\mathit{X}}_{\mathit{k}}}{\mathrm{d}\mathit{e}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}& [\mathrm{i}\mathit{k}\sin \mathit{E}\left(\cos \mathit{E}\mathrm{-}\mathit{e}\mathrm{+}\mathrm{i}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\sin \mathit{E}\right)\\ & \mathrm{-}\left(\mathrm{1}\mathrm{+}\mathrm{i}\frac{\mathit{e}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}}\sin \mathit{E}\right)]{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k}\mathrm{(}\mathit{E}\mathrm{-}\mathit{e}\sin \mathit{E}\mathrm{)}}\mathrm{d}\mathit{E}\mathit{.}\end{array}$(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₀. 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 $\begin{array}{ccc}{\mathit{V}}_{\min }& \mathrm{=}& \mathit{K}\mathrm{(}\mathit{e}\cos \mathit{\omega }\mathrm{-}\mathrm{1}\mathrm{)}\mathit{,}\\ {\mathit{V}}_{\max }& \mathrm{=}& \mathit{K}\mathrm{(}\mathit{e}\cos \mathit{\omega }\mathrm{+}\mathrm{1}\mathrm{)}\mathit{,}\\ {\mathit{v}}_{\min }& \mathrm{=}& \mathit{\pi }\mathrm{-}\mathit{\omega ,}\\ {\mathit{v}}_{\max }& \mathrm{=}& \mathrm{-}\mathit{\omega }\mathit{.}\end{array}$We have (see Eqs. (C.1), (C.2)) $\begin{array}{ccc}& & \mathit{K̂}\mathrm{=}\frac{{\mathit{V}}_{\max }\mathrm{-}{\mathit{V}}_{\min }}{\mathrm{2}}\mathit{,}\\ & & \mathit{ê}\cos \mathit{\omega ̂}\mathrm{=}\frac{{\mathit{V}}_{\max }\mathrm{+}{\mathit{V}}_{\min }}{{\mathit{V}}_{\max }\mathrm{-}{\mathit{V}}_{\min }}\mathrm{·}\end{array}$To determine the remaining parameters (esinω, M₀), we need to use the information contained in the timing of extrema. We have $\mathit{M}\mathrm{=}{\mathit{M}}_{\mathrm{0}}\mathrm{+}\mathit{nt,}$(C.7)thus $\begin{array}{ccc}{\mathit{M}}_{\min }& \mathrm{=}& {\mathit{M}}_{\mathrm{0}}\mathrm{+}\mathit{n}{\mathit{t}}_{\min }\mathit{,}\\ {\mathit{M}}_{\max }& \mathrm{=}& {\mathit{M}}_{\mathrm{0}}\mathrm{+}\mathit{n}{\mathit{t}}_{\max }\mathit{.}\end{array}$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) $\begin{array}{ccc}\mathit{M}& \mathrm{=}& \mathit{E}\mathrm{-}\mathit{e}\sin \mathit{E,}\\ \mathit{E}& \mathrm{=}& \mathrm{2}\arctan \left(\sqrt{\frac{\mathrm{1}\mathrm{-}\mathit{e}}{\mathrm{1}\mathrm{+}\mathit{e}}}\tan \left(\frac{\mathit{v}}{\mathrm{2}}\right)\right)\mathrm{·}\end{array}$Developing Eqs. (C.10), (C.11) at second order in eccentricity, we obtain $\mathit{M}\mathrm{=}\mathit{v}\mathrm{-}\mathrm{2}\mathit{e}\sin \mathit{v}\mathrm{+}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{e}}^{\mathrm{2}}\sin \mathrm{2}\mathit{v}\mathrm{+}\mathrm{O}\mathrm{\left(}{\mathit{e}}^{\mathrm{3}}\mathrm{\right)}\mathit{,}$(C.12)which gives

$\begin{array}{ccc}& & {\mathit{M}}_{\min }\mathrm{=}\mathit{\pi }\mathrm{-}\mathit{\omega }\mathrm{-}\mathrm{2}\mathit{e}\sin \mathit{\omega }\mathrm{-}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{e}}^{\mathrm{2}}\sin \mathrm{2}\mathit{\omega }\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{3}}\mathrm{)}\mathit{,}\\ & & {\mathit{M}}_{\max }\mathrm{=}\mathrm{-}\mathit{\omega }\mathrm{+}\mathrm{2}\mathit{e}\sin \mathit{\omega }\mathrm{-}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{e}}^{\mathrm{2}}\sin \mathrm{2}\mathit{\omega }\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{3}}\mathrm{)}\mathit{,}\\ & & {\mathit{M}}_{\max }\mathrm{-}{\mathit{M}}_{\min }\mathrm{=}\mathit{\pi }\mathrm{+}\mathrm{4}\mathit{e}\sin \mathit{\omega }\mathrm{+}\mathrm{O}\mathrm{(}{\mathit{e}}^{\mathrm{3}}\mathrm{)}\mathit{,}\end{array}$thus (see Eqs. (C.8), (C.9)) $\mathit{ê}\sin \mathit{\omega ̂}\mathrm{=}\frac{\mathit{n}}{\mathrm{4}}\mathrm{(}{\mathit{t}}_{\max }\mathrm{-}{\mathit{t}}_{\min }\mathrm{)}\mathrm{-}\frac{\mathit{\pi }}{\mathrm{4}}\mathrm{mod}\hspace{0.17em}\hspace{0.17em}\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{·}$(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₀ is straightforwardly derived from these values $\mathit{M̂}\mathrm{0}\mathrm{=}{\mathit{M}}_{\min }\mathrm{-}\mathit{n}{\mathit{t}}_{\min }\mathrm{=}{\mathit{M}}_{\max }\mathrm{-}\mathit{n}{\mathit{t}}_{\max }\mathit{.}$(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

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 (λ = M0 + ω), the semi-amplitude 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 Newton-Raphson correcting algorithm (see Sect. 3.4).

In the text

	Fig. 2 Radial-velocity data of GJ 3021 (top) and corresponding periodogram (bottom). The data are taken from Naef et al. (2001; 61 CORALIE-98 measurements).
In the text

	Fig. 3 Radial-velocity 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).
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 CORALIE-98 measurements (red), and 11 CORALIE-07 measurements (blue).
In the text

	Fig. 5 Phase-folded radial-velocity 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).
In the text

	Fig. 6 Radial-velocity data of HD 192310. The data are taken from Pepe et al. (2011; 139 HARPS measurements).
In the text

	Fig. 7 Periodogram of the radial-velocities of HD 192310 (A)), of the residuals of the one-Keplerian model (LM1 see Table 3, B)), and of the residuals of the two-Keplerians 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).
In the text

	Fig. 8 Phase folded radial-velocity curves of HD 192310b (top), and c (bottom), superimposed with our final two-Keplerians model (LM2). See Table 3 for the corresponding orbital elements.
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 CORALIE-98 measurements (red) and 95 CORALIE-07 measurements (blue).
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).
In the text

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

	Fig. A.1 Convergence of the numerical estimation of X1(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.
In the text