Dynamical stability analysis of the HD 202206 system and constraints to the planetary orbits

J. Couetdic; J. Laskar; A. C. M. Correia; M. Mayor; S. Udry

doi:10.1051/0004-6361/200913635

Home

All issues

Volume 519 (September 2010)

A&A, 519 (2010) A10

Full HTML

Free Access

Issue		A&A Volume 519, September 2010


Article Number		A10
Number of page(s)		14
Section		Planets and planetary systems
DOI		https://doi.org/10.1051/0004-6361/200913635
Published online		07 September 2010

A&A 519, A10 (2010)

Dynamical stability analysis of the HD 202206 system and constraints to the planetary orbits

J. Couetdic¹ - J. Laskar¹ - A. C. M. Correia^1,2 - M. Mayor³ - S. Udry³

1 - Astronomie et Systèmes Dynamiques, IMCCE-CNRS UMR 8028, Observatoire de Paris, UPMC, 77 Avenue Denfert-Rochereau, 75014 Paris, France
2 - Departamento de Física, Universidade de Aveiro, Campus de Santiago, 3810-193 Aveiro, Portugal
3 - Observatoire de Genève, 51 ch. des Maillettes, 1290 Sauverny, Switzerland

Received 10 November 2009 / Accepted 24 April 2010

Abstract
Context. Long-term, precise Doppler measurements with the CORALIE spectrograph have revealed the presence of two massive companions to the solar-type star HD 202206. Although the three-body fit of the system is unstable, it was shown that a 5:1 mean motion resonance exists close to the best fit, where the system is stable. It was also hinted that stable solutions with a wide range of mutual inclinations and low O-C were possible.
Aims. We present here an extensive dynamical study of the HD 202206 system, aiming at constraining the inclinations of the two known companions, from which we derive possible value ranges for the companion masses.
Methods. We consider each inclination and one of the longitudes of ascending node as free parameters. For any chosen triplet of these parameters, we compute a new fit. Then we study the long-term stability in a small (in terms of O-C) neighborhood using Laskar's frequency map analysis. We also introduce a numerical method based on frequency analysis to determine the center of libration mode inside a mean motion resonance.
Results. We find that acceptable coplanar configurations (with low $\chi ^2$ stable orbits) are limited with respect to inclinations to the line of sight between $30\hbox{$^\circ$ }$ and $90\hbox {$^\circ $ }$ . This limits the masses of both companions to roughly twice the minimum: $m_{\rm b} \in [16.6~M_{\rm Jup};~33.5~M_{\rm Jup}]$ and $m_{\rm c} \in [2.2~M_{\rm Jup};~4.4~M_{\rm Jup}]$ . Non-coplanar configurations are possible for a wide range of mutual inclinations from $0\hbox {$^\circ $ }$ to $90\hbox {$^\circ $ }$ , although $\Delta\Omega = 0 [\pi]$ configurations seem to be favored. We also confirm the 5:1 mean motion resonance to be most likely. In the coplanar edge-on case, we provide a very good stable solution in the resonance, whose $\chi ^2$ does not differ significantly from the best fit. Using our method for the determination of the center of libration, we further refine this solution to obtain an orbit with a very low amplitude of libration, as we expect that dissipative effects have dampened the libration.

Key words: stars: individual: HD 202206 - planetary systems - methods: numerical - techniques: radial velocities - celestial mechanics

1 Introduction

The CORALIE planet-search program in the southern hemisphere has found two companions around the HD 202206 star. The first one is a very massive body with $17.5~M_{\rm jup}$ minimum mass (Udry et al. 2002), while the second companion is a $2.4~M_{\rm jup}$ minimum mass planet (Correia et al. 2005). The parent star has a mass of 1.044 solar masses and is located 46.3 pc from the Solar System. The HD 202206 planetary system is an interesting case for investigating the brown dwarf desert since the more massive companion can be either a huge planet (formed in the circumstellar disk) or a low-mass brown dwarf candidate.

Correia et al. (2005) found that the orbital parameters obtained with the best fit for the two planets leads to catastrophic events in a short time (two Keplerians fit and full three-body fit alike). This was not completely unexpected given the very high eccentricities (0.435 and 0.267) and masses of the two planets. Using frequency analysis (Laskar 1993,1990,1999), they performed a study of the global dynamics around the best fit and found that the strong gravitational interactions with the first companion made the second planet evolution very chaotic, except for initial conditions in the 5:1 mean motion resonance. Since the associated resonant island actually lies close to the $\chi ^2$ minimum value of the best fit, they concluded that the system should be locked in this 5:1 resonance. Later on, Gozdziewski and co-workers also looked for stable solutions in this system using their GAMP algorithm (Gozdziewski et al. 2006). They provide two possible solutions (among many others), one coplanar and one with a very high mutual inclination.

Since then, new data have been acquired using the CORALIE spectrograph. A new reduction of the data changed some parameters, including the mass of the HD 202206 star. A new fit assuming a coplanar edge-on configuration was derived from the new set of radial velocity data. This new solution still appears to be in the 5:1 mean motion resonance, but is also still unstable. The most striking difference from Correia et al. (2005) is the lowerer eccentricity of HD 202206 c.

In the present work, we continue the dynamical study started in Correia et al. (2005) in more detail, using the new fit as a starting point. We also aim to find constraints on the orbital parameters of the two known bodies of this system, in particular the inclinations (and thus the real masses of the planets). Gozdziewski et al. (2006) have already shown that stable fits could be obtained with different inclinations, using a particular fitting genetic algorithm that adds stability computation to select its populations (GAMP). Although very effective in finding a stable fit, this algorithm cannot find all possible solutions. We prefer an approach here that separates the fitting procedure from dynamical considerations, bacause it allows for a better assessment of the goodness of the fit and of whether the model is a good description of the available data. The trade-off is a more difficult handling of the large number of parameters.

We briefly present the new set of data in Sect. 2 and the numerical methodology in Sect. 3. We review in detail the dynamics in the coplanar edge-on case in Sect. 4. We then release the constraint on the inclination of the system from the line of sight in Sect. 5, and finally briefly investigate the mutually inclined configurations in Sect. 6.

2 New orbital solution for the HD 202206 system

The CORALIE observations of HD 202206 started in August 1999 and the last point acquired in our sample dates from September 2006, corresponding to about seven years of observations and 92 radial-velocity measurements. Using the iterative Levenberg-Marquardt method (Press et al. 1992), we fit the observational data using a 3-body Newtonian model, assuming co-planar motion perpendicular to the plane of the sky, similar to what has been done in (Correia et al. 2009,2005). We changed the reference date with respect to the solution in Correia et al. (2005). The mass of the star has also been updated to $1.044~M_\odot$ (Sousa et al. 2008). This fit yields two planets with an adjustment of $\sqrt{\chi_{\rm r}^2}=1.41$ and $rms=7.45~{\rm m/s}$ , slightly above the photon noise of the instrument, which is around $6.69~{\rm m/s}$ . We confirm the already detected planets (Udry et al. 2002; Correia et al. 2005) with improved orbital parameters, one at P=254.8 day, e=0.431, and a minimum mass of $16.6~M_{\rm Jup}$ , and the other at P=1397 day, e=0.104, and a minimum mass of $2.18~M_{\rm Jup}$ (Table 1). In Fig. 1 we plot the observational data superimposed on the best fitted solution.

Table 1: Best Newtonian fit S1 for the HD 202206 system assuming $i_{\rm b}=i_\sc=90\hbox{$^\circ$ }$ and $\Delta \Omega =0\hbox {$^\circ $ }$ .

$\begin{figure} \par\includegraphics[angle=270,width=8.8cm,clip]{figs/13635_fig01} \end{figure}$	Figure 1: CORALIE radial velocities for HD 202206 superimposed on a 3-body Newtonian orbital solution (Table 1).
Open with DEXTER

We also fitted the data with a 3-body Newtonian model for which the inclination of the orbital planes, as well as the node of the outer planet orbit were free to vary. We were able to find a wide variety of configurations, some with low inclination values for one or both planets, which slightly improved our fit to a minimum $\sqrt{\chi_{\rm r}^2}= 1.30$ and rms = $7.08~{\rm m/s}$ . However, all of these determinations remain uncertain, and since we also increased the number of free parameters by three, we cannot say that there has been an improvement over the solution presented in Table 1.

3 Numerical set-up

3.1 Conventions

In this paper, the subscripts ${\rm b}$ and $\sc$ respectively refer to the bodies with the shortest and longest orbital periods (inner and outer). The initial conditions ill-constrained by the radial velocity data are $i_{\rm b},i_\sc,\Omega_{\rm b},\Omega_\sc$ (the inclinations and longitude of ascending nodes). We are using the observers' convention that sets the plane of sky as the reference plane (see Fig. 2). As a consequence, the nodal line is in the plane of sky and has no cinematic impact on the radial velocities. From the dynamical point of view, only the difference between the two lines of nodes $\Delta \Omega$ matters. In particular, the mutual inclination depends on this quantity, through

$\begin{displaymath}% \cos{I} = \cos{i_{\rm b}}\cos{i_\sc} + \sin{i_{\rm b}}\sin{i_\sc}\cos{\Delta\Omega}. \end{displaymath}$

(1)

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig02} \vspace*{5mm} \end{figure}$	Figure 2: Angles defining the orbit's orientation in space. We follow the observers' convention that sets the plane of sky as the reference frame, for which the edge-on coplanar configuration is $i_{\rm b}=i_\sc=90\hbox{$^\circ$ }$ and $\Delta \Omega =0\hbox {$^\circ $ }$ . The observer is limited to the velocity projected on the k axis.
Open with DEXTER

Thus, $\Omega_{\rm b}$ can always be set to $0\hbox {$^\circ $ }$ in the initial conditions, which leads to $\Delta\Omega = \Omega_\sc$ , and only three parameters are left free $(i_{\rm b},i_\sc,\Omega_\sc)$ . They are connected to the three interesting unknowns of the system, namely the mutual inclination I (Eq. (1)), and the two planetary masses $m_{\rm b}$ , $m_\sc$ as

$\begin{displaymath}% \left\{ \begin{array}{ll} \displaystyle\frac{m_{\rm b}(m_\a... ...rt{1-e_\sc^2}}{(2\pi G)^{1/3} \sin{i_\sc}}, \end{array}\right. \end{displaymath}$

(2)

where $m_\ast$ is the star mass, $K_\sigma$ the amplitude of the radial velocity variations, $P_\sigma$ the orbital period, $e_\sigma$ the eccentricity, and G the gravitational constant. Basically, choosing the values of the inclinations is akin to setting the two companions masses. And for given values of $(i_{\rm b},i_\sc)$ , we can control the mutual inclination with $\Omega_\sc$ through Eq. (1).

The denominator on the left hand sides in Eq. (2) is the total mass of the system. This term comes from the transformation to the barycentric coordinates system. Indeed, the elliptical elements are astrocentric elliptical elements, while the observed variations in radial velocity (and thus $K_\sigma$ ) are considered in the barycentric reference frame. As a consequence, the two equations are coupled. Of course, we can usually neglect the companion masses in this term, and decouple the equations; however, when we change the inclinations, the planetary masses will grow to a point where this approximation is no longer valid. Here we always solve the complete equations, regardless of the inclinations.

As long as the companions' masses are small compared to the primary, they are scaled to $1/\sin{i}$ to a good approximation. We can define two factors $k_{\rm b}= 1/\sin{i_{\rm b}}$ and $k_\sc = 1/\sin{i_\sc}$ . With $m_{\rm b}^{(0)}$ and $m_\sc^{(0)}$ the minimum masses obtained for the edge-on coplanar case (see Sect. 4 and Table 1), we can write

	$\displaystyle m_{\rm b}\approx k_{\rm b}m_{\rm b}^{(0)} = \frac{1}{\sin{i_{\rm b}}} m_{\rm b}^{(0)},$
	$\displaystyle m_\sc \approx k_\sc m_\sc^{(0)} = \frac{1}{\sin{i_\sc}} m_\sc^{(0)}.$	(3)

For a given factor k, two values of the inclination are possible: x and $\pi - x$ (where $x\in[0;\pi/2]$ ). For instance, $i=30\hbox {$^\circ $ }$ and $i=150\hbox{$^\circ$ }$ give k=2.

Additionally, for a given pair $(i_{\rm b},i_\sc)$ , the accessible mutual inclinations I are limited. Since the inclinations $i_{\rm b}$ and $i_\sc$ are angles between 0 and $\pi$ (excluded), prograde coplanar configurations are only possible for $i_{\rm b}= i_\sc$ and $\Delta\Omega = 0$ (Eq. (1)). Similarly, retrograde coplanar configurations are only possible for $i_{\rm b}+ i_\sc = \pi$ and $\Delta\Omega = \pi$ . This means that $k_{\rm b}= k_\sc$ in both cases. More generally, $I \ge \vert i_{\rm b}- i_\sc\vert$ , and

	$\displaystyle i_{\rm b}+ i_\sc \le \pi \ \Longrightarrow \ I \le i_{\rm b}+ i_\sc,$
	$\displaystyle i_{\rm b}+ i_\sc > \pi \ \Longrightarrow \ I \le 2\pi - i_{\rm b}- i_\sc.$	(4)

For any value of I, two values of $\Delta \Omega$ are possible: x and -x, since $\Delta \Omega$ appears through its cosine in Eq. (1). The extrema are obtained for $\Delta\Omega = 0$ (maximum), and $\Delta\Omega = \pi$ (minimum). Finally we notice that one can restrict one of the two inclinations to the line of sight to $]0;\pi/2]$ , which will be the case for $i_{\rm b}$ .

3.2 Fitting procedure

The influence of $i_{\rm b}$ , $i_{\rm c}$ , and $\Omega_{\rm c}$ on the radial velocity data is usually very small, and those perturbations that depend on them have very long time scales. This makes any attempt to fit those parameters virtually impossible at present. Only with very strong mean motion resonances, such as observed in the GJ 876 system (Correia et al. 2009; Laughlin & Chambers 2001), can one hope to fit the inclinations. In this case the mean motion resonance introduces important short time-scale terms (compared with the precision and time span of the observations). For the HD 202206 system, the set of radial velocity data does not cover a long enough period of time.

Since the three parameters $i_{\rm b}$ , $i_\sc$ , and $\Omega_\sc$ are very poorly constrained by the radial velocity data, we cannot fit the data with a model that includes them. Instead, for any chosen ( $i_{\rm b}$ , $i_\sc$ , $\Omega_\sc$ ) set, we compute a new best fit using a Levenberg-Marquard minimization (Press et al. 1992) and a three-body model, but with eleven free parameters: the center of mass velocity $\gamma$ , and for each planet, the semi-amplitude of radial velocity K, the period P, eccentricity e, mean anomaly M, and periastron $\omega$ , all given at the initial epoch. Throughout the paper the initial conditions are given at the same initial epoch $T_0=2~453~000\ {\rm JD}$ .

At this point we have a complete description of the system with the mass of the hosting star and 14 parameters (7 for each planets) as follows:

4 chosen: $i_\sigma$ and $\Omega_\sigma$ ;
and 10 fitted: $K_\sigma$ , $P_\sigma$ , $e_\sigma$ , $M_\sigma$ , and $\omega_\sigma$ .

However, masses and elliptic coordinates are easier to manipulate for a dynamical study. Using Eq. (2), we obtain a system of two equations with two unknowns ( $m_{\rm b}$ and $m_\sc$ ), which is easily solved with a Newton algorithm. The semi-major axis are then obtained using Kepler's third law.

3.3 Numerical integrations

For the numerical integrations, we use the Newtonian equations with secular corrections for the relativity. The Newtonian part of the integration is carried out by the symplectic integrator SABAC4 of Laskar & Robutel (2001) with a step size of 0.02 year. The secular corrections for the relativity are computed from the perturbation formulae given in Lestrade & Bretagnon (1982)

	$\displaystyle \displaystyle\frac{{\rm d}a}{{\rm d}t} = \mu\frac{2na^2e\sin{v}}{r^2\eta} \left[10\frac{a}{r} - 3\right],$
	$\displaystyle \displaystyle\frac{{\rm d}\lambda}{{\rm d}t} = n + \mu \displayst... ...r^2} \left\{8\displaystyle\frac{a\eta^2}{r}+6\displaystyle\frac{r}{a}-20\right.$
	$\displaystyle \qquad \left. -\displaystyle\frac{\eta^2-\eta}{e^2} \left[ -10\di... ...-\displaystyle\frac{a\eta^2}{r} + 18\displaystyle\frac{a}{r}-7 \right]\right\},$
	$\displaystyle \displaystyle\frac{{\rm d}\varpi}{{\rm d}t} = \mu \displaystyle\f... ...}}{r^2} -\displaystyle\frac{a\eta^2}{r} + 18\displaystyle\frac{a}{r}-7 \right],$
	$\displaystyle \displaystyle\frac{{\rm d}e}{{\rm d}t} = \mu \displaystyle\frac{na\eta\sin{v}}{r^2} \left[ 10 \displaystyle\frac{a}{r} - 7 \right],$	(5)

where $\mu = Gm_\ast/c^2$ , n is the mean motion, and $\eta = \sqrt{1-e^2}$ . These equations are averaged to obtain the following first-order secular perturbations

	$\displaystyle \langle \dot{a} \rangle = 0,$
	$\displaystyle \langle \dot{e} \rangle = 0,$
	$\displaystyle \langle \dot{M} \rangle = \displaystyle\frac{\mu n}{a} \left[ 6-\displaystyle\frac{15}{\sqrt{1-e^2}} \right],$
	$\displaystyle \langle \dot{\varpi} \rangle = \displaystyle\frac{3\mu n}{a(1-e^2)}\cdot$	(6)

These corrections are computed every 100 steps, that is, every two years with the current values at the given step for e, a, and n, for each planet. These approximated equations have been successfully tested by comparison with INPOP (Fienga et al. 2008).

3.4 Stability threshold

To study the stability of any given orbit, we use Laskar's frequency map analysis (Laskar 1993,1990). Using a numerical integration of the orbit over a time interval of length T, we compute a refined determination (in $\hbox{$^\circ$ }/{\rm yr}$ ) of the mean motions n₁, $n_1^\prime$ obtained over two consecutive time intervals of length T₁ = T/2. The stability index $D_1 = \vert n_1^\prime - n_1\vert$ provides a measure of the chaotic diffusion of the trajectory. Low values close to zero correspond to a regular solution, while high values are synonymous with strong chaotic motion (Laskar 1993).

In this paper we looked at many different orbits for many different initial conditions to detect stable regions. This calls for a way to automatically calibrate a threshold for stability $D_{\rm lim}$ . To that end, we used a second stability index D₂. Using the same numerical integration, we computed two new determinations of the mean motion n₂ and $n_2^\prime$ over two consecutive time intervals of length T₂ = T₁/k, where k > 1. In the case of quasi-periodic motion, the diffusion should be close to zero but it is limited by the precision in determining the frequencies. Since D₁ is computed over longer time intervals, the frequencies are better determined, and thus, D₁ should be, on average, smaller than D₂. On the contrary, for chaotic trajectories, the diffusion will increase on average for longer time intervals. We can then determine an approximated value $D_{\rm lim}$ for which

$\begin{displaymath}% \left\{\begin{array}{ll} D_1 \gg D_{\rm lim} &\Longrightar... ...D_{\rm lim} &\Longrightarrow D_1 < D_2.\\ \end{array}\right. \end{displaymath}$

(7)

We are then assured that the first kind of orbits is in general stable, while the latter is considered chaotic.

This approach is best used statistically over a grid of initial conditions, especially when we try to use a small integration time. To determine $D_{\rm lim}$ for a particular diffusion grid, we look at the distribution of D₁ < D₂ trajectories as a function of D₁. We actually work with smoothed values $D_1^{\rm s}$ and $D_2^{\rm s}$ of D₁ and D₂ to reduce the influence of the chaotic orbits whose mean motion diffusions are small by mere chance. They appear as low diffusion orbits inside a high diffusion region. The smoothing function is a simple geometric mean over the closest neighbors. Other functions, such as a convolution with 2D Gaussian, have been tested, but do not yield significantly better results.

We bin the $\log{D_1^{\rm s}}$ data in 0.5 wide bins, and compute the percentage of $D_1^{\rm s} < D_2^{\rm s}$ orbits for each bin. Figure 3 shows a typical distribution obtained from a diffusion grid (in this case the top panel of Fig. 6). It reproduces the behavior expected from Eq. (7): low diffusion orbits tend, for the majority, to have their diffusion index diminish when time increases.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig03} \end{figure}$

Figure 3:

Distribution of D₁ < D₂ trajectories from the top panel of Fig. 6. Each integrated trajectory is binned with respect to its diffusion index $\log{D_1}$ , after the diffusion grid has been smoothed. For each bin, we computed the proportion of trajectories with a decreasing diffusion index over time: D₁ < D₂. The histogram shows the results for 0.5 wide bins. As expected from Eq. (7), orbits with a high diffusion index D₁, and D₁ < D₂ are nearly non existent ( $\log{D_1} > -1$ ), while we observe the opposite situation for low-diffusion index orbits ( $\log{D_1} < -3$ ).

Open with DEXTER

We choose to define $D_{\rm lim}$ as the D₁ value for which $99\%$ of the trajectories exhibit D₁ < D₂. Graphically, it is the abscissa for which the curve in Fig. 3 crosses y = 0.99. In this example we get $\log{D_{\rm lim}} \approx -2.87$ .

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig04} \end{figure}$

Figure 4:

Percentage of orbits wrongly flagged asstable (false stable) or unstable (false unstable). Once $D_{\rm lim}$ is chosen using a given percentage threshold (see Fig. 3), we have a criterion for stable ( $D_1 < D_{\rm lim}$ ) and unstable orbits ( $D_1 > D_{\rm lim}$ ). We compared the results with a diffusion grid computed over a longer time interval (in this case 2 $\times$ 40 000 years) taken as a reference. A few orbits deemed stable from the reference diffusion grid were thought to be unstable and vice versa. The solid line traces the number of those faulty orbits for different values of the threshold. It appears to be lowest between 0.95 and 1. The dotted line traces the false positives, and the dashed line the false negatives. The solid one is the sum of those two.

Open with DEXTER

The $99\%$ threshold is actually a compromise that works in the majority of encountered cases: it minimizes the number of orbits wrongly flagged as stable (false stable) or unstable (false unstable). This is also approximately the value for which the number of false stables and false unstables is equivalent. For percentages lower than $95\%$ , the number of false unstables is nearly zero. This is easy to understand, since a higher percentage threshold implies a lower value for $D_{\rm lim}$ , which in turn leads to flagging very few actually unstable orbits, while we might miss several stable ones, and vice-versa. To estimate the number of false stables and unstables, we recomputed the diffusion grid on a longer integration time, 2 $\times$ 40 000 years, and used this grid as a reference (Fig. 4).

4 Review of the coplanar edge-on case

4.1 Global dynamics

Following Correia et al. (2005), we study in more detail the dynamics in the neighborhood of the 3-body fit obtained in the case of coplanar orbits with $\sin i_{\rm b}= \sin i_\sc = 1$ (that is, the system seen edge-on). The best fit to the radial velocity data for this particular configuration is given in Table 1. It is different from the solution S4 in Correia et al. (2005, Table 4) as explained in Sect. 2. For the dynamical aspect of the system, the important change is the decrease in planet c's eccentricity. As a consequence regions outside resonances are expected to be more stable, and the environment of the fit should be less chaotic. However this new solution is still unstable and the outer planet is lost shortly after about 150 million years.

We look for possible nearby stable zones, keeping HD 202206 b parameters constant since they are much better constrained, with small standard errors. We assume for now that the system is coplanar and seen edge-on, that is, with both inclination at $90\hbox {$^\circ $ }$ and $\Omega_\sc = 0\hbox{$^\circ$ }$ . We let $a_\sc$ , $e_\sc$ , and $\omega_\sc$ vary. We always keep $M_\sc + \omega_\sc$ constant, as it is much better constrained by the radial velocity data. This implies that when we change $\omega_\sc$ , the mean anomaly $M_\sc$ varies accordingly. In the particular case where $\Omega_\sc = 0\hbox{$^\circ$ }$ , this means that the initial mean longitude $\lambda_\sc$ is kept constant. For each initial conditions, we compute the diffusion index $\log{D_1}$ and the square root of the reduced $\chi ^2$ .

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig05} \end{figure}$

Figure 5:

Global view of the dynamics of HD 202206 for variations of the semi-major axis and periastrum of the outer planet ( bottom panel) or semi-major axis and eccentricity ( top panel). The step sizes for $a_\sc$ , $\omega_\sc$ , $e_\sc$ are $0.005~{\rm AU}$ , $2\hbox {$^\circ $ }$ , and 0.004, respectively. The other parameters were kept constant and taken from the fit S1 (Table 1). The color scale is the stability index $\log{D_1}$ obtained through a frequency analysis of the longitude of the outer planet over two consecutive time intervals of 8000 years. The level curves give the $\sqrt {\chi _{\rm r}^2}$ value computed for each choice of parameters. The two horizontal black lines mark the intersection of the two grids. Most of the orbits are chaotic (yellow to red dots), but regular orbit zones exist (blue dots). One of these regular orbits regions lies inside the low $\sqrt {\chi _{\rm r}^2}$ region: the dark area around $a_\sc=2.5$ AU and $\omega_\sc=50\hbox{$^\circ$ }$ . It corresponds to the stable island of the 1/5 mean motion resonance. The cross marks the S1 orbital solution.

Open with DEXTER

Figure 5 shows a global picture of the dynamics around the fit, in the planes $(a_\sc,e_\sc)$ and $(a_\sc,\omega_\sc)$ of initial conditions. The step sizes for $a_\sc$ , $\omega_\sc$ , $e_\sc$ are $0.005~{\rm AU}$ , $2\hbox {$^\circ $ }$ , and 0.004 respectively. The other parameters were kept constant and taken from the fit S1 (Table 1). The level curves give the $\sqrt {\chi _{\rm r}^2}$ value computed for each set of initial conditions. The color scale gives the diffusion index $\log{D_1}$ . The yellow to red areas are very chaotic, mainly owing the high eccentricities and masses of both planets.

The orbital solution S1 lies inside the $\sqrt {\chi _{\rm r}^2}= 1.5$ level curves, at the coordinates marked by a cross, inside a high-diffusion (green) area. Several low-diffusion (blue and dark blue) zones exist for which the orbits are stabilized either by mean motion resonances or by locking of $\Delta\varpi$ around $0\hbox {$^\circ $ }$ . Orbits stabilized by the corotation of the apsidal lines are the blue to black zones around $\omega_\sc = 190\hbox{$^\circ$ }$ (bottom panel). The width (in the $\omega_\sc$ direction) increases with the semi-major axis $a_\sc$ , from 90 degrees at 2 AU, to nearly 360 degrees at 4 AU, since a wider libration of $\Delta\varpi$ around $0\hbox {$^\circ $ }$ is possible without close encounters when the distance between the two planets increases. The red to green more or less vertical stripes cutting through these zones mark mean motion resonances, which for the most part have a destabilizing effect in the two sections of the phase space considered in Fig. 5. However, the stronger resonances also have stable orbits:

the 1/4 MMR at 2.2 AU with two stable islands around $\omega_\sc = 260\hbox{$^\circ$ }$ and $\omega_\sc=50\hbox{$^\circ$ }$ ;
the 5:1 MMR at 2.5 AU with a notable stable island around $\omega_\sc = 70\hbox{$^\circ$ }$ where the best fit is located;
the 1/6 MMR at 2.8 AU with a stable island around $\omega_\sc = 0\hbox{$^\circ$ }$ or high eccentricity.

4.2 Stable fit

We now take a closer look at the 5:1 mean motion resonance island around $a_\sc=2.5$ AU and $\omega_\sc=50\hbox{$^\circ$ }$ , where we believe the system is presently locked. Figure 6 was constructed the same way as Fig. 5, but the step sizes are now 0.0025 AU for $a_\sc$ , $0.5\hbox {$^\circ $ }$ for $\omega_\sc$ , and 0.002 for $e_\sc$ .

For eccentricities higher than 0.2, orbits are very chaotic (red dots), as the outer planet undergoes close encounters with the inner body. At lower eccentricities we notice some lower diffusion orbits for $\omega_\sc > 120\hbox{$^\circ$ }$ . Those orbits lie far outside the resonance, but may be stable because of the low eccentricity of the outer planet and the apsidal locking mechanism. They are however too far from the best fit $(\!\!\sqrt{\chi_{\rm r}^2} > 1.5)$ and are less likely to be a good guess of the actual configuration of HD 202206 system.

A very noticeable feature of this resonant island appearing in both panels is the existence of two distinct stable regions, separated by chaotic orbits inside the resonance itself. The two stable regions actually correspond to two different critical arguments: $\lambda_{\rm b}- 5\lambda_\sc + 4\varpi_\sc$ in the structure on the rim, and $\lambda_{\rm b}- 5\lambda_\sc + \varpi_{\rm b}+ 3\varpi_\sc$ in the center.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig06} \end{figure}$

Figure 6:

Global view of the dynamics of HD 202206 for variations of the semi-major axis and periastrum ( bottom), and semi-major axis and eccentricity ( top) of the outer planet. The step sizes are respectively 0.0025 AU, $0.5\hbox {$^\circ $ }$ , and 0.002 for the eccentricity. The color scale is the stability index ( $\log{D_1}$ ), and the level curves give the $\sqrt {\chi _{\rm r}^2}$ values (see Fig. 5). The cross marks the best fit S1, and the horizontal line in each panel, the orbits common to both maps (i.e. the intersection of each map). The best fit lies very close to a stable resonant island where we can pick a stable solution with a $\chi ^2$ only marginally higher than that of the best fit. As an example, we picked such a solution (S2), marked by a white filled circle, by slightly increasing $a_{\rm c}$ from 2.4832 AU to 2.49 AU. The complete set of orbital elements for S2 is given in Table 2.

Open with DEXTER

The orbital solution S1 lies very close to the latter, in a chaotic region (green and yellow dots), between the two stable parts of the resonant island. In fact, one can pick stable orbits with a $\sqrt {\chi _{\rm r}^2}$ not significantly worse than the best fit. For instance, the orbital solution S2 given in Table 2 is stable and has a $\sqrt {\chi _{\rm r}^2}$ of 1.4136 (marked by a filled white circle). The orbital elements are the same as S1, except for $a_\sc$ which was adjusted from 2.4832 AU to 2.49 AU.

Table 2: Stable orbital parameters S2 for the HD 202206 system for $i_{\rm b}=i_\sc=90\hbox{$^\circ$ }$ and $\Delta \Omega =0\hbox {$^\circ $ }$ .

4.3 Resonant and secular dynamics

The orbital solution S2 was integrated over $5~{\rm Gyr}$ . It remained stable and displays a regular behavior during the whole time. Using frequency analysis on an integration over $1~{\rm Myr}$ , we determined its fundamental frequencies (Table 3). Following our notation, $n_{\rm b}$ and $n_\sc$ are the mean motions. The secular frequencies are denoted g₁ and g₂. Finally $l_\theta$ is the frequency associated with the resonance's critical angle $\theta$ . The fundamental secular frequencies g₁ and g₂, related to the periastron of the inner and outer planets correspond to the periods $P_1 \approx 14.1$ $\times$ 10³ yr and $P_2 \approx - 791$ yr (the periastron of the outer planet is retrograde).

Due to the mean motion resonance, a linear relation links the first four fundamental frequencies in Table 3: $n_{\rm b}- 5n_\sc + g_1 + 3g_2 = 0$ . As a consequence, one of them is superfluous. A new fundamental frequency associated to the resonance $l_\theta$ replaces it.

Table 3: Fundamental frequencies for S2.

The solution S2 is trapped in the 5:1 mean motion resonance with the following main resonant argument

$\begin{displaymath}% \theta=\lambda_{\rm b}- 5\lambda_\sc + \varpi_{\rm b}+ 3\varpi_\sc. \end{displaymath}$

(8)

The variations of $\theta$ versus time are plotted in Fig. 7 and exhibit a libration around $\theta_0 = 0\hbox{$^\circ$ }$ . We nonetheless observe that this libration is the modulation of several different terms with similar amplitudes of approximately $40\hbox{$^\circ$ }$ , but on different time scales. This leads to a libration with an amplitude that can be higher than $180\hbox{$^\circ$ }$ .

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig07} \end{figure}$

Figure 7:

Time variation of the resonant argument $\theta = \lambda_{\rm b}- 5\lambda_\sc + \varpi_{\rm b}+ 3\varpi_\sc$ for the orbital solution S2 (green line). $\theta$ is in libration around $\theta _0 = 0^\circ$ , with a modulation of two terms of period $P_{\theta } = 80.13$ years, and $P_{\Delta \varpi }=749.10$ years, amplitudes of about 35 and 50 degrees, respectively. The black line shows the first secular term contribution.

Open with DEXTER

In order to describe the behavior of $\theta$ more accurately, we searched for a quasi-periodic decomposition of $\theta(t)$ . We started with a frequency decomposition using frequency analysis (Laskar 2005) as

$\begin{displaymath}% \theta(t) = \sum_j{A_j\cos(\nu_j t + \phi_j)}. \end{displaymath}$

(9)

We then decomposed each frequency $\nu_j$ on the four fundamental frequencies ( $n_{\rm c}$ , g₁, g₂, $l_\theta$ )

$\begin{displaymath}% \nu_j = \alpha_j n_\sc + \beta_j g_1 + \delta_j g_2 + \gamma_j l_\theta, \end{displaymath}$

(10)

where $\alpha_j,\beta_j,\delta_j$ , and $\gamma_j$ are integers. Each term of the decomposition follows a D'Alembert-like relationship expressed by

$\begin{displaymath}% \alpha_j+\beta_j+\delta_j=0, \end{displaymath}$

(11)

which can be used to simplify the expression of $\nu_j$ . Indeed, by rearranging the righthand side of Eq. (10), we can write

$\begin{displaymath}% \nu_j = \alpha_j n_\sc + \beta_j (g_1 - g_2) + (\beta_j+\delta_j)g_2 + \gamma_j l_\theta, \end{displaymath}$

(12)

and using Eq. (11), we get

$\begin{displaymath}% \nu_j = \alpha_j (n_\sc - g_2) + \beta_j (g_1 - g_2) + \gamma_j l_\theta. \end{displaymath}$

(13)

Table 4: Quasi-periodic decomposition of the resonant angle $\theta = \lambda_{\rm b}- 5\lambda_\sc + \varpi_{\rm b}+ 3\varpi_\sc$ for an integration of the orbital solution S2 over 1 million years.

The first thirty terms of the decomposition can be found in Table 4. The decrease in the amplitudes A_j is slow owing the strong perturbations.

$\begin{figure} \par\includegraphics[width=18cm,clip]{figs/13635_fig08} \end{figure}$

Figure 8:

Mean motion diffusion of test particles. integration S2 solution with massless particles over 16 000 yr. We computed two determination n and $n^\prime$ of each particules mean motion over two consecutive time intervals of 8000 yr. The color scale codes the stability index $\log{\vert n-n^\prime\vert}$ . Initial conditions for the massless particles are $(i_{\rm p},M_{\rm p},\omega _{\rm p},\Omega _{\rm p}) = (90\hbox {$^\circ $ },0\hbox {$^\circ $ },0\hbox {$^\circ $ },0\hbox {$^\circ $ })$ . In both panels, we vary the initial eccentricity of the test particules from 0 to 0.8 with a step size of 0.004. In the left panel, the semi-major axis varies from 0.05 AU to 0.5 AU with a step size of 0.0025 AU, and in the right panel it varies from 0.5 AU to 10 AU with a step size of 0.05. The white crosses mark the position of the two planets, and the collision lines are traced with white lines.

Open with DEXTER

The first term (j=1) is responsible for the long-term oscillations with frequency g₁ - g₂. The period corresponding to g₁ - g₂, associated to the angle $\Delta\varpi=\varpi_{\rm b}- \varpi_\sc$ , is $P_{\Delta\varpi} \approx 749.09$ years. It is worth mentioning that the angle $\Delta\varpi$ is not in libration in this system, as opposed to what has been observed in planetary systems locked in a 1/2 mean motion resonance, such as GJ 876 (Laughlin & Chambers 2001; Ji et al. 2002; Correia et al. 2010; Lee & Peale 2002) or HD 82943 (Gozdziewski & Maciejewski 2001), or in a 3/2 mean motion resonance such as HD 45364 (Correia et al. 2009). The second one (j=2) introduces the short-term oscillations of frequency $l_\theta$ , corresponding to a period $P_{\theta} \approx 101.49$ years. More precisely, the libration of $\theta$ is made of three different kinds of contributions:

secular terms: $\alpha_j=\gamma_j=0$ , hence of the form $\beta_j(g_1-g_2)$ ;
resonant terms: $\alpha_j=0$ and $\gamma_j \ne 0$ , of the form $\beta_j(g_1-g_2)+\gamma_j g_\theta$ ;
and short-period terms: $\alpha_j \ne 0$ .

We plotted also in Fig. 7 the contribution of the term j=2, the first secular contribution.

4.4 Test particles

In this section, we test the possibilities for a third body around HD 202206. To that end we integrate the S2 solution with added massless particles over 16 000 years. As explained in the previous sections, we compute two determinations n and $n^\prime$ of each particle mean motion over two consecutive time intervals of 8000 years. We then obtain a stability index $D = \log{\vert n-n^\prime\vert}$ for each particle, which we plotted with a color code in Fig. 8. Both panels are (a,e) grids of initial conditions of the test particles. We vary the semi-major axis from 0.05 AU to 0.5 AU with a 0.0025 step size in the left panel, and from 0.5 AU to 10 AU with a 0.05 step size in the right panel. We span eccentricities of the test particles from 0 to 0.9 with a 0.005 step size. Since the particle mean motions are higher in the left panel, they were integrated with a time step of 10^-3 years instead of 2 $\times$ 10^-3.

Due to the very large eccentricities of HD 202206 b and HD 202206 c, the dynamical environment between and around them is very unstable. As a result, we do not expect any viable planets with a semi-major axis between approximately 0.12 AU and 6.5 AU. Most of these particles were actually lost before the end of the integration, either through collision or by having their eccentricity increased higher than 1. The same computation with particles of one Earth mass yields very similar results. Assuming that S2 is a good representation of the HD 202206 planetary system, we can use these results to put constraints on hypothetical and still undetected additional companions. There are clearly two possible regions for new planets: either close to the star (a<0.12 AU) or outside HD 202206 c (a>6.5 AU).

In the first case, any planet massive enough should have already been detected, since many full period are available in the data. Assuming a low eccentricity for the hypothetical companion and a 6 m/s instrumental precision, we find that planets bigger than 24 earth masses should have already been detected. A Neptune-sized planet can exist anywhere between 0.06 AU and 0.12 AU, and a 10 earth mass planet anywhere between 0.02 AU and 0.12 AU.

In the second case (a>6.5 AU), the period is longer than 16 years, meaning that we have only covered approximately half an orbit at best. However, a massive enough planet would create at least a detectable trend in the data. At 6.5 AU we can rule out the existence of a planet with more than half a Jupiter mass. At 10 AU we can rule out a planet between 1 ${M}_{\rm Jup}$ and 3 ${M}_{\rm Jup}$ , depending on the phase. We conclude that a yet undetected planet smaller than half a Jupiter mass could exist at semi-major axis greater than 6.5 AU.

4.5 Finding the center of libration inside a resonance through frequency analysis

4.5.1 Center of libration

We consider now the planar three-body problem and, more particularly, the planetary problem with a p:p+1 mean-motion resonance. The problem is to find the orbit center of libration starting from a quasi-periodic orbit in the resonance.

In the restricted case, this center of libration is a well-defined periodic orbit, such as the Lagrangian points in the 1:1 resonance. However, in the general problem it is not so easy to define properly this orbit and to find it. The system has 4 degrees of freedom, and its quasi-periodic orbits live on 4-torus in the phase space, although it can be restricted to 3 degrees of freedom using the angular momentum reduction. Each dimension of the 4-torus is associated to one of the four fundamental frequencies (f₀,f₁,f₂,f₃) of the orbit. Supposing that f₀ is the frequency of the resonant mode, the orbit equivalent to the center of libration is living on a 3-torus, depending only on (f₁,f₂,f₃). In other words, it is a torus where the fourth dimension associated to the resonance has zero amplitude. Intuitively, we can represent it as the center of the 4-torus in the dimension associated with the resonance.

The analysis of resonant orbits at the center of libration, or in its vicinity, has already been performed by several authors for the 2:1 resonant systems GJ 876 (Beaugé et al. 2003; Lee & Peale 2002) or for the system HD 82943 that as been also proposed as a 2:1 resonance candidate (Lee et al. 2006; Beaugé et al. 2008). But in both cases, the resonant orbit at the center of libration was assumed to be a periodic orbit, resulting from the double resonance in longitude and perihelions. Here, it should be noted that the configuration is more complex, as the perihelions are not supposed to be locked into resonance. The orbit at the center of libration is thus not periodic, but only quasiperiodic.

4.5.2 Quasi-periodic decomposition

If this orbit was periodic we could use a simple Newton algorithm to find it (provided that we start within the convergence radius), as it is a fixed point on a Poincaré map. This method has been extensively used in numerical searches of periodic orbit families of the three-body problem (see for instance Henon 1997,1974).

We present here a new numerical method of finding the quasi-periodic center of libration, using the fact that we can get an accurate quasi-periodic decomposition of a numerically integrated quasi-periodic orbit with frequency map analysis (Laskar 2005). Let ${\vec X}(t)$ be a state vector of such an orbit. Using frequency analysis, we can obtain a quasi-periodic representation for any component x of the vector $\vec{X}$

$\begin{displaymath}% x(t) = \sum_{(k)}{A_{(k)}{\rm E}^{{\rm i} \langle k,f \rangle t}} \end{displaymath}$

(14)

for each coordinate of the vector $\vec{X}$ , with $k = (k_0, k_1, k_2, k_3) \in \mathbb{Z}^4$ , f = (f₀, f₁, f₂, f₃), and $A_{(k)} \in \mathbb{C}$ . We can separate these sums in two parts: x(t) = u(t) + v(t) where u does not depend on the resonant frequency f₀, and v has all the terms depending on f₀.

$\begin{displaymath}% \left\{\begin{array}{ll} u(t) &= \sum_{(k), k_0 = 0}{A_{(k)... ...)}{\rm E}^{{\rm i} \langle k,f \rangle t}}. \end{array}\right. \end{displaymath}$

(15)

The quasi-periodic orbit described by ${\vec U}(t)$ is precisely lying on a 3-torus that has the characteristics we are looking for. However, it is probably not a solution of the equations of motion. But we can assume that it is close to one. Hence, we can use ${\vec U}(0)$ as a new initial condition, and obtain a new resonant quasi-periodic orbit ${\vec X}^\prime(t)$ with

$\begin{displaymath}% {\vec X}^\prime(0) = {\vec U}(0). \end{displaymath}$

(16)

The amplitude of the resonant terms in this new orbit will be smaller. In other words, it lives on 4-torus closer to the 3-torus of the center of libration. We can then iterate this procedure to suppress all the terms with f₀.

$\begin{figure} \par\includegraphics[width=6cm,clip]{figs/13635_fig09} \end{figure}$	Figure 9: Schematic representation of the 3-torus of the center of libration. We represent the 4-torus of a given resonant orbit by a 2-torus (a doughnut) as it is not possible to represent it otherwise. The center of libration 3-torus is then represented by the circle in the center of the interior of the doughnut.
Open with DEXTER

4.5.3 Application to HD 202206

We present here its application to the HD 202206 system, using our orbital solution S2 as a starting resonant orbit. We work with the $\vec{X} = (a_{\rm b}, \lambda_{\rm b}, z_{\rm b}, a_\sc,\lambda_\sc, z_\sc)$ state vector, where $z_{\rm b}= e_{\rm b}{\rm E}^{{\rm i}\varpi_{\rm b}}$ and $z_\sc = e_\sc {\rm E}^{{\rm i}\varpi_\sc}$ . Each step p of the method is decomposed as follows:

numerical integration of ${\vec X}^{(p)}$ ;
determination of (f₀,f₁,f₂,f₃) using frequency analysis;
quasi-periodic decomposition: $x^{(p)}(t) \approx u^{(p)}(t) + v^{(p)}(t)$ ;
new initial conditions: x^(p+1)(0) = x^(p)(0) - v^(p)(0).

The initial conditions are not computed exactly following Eq. (16). Since we work with a finite number of terms and the amplitudes of the terms in v(t) are supposed to become small, the error due to the limited number of terms will decrease.

The convergence proved to be fast as we reduced the amplitude of the most important resonant terms by 2 orders of magnitude in just 4 steps (Fig. 11). We show graphically the decreasing amplitude of libration at each step in Fig. 10 where we plot approximated sections of the successive torus projected in the $(a_{\rm c},\varpi_{\rm c})$ plane. No resonant terms are found in the first 100 terms of the quasi-periodic decomposition of each variable at the last step (with the exception of $\lambda_{\rm b}$ ). In fact, in half of the variables, there are no resonant terms left in the first 300 terms (Fig. 11).

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig10.eps} \end{figure}$	Figure 10: Projection of a section of the 4-torus of each step's trajectory in the $(a_\sc,\omega_\sc)$ plane. Each torus is approximated using a (truncated) quasi-periodic decomposition of the trajectory.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig11} \end{figure}$

Figure 11:

Evolution of the amplitude of the libration mode in the quasi-periodic decomposition of each variable at each step. In the top panel we plot the relative amplitude of the first term depending on f₀ compared to the first non-constant term. In the bottom panel we plot the position of this first term in the decomposition. The first step (abscissa 0) is the S2 orbital solution, and the last step (abscissa 6) is the orbital solution S3 (Table 5).

Open with DEXTER

Once we find an orbit with a zero amplitude libration mode to a good approximation, we also get the 3-torus it is living on since its quasi-periodic decomposition gives us a parametrization of the torus. The three angular variables of the torus appears in the decomposition

$\begin{displaymath}% \left\{\begin{array}{ll} \phi_1 &= f_1 t + \psi_1 \\ \ph... ... t + \psi_2 \\ \phi_3 &= f_3 t + \psi_3 \end{array} \right. \end{displaymath}$

(17)

where $\psi_1$ , $\psi_2$ , and $\psi_3$ are initial phases. With $\xi_k = {\rm E}^{i\phi_k}$ , we can rewrite Eq. (14) to reveal the underlying torus

$\begin{displaymath}% x(\phi_1,\phi_2,\phi_3) = \sum_{(k)}{A_{(k)} \xi_1^{\alpha_{(k)}} \xi_2^{\beta_{(k)}} \xi_3^{\gamma_{(k)}}}. \end{displaymath}$

(18)

Of all the orbits living on this torus we can choose the closest to the radial velocity data. This can be done by minimizing the $\chi ^2$ in the three initial phases $\psi_1$ , $\psi_2$ , and $\psi_3$ .

We denote S3 as the orbit obtained that way at step 6. We give initial conditions for S3 in Table 5. This solution yields a square root of reduced $\chi ^2$ equal to 1.55. If the system is locked in the 5:1 mean motion resonance, the libration mode is likely to be dampened through dissipative processes. In that regard, solution S2 is unlikely as it is on the edge of the resonant island, and exhibits a high resonant mode amplitude ( $40\hbox{$^\circ$ }$ for the resonant critical angle). We expect that the real solution will be closer to S3.

Table 5: Orbital parameters of an orbit close to the center of libration of the 5:1 mean motion resonance.

5 Coplanar orbits

5.1 Stability and low $\chi ^2$ orbits

In this section we investigate the system behavior when both planets remain in the same orbital plane ( $I = 0\hbox{$^\circ$ }$ ) and are prograde, but the inclination to the line of sight is lower than $90\hbox {$^\circ $ }$ :

$i_{\rm b}=i_\sc=i \leq 90\hbox{$^\circ$ }$ ;
$\Omega_\sc = 0\hbox{$^\circ$ }$ .

For each inclination value, we compute a new best fit with a Levenberg-Marquard algorithm. The $\sqrt {\chi _{\rm r}^2}$ value obtained at different inclinations is plotted in Fig. 13. Interestingly enough, we obtain better fits for lower inclination down to $15\hbox{$^\circ$ }$ .

In this configuration, the masses of the two companions grow approximately in proportion with $1/\sin{i}$ when the inclination i diminishes. There are thus only small changes between $90\hbox {$^\circ $ }$ and $50\hbox{$^\circ$ }$ , as seen in Fig. 12 (two top leftmost panels). For lower inclinations, because of the increased masses, mutual perturbations become stronger and less orbits are stable. However, among the lowest $\sqrt {\chi _{\rm r}^2}$ orbits, the ones in the 5:1 mean motion resonance remain stable for inclinations up to $15\hbox{$^\circ$ }$ . For $i=10\hbox{$^\circ$ }$ , no stable orbits are left (bottom rightmost panel). This puts clear limits on the inclination of the system and on the masses of the two companions: $m_{\rm b}< 95.5~M_{\rm jup}$ and $m_\sc < 12.5~M_{\rm jup}$ .

$\begin{figure} \par\includegraphics[width=17cm,clip]{figs/13635_fig12} \vspace*{5mm} \end{figure}$	Figure 12: Dynamics of a coplanar HD 202206 system for different values of the inclination i. Each panel is a diffusion map in the $(a_\sc,e_\sc)$ plane of initial conditions constructed the same way as the top panel of Fig. 5 ( the top-left panel is in fact a copy). For each value of inclination i, a new fit that assumes a coplanar system is computed.
Open with DEXTER

An interesting property shown in Fig. 12 is that the lowest $\chi ^2$ orbits are always in the vicinity of the stable island of the 5:1 mean motion resonance for all inclinations. We believe that it is a strong point supporting the hypothesis that the HD 202206 system is locked in this resonance. It appears, however, that after $i=50\hbox {$^\circ $ }$ , and for lower inclinations, the low $\chi ^2$ region is slightly shifted towards $a_\sc$ values that are lower than the resonance island location.

To have a more precise picture, we spanned i values from $10\hbox {$^\circ $ }$ to $90\hbox {$^\circ $ }$ with a step size of $5\hbox{$^\circ$ }$ . For each inclination value, we started by computing a new best fit with the Levenberg-Marquard algorithm. The $\sqrt {\chi _{\rm r}^2}$ values and stability indexes D₁ and D₂ are computed in the $(a_\sc,e_\sc)$ plane of initial conditions, around the best fit. The size of each grid is 101 $\times$ 161 dots, with step sizes of 0.0025 AU and 0.002. For each of these maps, we computed $D_{\rm lim}$ (see Sect. 3.4), and detected stable regions and their relative positions to the observations in terms of $\sqrt {\chi _{\rm r}^2}$ .

To get a synthetic vision of all this data, we can get an estimation of the lowest stable $\sqrt {\chi _{\rm r}^2}$ for each inclination (Fig. 13). We plotthe lowest $\sqrt {\chi _{\rm r}^2}$ of $D^{\rm s}_1 < D_{\rm lim}$ orbits with respect to i (broken curve) along with the best fit value (solid curve). It mainly confirms that the distance between stable orbits and low $\sqrt {\chi _{\rm r}^2}$ orbits is really small from $90\hbox {$^\circ $ }$ to $50\hbox{$^\circ$ }$ . Is is actually slightly better between $60\hbox{$^\circ$ }$ and $50\hbox{$^\circ$ }$ .

To have a better idea of the correlation between low $\sqrt {\chi _{\rm r}^2}$ orbits and stable orbits, we looked at the percentage of stable orbits inside a given $\sqrt {\chi _{\rm r}^2}$ level curves (Fig. 14). It gives a synthetic representation of the overlap between low $\sqrt {\chi _{\rm r}^2}$ regions and stable regions. It is very clear that inclinations lower than $30\hbox{$^\circ$ }$ are very unlikely, although possible. It also appears that inclinations between $90\hbox {$^\circ $ }$ and $50\hbox{$^\circ$ }$ are the most probable.

If we limit the acceptable inclinations to $30\hbox{$^\circ$ }\le i \le 90\hbox{$^\circ$ }$ , we can derive limits for the masses of the inner and outer planets (assuming a coplanar prograde configuration):

$1 \le m_{\rm b} / m_{\rm b}^{(0)} \le 2$ ;
$1 \le m_{\rm c} / m_{\rm c}^{(0)} \le 2$ .

It is interesting to estimate the time when the radial velocity data will allow a more accurate estimation of the inclination i and masses for the system. Assuming this model is close enough to the true system, we look at the differences of the radial velocities in the $i=90\hbox {$^\circ $ }$ case, and $i=50\hbox {$^\circ $ }$ or $i=30\hbox {$^\circ $ }$ . With the hypothesis of an instrumental precision of 7 m/s, we find (Fig. 15) that we will have to wait until about 2015 to discriminate between $i=90\hbox {$^\circ $ }$ and $i=30\hbox {$^\circ $ }$ , which corresponds to approximately a factor 2 in the masses. For $i=50\hbox {$^\circ $ }$ , five more years are needed for the differences between the two configurations to be greater than the measurements precision.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig13} \end{figure}$

Figure 13:

Evolution of the best fit $\sqrt {\chi _{\rm r}^2}$ (solid curve) and an estimation of the best stable orbits $\sqrt {\chi _{\rm r}^2}$ as a function of the inclination i. For each value of i between $90\hbox {$^\circ $ }$ and $10\hbox {$^\circ $ }$ , a new fit is computed, assuming a coplanar configuration. We then look for stable orbits around each of these orbital solutions in a $(a_\sc,e_\sc)$ plane of initial conditions.

Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig14} \end{figure}$	Figure 14: Percentage of stable orbits inside a given $\sqrt {\chi _{\rm r}^2}$ level curve against i (see Fig. 13).
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig15} \end{figure}$

Figure 15:

Differences between the radial velocity of a stable solution with $i=90\hbox {$^\circ $ }$ , and $i=50\hbox {$^\circ $ }$ (red curve), or $i=30\hbox {$^\circ $ }$ (green curve). The starting epoch of the three integrations is $2~453~000.00~{\rm JD}$ . The blue horizontal lines represents the precision of the current set of data, obtained with the CORALIE instrument. The vertical blue line marks the beginning of 2009.

Open with DEXTER

5.2 Resonant and secular behavior

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig16} \end{figure}$	Figure 16: Libration and secular frequencies. For each inclination i, we pick a stable orbit with a low $\chi ^2$ and integrate it over 1 million years. The fundamental frequencies are then determined through frequency analysis. The solid curves represent, for each frequency, a $1/\sin^2{i}$ law.
Open with DEXTER

We end this study of the coplanar configurations with a quick look at the dependence of the resonant and secular dynamics on the inclination. For each inclination, we picked a stable orbit with a low $\sqrt {\chi _{\rm r}^2}$ value, and plotted its fundamental frequencies $l_\theta$ , g₁, and g₂ in Fig. 16. As expected from the perturbation theory, when the masses increase, the secular frequencies also increase (in absolute value). We verify that it follows a rule in $1/\sin^2{i}$ . This is a consequence of the fact that the most important terms responsible for the secular dynamics are of order two with respect to the masses.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig17} \end{figure}$

Figure 17:

Stable configurations for $i_{\rm b} = 90\hbox {$^\circ $ }$ (i.e. for the minimum mass of planet b). For each value of $i_{\rm c}$ in $(90\hbox {$^\circ $ }, 30\hbox {$^\circ $ }, 11\hbox {$^\circ $ })$ , and $\Delta \Omega$ between $0\hbox {$^\circ $ }$ and $360\hbox {$^\circ $ }$ we compute the diffusion index over grids in the $(a_{\rm c},e_{\rm c})$ plane of initial conditions with step sizes of 0.0025 AU and 0.004, respectively. Each grid is also centered on the minimum $\chi ^2$ computed for each pair $(i_{\rm c}, \Delta \Omega )$ . We plot against $\Delta \Omega$ the square root of the minimum $\chi ^2$ in the top panel, and the proportion of stable orbits inside $\sqrt {\chi _{\rm r}^2}= 1.5$ in the middle panel. The mutual inclination corresponding to each triplet $(i_{\rm b}, i_{\rm c}, \Delta \Omega )$ is plotted in the bottom panel.

Open with DEXTER

6 Mutually inclined orbits

In this section we drop the coplanarity constraint. We allow the inclinations $i_{\rm b}$ and $i_\sc$ to vary independently and allow variations in $\Omega_\sc$ , the longitude of the ascending node of c. To span the possible values for $i_{\rm b}$ , $i_\sc$ , and $\Omega_\sc$ in an efficient manner, we restrict ourselves to two mass ratios $k_{\rm b}$ (Eq. (3)) for planet b

$k_{\rm b} = 1$ ( $i_{\rm b} = 90\hbox {$^\circ $ }$ );
$k_{\rm b} = 2$ ( $i_{\rm b} = 30\hbox {$^\circ $ }$ );

and three mass ratios $k_{\rm c}$ for planet c

$k_{\rm c} = 1$ ( $i_{\rm c} = 90\hbox{$^\circ$ }$ );
$k_{\rm c} = 2$ ( $i_{\rm c} = 30\hbox{$^\circ$ }$ or $i_{\rm c} = 150\hbox{$^\circ$ }$ );
$k_{\rm c} \simeq 5$ ( $i_{\rm c} = 11\hbox{$^\circ$ }$ or $i_{\rm c} = 169\hbox{$^\circ$ }$ ).

For each pair of inclinations $i_{\rm b}$ and $i_\sc$ , we let $\Delta \Omega$ vary between $0\hbox {$^\circ $ }$ and $360\hbox {$^\circ $ }$ with a $10\hbox {$^\circ $ }$ step size, and for each triplet $(i_{\rm b},i_\sc,\Delta\Omega)$ , we perform a fit with the Levenberg-Marquard minimization and compute a diffusion grid in the $(a_\sc,e_\sc)$ plane of initial conditions around the fit. The step sizes of the grids are respectively 0.0025 AU and 0.004. We plotted in Figs. 17c and 18c the mutual inclination I as a function of $\Delta \Omega$ for reference. For each configuration we look at the proportion of stable orbits inside the $\sqrt {\chi _{\rm r}^2}= 1.5$ level curve (Figs. 17b and 18b). We assume that stable zones with $\sqrt{\chi_{\rm r}^2}\le 1.5$ orbits harbor potential solutions of the system. We keep the initial inclination $i_{\rm b}$ in $]0;\pi/2]$ as opposed to $i_{\rm c}$ (see Sect. 3.1). Also for any given values $\alpha \in ]0;\pi/2]$ and $\beta \in [0;2\pi[$ , assuming all the other orbital elements identical, the configuration ( $i_{\rm b}=\pi/2$ , $i_{\rm c}=\alpha$ , $\Omega_{\rm c}=\beta$ ) is symmetric to the configuration ( $i_{\rm b}=\pi/2$ , $i_{\rm c}=\pi-\alpha$ , $\Omega_{\rm c}=-\beta$ ) with respect to the plane $(\vec{i},\vec{k})$ (see Fig. 2). Since this plane contains the line of sight, the two configurations are indistinguishable using radial velocity measurements. Hence for $k_{\rm b} = 1$ ( $i_{\rm b} = 90\hbox {$^\circ $ }$ ), we do not look at $i_{\rm c} = 150\hbox{$^\circ$ }$ and $169\hbox{$^\circ$ }$ .

6.1 k_b = 1

For $k_{\rm b} = 1$ ( $i_{\rm b} = 90\hbox {$^\circ $ }$ , Fig. 17) we find that significant stable zones inside $\chi^2 = 1.5$ exist mostly for aligned and anti-aligned ascending nodes (i.e. $\Delta\Omega \simeq 0\hbox{$^\circ$ }$ and $180\hbox{$^\circ$ }$ in Fig. 17b). For $k_{\rm c} = 1$ (red curves) the aligned configurations ( $\Delta \Omega =0\hbox {$^\circ $ }$ $\pm$ $10\hbox {$^\circ $ }$ ) are coplanar and prograde ( $I \simeq 0\hbox{$^\circ$ }$ ), and the anti-aligned configurations ( $\Delta\Omega = 180\hbox{$^\circ$ }$ $\pm$ $10\hbox {$^\circ $ }$ ) are coplanar and retrograde. Outside those two particular cases, we find no significant stable zones for $\sqrt{\chi_{\rm r}^2}\le 1.5$ . Indeed, while the resonant island is roughly centered on the lowest $\chi ^2$ level curve, it is not stable outside the coplanar configurations. We find that an extended zone of stability exists outside the resonance for $\Delta\Omega \simeq 30\hbox{$^\circ$ }$ , but it lies just outside the 1.5 level curve. Due to symmetries, the situation for $\Delta\Omega > 180\hbox{$^\circ$ }$ mirrors that of $\Delta\Omega < 180\hbox{$^\circ$ }$ .

For $k_{\rm c} = 2$ (green curves) we again find stable zones for $\Delta\Omega \simeq 0\hbox{$^\circ$ }$ , but not around $180\hbox{$^\circ$ }$ . However stable regions with potential solutions exist for $\Delta \Omega$ values up to $90\hbox {$^\circ $ }$ , corresponding to mutual inclinations between $60\hbox{$^\circ$ }$ and $90\hbox {$^\circ $ }$ . This is mostly due to the minimum $\chi ^2$ getting smaller up to $\Delta\Omega = 90\hbox{$^\circ$ }$ (see the green curve in Fig. 17a). While the stable regions, both in the resonance and outside, shrink when $\Delta \Omega$ augments, the $\chi^2 = 1.5$ level curves encompass a larger area. Finally, for $k_{\rm c} = 5$ (blue curves), no significant stable zones are found at low $\chi ^2$ values.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig18}\vspace*{-0.25mm} \end{figure}$	Figure 18: Stable configurations for $i_{\rm b} = 30\hbox {$^\circ $ }$ (i.e. for approximately twice the minimum mass: $m_{\rm b} \simeq 33.5~M_{{\rm jup}}$ ). See Fig. 17 caption for more details.
Open with DEXTER

To summarize, potential solutions (stable orbit with a low $\chi ^2$ ) for $i_{\rm b} = 90\hbox {$^\circ $ }$ mainly exists for coplanar configurations, inside the 5:1 mean motion resonance. If planet c's mass is kept low ( $k_{\rm b} < 5$ ), stable regions do exist for non coplanar configurations, but they are located outside the lowest $\chi ^2$ region. A noteworthy exception is the retrograde configuration, where stable orbits are only found close to the coplanar case, which only happens for $k_{\rm c} = 1$ and $\Delta\Omega \simeq 180\hbox{$^\circ$ }$ .

6.2 k_b = 2

When we double the mass of planet b ( $i_{\rm b} = 30\hbox {$^\circ $ }$ , Fig. 18), once again retrograde potential solutions can only occur for coplanar orbits: $i_{\rm c} = 150\hbox{$^\circ$ }$ and $\Delta\Omega \simeq 180\hbox{$^\circ$ }$ (green dotted curve). Concerning prograde orbits, potential stable solutions exist for a mutual inclination of $I \simeq 60\hbox{$^\circ$ }$ :

$i_{\rm c} = 90\hbox{$^\circ$ }$ and $\Delta \Omega =0\hbox {$^\circ $ }$ (red curve);
$i_{\rm c} = 30\hbox{$^\circ$ }$ and $\Delta\Omega = 180\hbox{$^\circ$ }$ (green solid curve).

6.3 Conclusions

To sum up, we find that the configurations that have a significant stable zone at low $\chi ^2$ values are mostly found when the two lines of nodes are aligned. That is for $\Delta \Omega$ close to $0\hbox {$^\circ $ }$ or $180\hbox{$^\circ$ }$ . In addition, stable orbits with the lowest $\chi ^2$ are all in the 5:1 mean motion resonance except for nearly coplanar retrograde configurations, where they can also be close to the commensurability. However, we can also find stable orbits outside the resonance in the prograde resonance, usually with $\sqrt {\chi _{\rm r}^2}$ higher than 1.45. Retrograde configurations seem to be limited to nearly coplanar orbits with anti-aligned ascending nodes. Other than that, we could not find any clear correlation with mutual inclination.

That we find non-resonant stable solutions for retrograde configurations is consistent with Smith & Lissauer (2009). They show that retrograde configurations allow more closely packed systems than prograde configurations. It is also suggested by Gayon & Bois (2008) that retrograde configurations are likely alternatives both from the radial velocity data and the long-term stability point of view. This hypothesis has been reinforced by the recent detection of WASP-17b, an ultra-low density planet in a proposed retrograde orbit (Anderson et al. 2010). Forming such a system however does remain difficult, so we do not favor this hypothesis.

7 Discussion and conclusion

By assuming that the system is coplanar, we performed a systematic study of the dynamics of the system for different inclinations to the line of sight. We were able to find constraints for the inclination to the line of sight: $30\hbox{$^\circ$ }\le i \le 90\hbox{$^\circ$ }$ . This means that the companions' masses are most likely not greater than twice their minimum values:

$1 \le m_{\rm b} / m_{\rm b}^{(0)} \le 2$ ;
$1 \le m_{\rm c} / m_{\rm c}^{(0)} \le 2$ .

We also studied the influence of the planets mutual inclination for two different inclination values of planet b ( $i_{\rm b} = 90\hbox {$^\circ $ }$ and $i_{\rm b} = 30\hbox {$^\circ $ }$ ), but did not find any clear correlation other than that a retrograde potential stable solution consistent with the radial velocity data seems to be limited to mutual inclinations close to $180\hbox{$^\circ$ }$ (i.e. nearly coplanar orbits). Like Gozdziewski et al. (2006), we find possible stable solutions with low $\chi ^2$ for a wide range of mutual inclinations between $0\hbox {$^\circ $ }$ and $90\hbox {$^\circ $ }$ . The current data cannot yield more precise constraints. Also the masses determination is dependent on the stellar mass which is not well established.

Although all published dynamical studies of HD 202206 suggest that b and c are in a 5:1 mean motion resonance, it is still a debated question. For instance, Libert & Henrard (2007) assume that it is just close to the commensurability. Libration of $\Delta\varpi$ occurs for particular initial values of this angle, providing a stabilizing mechanism outside the mean-motion resonance, not far from the best fits. In most cases other than retrograde coplanar configurations, those orbits in near-commensurability are worse solutions than the ones in resonance, but they could be more probable if the eccentricities are overestimated (especially for b). We find that all significant stable zones with the best O-C are in the 5:1 mean motion resonance. In fact, the minimum $\chi ^2$ is almost always in the resonance or very close to it, and stable orbits in the resonance can be found with $\chi ^2$ not significantly higher than the best fit. In addition the O-C level curves tend to follow the resonant island roughly, even though the agreement is not as perfect as for the HD 45364 system (Correia et al. 2009). This is an improvement over Correia et al. (2005), where the best fit lay outside the resonant island and the $\chi ^2$ had to be degraded to find a stable solution. We thus believe that the resonant configuration is the most probable one. We provide a stable solution (S2, Table 2) in the coplanar edge-on case. This solution shows a high amplitude resonant mode in the libration of the critical angle. We believe that this resonant mode is probably dampened by dissipative processes. We used frequency analysis to find a tore on which such orbits exist. Although the specific orbit we give in Table 5 does not have a very low $\sqrt {\chi _{\rm r}^2}$ at 1.55, we expect that the true orbit will be close to it with a low libration amplitude.

For retrograde configurations, the picture is quite different. The best fit lies in a very stable region just outside the mean motion resonance. While these orbits are valid candidates from the dynamical and the observational points of view, we do not favor them because the formation of these systems is hard to explain.

We investigated the possibility of undetected companions. We found that planets with masses lower than approximately one Neptune mass can exist for a semi-major axis lower than 0.12 AU, and that $0.5~M_{\rm Jup}$ planets are also possible beyond 6.5 AU. No planets are possible between 0.12 AU and 6.5 AU as they would be unstable. The two planets model may prove to be wrong in the future, but these hypothetical new companions should not have a big impact on the already detected ones.

Acknowledgements

We acknowledge support from the Swiss National Research Found (FNRS), French PNP-CNRS, Fundação para a Ciência e a Tecnologia, Portugal (PTDC/CTE-AST/098528/2008), and Genci/CINES.

References

Anderson, D. R., Hellier, C., Gillon, M., et al. 2010, ApJ, 709, 159 [NASA ADS] [CrossRef] [Google Scholar]
Beaugé, C., Ferraz-Mello, S., & Michtchenko, T. A. 2003, ApJ, 593, 1124 [Google Scholar]
Beaugé, C., Giuppone, C. A., Ferraz-Mello, S., & Michtchenko, T. A. 2008, MNRAS, 385, 2151 [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. 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]
Fienga, A., Manche, H., Laskar, J., & Gastineau, M. 2008, A&A, 477, 315 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gayon, J., & Bois, E. 2008, A&A, 482, 665 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gozdziewski, K., & Maciejewski, A. J. 2001, ApJ, 563, L81 [NASA ADS] [CrossRef] [Google Scholar]
Gozdziewski, K., Konacki, M., & Maciejewski, A. J. 2006, ApJ, 645, 688 [NASA ADS] [CrossRef] [Google Scholar]
Henon, M. 1974, Celest. Mech., 10, 375 [NASA ADS] [CrossRef] [Google Scholar]
Henon, M. 1997, Generating Families in the Restricted Three-Body Problem, ed. M. Henon [Google Scholar]
Ji, J., Li, G., & Liu, L. 2002, ApJ, 572, 1041 [NASA ADS] [CrossRef] [Google Scholar]
Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
Laskar, J. 1993, Phys. D, 67, 257 [Google Scholar]
Laskar, J. 1999, in NATO ASI Hamiltonian Systems with Three or more Degrees of Freedom, ed. C. Simo (Kluwer), 134 [Google Scholar]
Laskar, J. 2005, in Hamiltonian systems and Fourier analysis, ed. D. Benest, C. Froeschle, & E. Lega (Cambridge Scientific Publishers), 99 [Google Scholar]
Laskar, J., & Robutel, P. 2001, Celest. Mech. Dynam. Astron., 80, 39 [NASA ADS] [CrossRef] [Google Scholar]
Laughlin, G., & Chambers, J. E. 2001, ApJ, 551, L109 [NASA ADS] [CrossRef] [Google Scholar]
Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 [NASA ADS] [CrossRef] [Google Scholar]
Lee, M. H., Butler, R. P., Fischer, D. A., Marcy, G. W., & Vogt, S. S. 2006, ApJ, 641, 1178 [NASA ADS] [CrossRef] [Google Scholar]
Lestrade, J.-F., & Bretagnon, P. 1982, A&A, 105, 42 [NASA ADS] [Google Scholar]
Libert, A.-S., & Henrard, J. 2007, A&A, 461, 759 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN, The art of scientific computing, 2nd edn. (Cambridge: University Press) [Google Scholar]
Smith, A. W., & Lissauer, J. J. 2009, Icarus, 201, 381 [NASA ADS] [CrossRef] [Google Scholar]
Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
Udry, S., Mayor, M., Naef, D., et al. 2002, A&A, 390, 267 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

Footnotes

... orbits: The CORALIE radial velocity measurements discussed in this paper are only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/519/A10
... measurements: CORALIE data were not taken during 2007-2008. There has nevertheless been a new reduction of the data since (Correia et al. 2005) and some previous observations were excluded here. The new data set corresponds to 14 additional measurements.
... system: The chosen value $\sqrt {\chi _{\rm r}^2}= 1.5$ corresponds to ${\approx}2\sigma$ on the variation in the parameters of the best-fit solution S1.

All Tables

Table 1: Best Newtonian fit S1 for the HD 202206 system assuming $i_{\rm b}=i_\sc=90\hbox{$^\circ$ }$ and $\Delta \Omega =0\hbox {$^\circ $ }$ .

Table 2: Stable orbital parameters S2 for the HD 202206 system for $i_{\rm b}=i_\sc=90\hbox{$^\circ$ }$ and $\Delta \Omega =0\hbox {$^\circ $ }$ .

Table 3: Fundamental frequencies for S2.

Table 4: Quasi-periodic decomposition of the resonant angle $\theta = \lambda_{\rm b}- 5\lambda_\sc + \varpi_{\rm b}+ 3\varpi_\sc$ for an integration of the orbital solution S2 over 1 million years.

Table 5: Orbital parameters of an orbit close to the center of libration of the 5:1 mean motion resonance.

All Figures

$\begin{figure} \par\includegraphics[angle=270,width=8.8cm,clip]{figs/13635_fig01} \end{figure}$	Figure 1: CORALIE radial velocities for HD 202206 superimposed on a 3-body Newtonian orbital solution (Table 1).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig02} \vspace*{5mm} \end{figure}$	Figure 2: Angles defining the orbit's orientation in space. We follow the observers' convention that sets the plane of sky as the reference frame, for which the edge-on coplanar configuration is $i_{\rm b}=i_\sc=90\hbox{$^\circ$ }$ and $\Delta \Omega =0\hbox {$^\circ $ }$ . The observer is limited to the velocity projected on the k axis.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig03} \end{figure}$	Figure 3: Distribution of D₁ < D₂ trajectories from the top panel of Fig. 6. Each integrated trajectory is binned with respect to its diffusion index $\log{D_1}$ , after the diffusion grid has been smoothed. For each bin, we computed the proportion of trajectories with a decreasing diffusion index over time: D₁ < D₂. The histogram shows the results for 0.5 wide bins. As expected from Eq. (7), orbits with a high diffusion index D₁, and D₁ < D₂ are nearly non existent ( $\log{D_1} > -1$ ), while we observe the opposite situation for low-diffusion index orbits ( $\log{D_1} < -3$ ).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig04} \end{figure}$	Figure 4: Percentage of orbits wrongly flagged asstable (false stable) or unstable (false unstable). Once $D_{\rm lim}$ is chosen using a given percentage threshold (see Fig. 3), we have a criterion for stable ( $D_1 < D_{\rm lim}$ ) and unstable orbits ( $D_1 > D_{\rm lim}$ ). We compared the results with a diffusion grid computed over a longer time interval (in this case 2 $\times$ 40 000 years) taken as a reference. A few orbits deemed stable from the reference diffusion grid were thought to be unstable and vice versa. The solid line traces the number of those faulty orbits for different values of the threshold. It appears to be lowest between 0.95 and 1. The dotted line traces the false positives, and the dashed line the false negatives. The solid one is the sum of those two.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig05} \end{figure}$	Figure 5: Global view of the dynamics of HD 202206 for variations of the semi-major axis and periastrum of the outer planet ( bottom panel) or semi-major axis and eccentricity ( top panel). The step sizes for $a_\sc$ , $\omega_\sc$ , $e_\sc$ are $0.005~{\rm AU}$ , $2\hbox {$^\circ $ }$ , and 0.004, respectively. The other parameters were kept constant and taken from the fit S1 (Table 1). The color scale is the stability index $\log{D_1}$ obtained through a frequency analysis of the longitude of the outer planet over two consecutive time intervals of 8000 years. The level curves give the $\sqrt {\chi _{\rm r}^2}$ value computed for each choice of parameters. The two horizontal black lines mark the intersection of the two grids. Most of the orbits are chaotic (yellow to red dots), but regular orbit zones exist (blue dots). One of these regular orbits regions lies inside the low $\sqrt {\chi _{\rm r}^2}$ region: the dark area around $a_\sc=2.5$ AU and $\omega_\sc=50\hbox{$^\circ$ }$ . It corresponds to the stable island of the 1/5 mean motion resonance. The cross marks the S1 orbital solution.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig06} \end{figure}$	Figure 6: Global view of the dynamics of HD 202206 for variations of the semi-major axis and periastrum ( bottom), and semi-major axis and eccentricity ( top) of the outer planet. The step sizes are respectively 0.0025 AU, $0.5\hbox {$^\circ $ }$ , and 0.002 for the eccentricity. The color scale is the stability index ( $\log{D_1}$ ), and the level curves give the $\sqrt {\chi _{\rm r}^2}$ values (see Fig. 5). The cross marks the best fit S1, and the horizontal line in each panel, the orbits common to both maps (i.e. the intersection of each map). The best fit lies very close to a stable resonant island where we can pick a stable solution with a $\chi ^2$ only marginally higher than that of the best fit. As an example, we picked such a solution (S2), marked by a white filled circle, by slightly increasing $a_{\rm c}$ from 2.4832 AU to 2.49 AU. The complete set of orbital elements for S2 is given in Table 2.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig07} \end{figure}$	Figure 7: Time variation of the resonant argument $\theta = \lambda_{\rm b}- 5\lambda_\sc + \varpi_{\rm b}+ 3\varpi_\sc$ for the orbital solution S2 (green line). $\theta$ is in libration around $\theta _0 = 0^\circ$ , with a modulation of two terms of period $P_{\theta } = 80.13$ years, and $P_{\Delta \varpi }=749.10$ years, amplitudes of about 35 and 50 degrees, respectively. The black line shows the first secular term contribution.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=18cm,clip]{figs/13635_fig08} \end{figure}$	Figure 8: Mean motion diffusion of test particles. integration S2 solution with massless particles over 16 000 yr. We computed two determination n and $n^\prime$ of each particules mean motion over two consecutive time intervals of 8000 yr. The color scale codes the stability index $\log{\vert n-n^\prime\vert}$ . Initial conditions for the massless particles are $(i_{\rm p},M_{\rm p},\omega _{\rm p},\Omega _{\rm p}) = (90\hbox {$^\circ $ },0\hbox {$^\circ $ },0\hbox {$^\circ $ },0\hbox {$^\circ $ })$ . In both panels, we vary the initial eccentricity of the test particules from 0 to 0.8 with a step size of 0.004. In the left panel, the semi-major axis varies from 0.05 AU to 0.5 AU with a step size of 0.0025 AU, and in the right panel it varies from 0.5 AU to 10 AU with a step size of 0.05. The white crosses mark the position of the two planets, and the collision lines are traced with white lines.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=6cm,clip]{figs/13635_fig09} \end{figure}$	Figure 9: Schematic representation of the 3-torus of the center of libration. We represent the 4-torus of a given resonant orbit by a 2-torus (a doughnut) as it is not possible to represent it otherwise. The center of libration 3-torus is then represented by the circle in the center of the interior of the doughnut.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig10.eps} \end{figure}$	Figure 10: Projection of a section of the 4-torus of each step's trajectory in the $(a_\sc,\omega_\sc)$ plane. Each torus is approximated using a (truncated) quasi-periodic decomposition of the trajectory.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig11} \end{figure}$	Figure 11: Evolution of the amplitude of the libration mode in the quasi-periodic decomposition of each variable at each step. In the top panel we plot the relative amplitude of the first term depending on f₀ compared to the first non-constant term. In the bottom panel we plot the position of this first term in the decomposition. The first step (abscissa 0) is the S2 orbital solution, and the last step (abscissa 6) is the orbital solution S3 (Table 5).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=17cm,clip]{figs/13635_fig12} \vspace*{5mm} \end{figure}$	Figure 12: Dynamics of a coplanar HD 202206 system for different values of the inclination i. Each panel is a diffusion map in the $(a_\sc,e_\sc)$ plane of initial conditions constructed the same way as the top panel of Fig. 5 ( the top-left panel is in fact a copy). For each value of inclination i, a new fit that assumes a coplanar system is computed.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig13} \end{figure}$	Figure 13: Evolution of the best fit $\sqrt {\chi _{\rm r}^2}$ (solid curve) and an estimation of the best stable orbits $\sqrt {\chi _{\rm r}^2}$ as a function of the inclination i. For each value of i between $90\hbox {$^\circ $ }$ and $10\hbox {$^\circ $ }$ , a new fit is computed, assuming a coplanar configuration. We then look for stable orbits around each of these orbital solutions in a $(a_\sc,e_\sc)$ plane of initial conditions.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig14} \end{figure}$	Figure 14: Percentage of stable orbits inside a given $\sqrt {\chi _{\rm r}^2}$ level curve against i (see Fig. 13).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig15} \end{figure}$	Figure 15: Differences between the radial velocity of a stable solution with $i=90\hbox {$^\circ $ }$ , and $i=50\hbox {$^\circ $ }$ (red curve), or $i=30\hbox {$^\circ $ }$ (green curve). The starting epoch of the three integrations is $2~453~000.00~{\rm JD}$ . The blue horizontal lines represents the precision of the current set of data, obtained with the CORALIE instrument. The vertical blue line marks the beginning of 2009.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig16} \end{figure}$	Figure 16: Libration and secular frequencies. For each inclination i, we pick a stable orbit with a low $\chi ^2$ and integrate it over 1 million years. The fundamental frequencies are then determined through frequency analysis. The solid curves represent, for each frequency, a $1/\sin^2{i}$ law.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig17} \end{figure}$	Figure 17: Stable configurations for $i_{\rm b} = 90\hbox {$^\circ $ }$ (i.e. for the minimum mass of planet b). For each value of $i_{\rm c}$ in $(90\hbox {$^\circ $ }, 30\hbox {$^\circ $ }, 11\hbox {$^\circ $ })$ , and $\Delta \Omega$ between $0\hbox {$^\circ $ }$ and $360\hbox {$^\circ $ }$ we compute the diffusion index over grids in the $(a_{\rm c},e_{\rm c})$ plane of initial conditions with step sizes of 0.0025 AU and 0.004, respectively. Each grid is also centered on the minimum $\chi ^2$ computed for each pair $(i_{\rm c}, \Delta \Omega )$ . We plot against $\Delta \Omega$ the square root of the minimum $\chi ^2$ in the top panel, and the proportion of stable orbits inside $\sqrt {\chi _{\rm r}^2}= 1.5$ in the middle panel. The mutual inclination corresponding to each triplet $(i_{\rm b}, i_{\rm c}, \Delta \Omega )$ is plotted in the bottom panel.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{figs/13635_fig18}\vspace*{-0.25mm} \end{figure}$	Figure 18: Stable configurations for $i_{\rm b} = 30\hbox {$^\circ $ }$ (i.e. for approximately twice the minimum mass: $m_{\rm b} \simeq 33.5~M_{{\rm jup}}$ ). See Fig. 17 caption for more details.
Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] Anderson, D. R., Hellier, C., Gillon, M., et al. 2010, ApJ, 709, 159 [NASA ADS] [CrossRef] [Google Scholar]

[2] Beaugé, C., Ferraz-Mello, S., & Michtchenko, T. A. 2003, ApJ, 593, 1124 [Google Scholar]

[3] Beaugé, C., Giuppone, C. A., Ferraz-Mello, S., & Michtchenko, T. A. 2008, MNRAS, 385, 2151 [Google Scholar]

[4] Correia, A. C. M., Udry, S., Mayor, M., et al. 2005, A&A, 440, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[5] Correia, A. C. M., Udry, S., Mayor, M., et al. 2009, A&A, 496, 521 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[6] Correia, A. C. M., Couetdic, J., Laskar, J., et al. 2010, A&A, 511, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[7] Fienga, A., Manche, H., Laskar, J., & Gastineau, M. 2008, A&A, 477, 315 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[8] Gayon, J., & Bois, E. 2008, A&A, 482, 665 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[9] Gozdziewski, K., & Maciejewski, A. J. 2001, ApJ, 563, L81 [NASA ADS] [CrossRef] [Google Scholar]

[10] Gozdziewski, K., Konacki, M., & Maciejewski, A. J. 2006, ApJ, 645, 688 [NASA ADS] [CrossRef] [Google Scholar]

[11] Henon, M. 1974, Celest. Mech., 10, 375 [NASA ADS] [CrossRef] [Google Scholar]

[12] Henon, M. 1997, Generating Families in the Restricted Three-Body Problem, ed. M. Henon [Google Scholar]

[13] Ji, J., Li, G., & Liu, L. 2002, ApJ, 572, 1041 [NASA ADS] [CrossRef] [Google Scholar]

[14] Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]

[15] Laskar, J. 1993, Phys. D, 67, 257 [Google Scholar]

[16] Laskar, J. 1999, in NATO ASI Hamiltonian Systems with Three or more Degrees of Freedom, ed. C. Simo (Kluwer), 134 [Google Scholar]

[17] Laskar, J. 2005, in Hamiltonian systems and Fourier analysis, ed. D. Benest, C. Froeschle, & E. Lega (Cambridge Scientific Publishers), 99 [Google Scholar]

[18] Laskar, J., & Robutel, P. 2001, Celest. Mech. Dynam. Astron., 80, 39 [NASA ADS] [CrossRef] [Google Scholar]

[19] Laughlin, G., & Chambers, J. E. 2001, ApJ, 551, L109 [NASA ADS] [CrossRef] [Google Scholar]

[20] Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 [NASA ADS] [CrossRef] [Google Scholar]

[21] Lee, M. H., Butler, R. P., Fischer, D. A., Marcy, G. W., & Vogt, S. S. 2006, ApJ, 641, 1178 [NASA ADS] [CrossRef] [Google Scholar]

[22] Lestrade, J.-F., & Bretagnon, P. 1982, A&A, 105, 42 [NASA ADS] [Google Scholar]

[23] Libert, A.-S., & Henrard, J. 2007, A&A, 461, 759 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[24] Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN, The art of scientific computing, 2nd edn. (Cambridge: University Press) [Google Scholar]

[25] Smith, A. W., & Lissauer, J. J. 2009, Icarus, 201, 381 [NASA ADS] [CrossRef] [Google Scholar]

[26] Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]

[27] Udry, S., Mayor, M., Naef, D., et al. 2002, A&A, 390, 267 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

Dynamical stability analysis of the HD 202206 system and constraints to the planetary orbits

1 Introduction

2 New orbital solution for the HD 202206 system

3 Numerical set-up

3.1 Conventions

3.2 Fitting procedure

3.3 Numerical integrations

3.4 Stability threshold

4 Review of the coplanar edge-on case

4.1 Global dynamics

4.2 Stable fit

4.3 Resonant and secular dynamics

4.4 Test particles

4.5 Finding the center of libration inside a resonance through frequency analysis

4.5.1 Center of libration

4.5.2 Quasi-periodic decomposition

4.5.3 Application to HD 202206

5 Coplanar orbits

5.1 Stability and low orbits

5.2 Resonant and secular behavior

6 Mutually inclined orbits

6.1 kb = 1

6.2 kb = 2

6.3 Conclusions

7 Discussion and conclusion

References

Footnotes

All Tables

All Figures

5.1 Stability and low $\chi ^2$ orbits

6.1 k_b = 1

6.2 k_b = 2