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/00046361/200913635  
Published online  07 September 2010 
Dynamical stability analysis of the HD 202206 system and constraints to the planetary orbits^{}
J. Couetdic^{1}  J. Laskar^{1}  A. C. M. Correia^{1,2}  M. Mayor^{3}  S. Udry^{3}
1  Astronomie et Systèmes Dynamiques, IMCCECNRS
UMR 8028, Observatoire de Paris, UPMC, 77 Avenue DenfertRochereau,
75014 Paris, France
2  Departamento de Física, Universidade de Aveiro, Campus de Santiago, 3810193 Aveiro, Portugal
3  Observatoire de Genève, 51 ch. des Maillettes, 1290 Sauverny, Switzerland
Received 10 November 2009 / Accepted 24 April 2010
Abstract
Context. Longterm, precise Doppler measurements with the
CORALIE spectrograph have revealed the presence of two massive
companions to the solartype star HD 202206. Although the
threebody 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 OC
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 longterm stability
in a small (in terms of OC) 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
stable orbits) are limited with respect to inclinations to the line of sight between
and
.
This limits the masses of both companions to roughly twice the minimum:
and
.
Noncoplanar configurations are possible for a wide range of mutual inclinations from
to
,
although
configurations
seem to be favored. We also confirm the 5:1 mean motion resonance to be
most likely. In the coplanar edgeon case, we provide a very good
stable solution in the resonance, whose
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 planetsearch program in the southern hemisphere has found two companions around the HD 202206 star. The first one is a very massive body with minimum mass (Udry et al. 2002), while the second companion is a 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 lowmass 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 threebody 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 minimum value of the best fit, they concluded that the system should be locked in this 5:1 resonance. Later on, Gozdziewski and coworkers 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 edgeon 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 tradeoff 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 edgeon 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 radialvelocity measurements^{}. Using the iterative LevenbergMarquardt method (Press et al. 1992), we fit the observational data using a 3body Newtonian model, assuming coplanar 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 (Sousa et al. 2008). This fit yields two planets with an adjustment of and , slightly above the photon noise of the instrument, which is around . 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 , and the other at P=1397 day, e=0.104, and a minimum mass of (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 and .
Figure 1: CORALIE radial velocities for HD 202206 superimposed on a 3body Newtonian orbital solution (Table 1). 

Open with DEXTER 
We also fitted the data with a 3body 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 and rms = . 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 setup
3.1 Conventions
In this paper, the subscripts
and
respectively refer to the bodies with the shortest and longest orbital
periods (inner and outer). The initial conditions illconstrained by
the radial velocity data are
(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
matters. In particular, the mutual inclination depends on this quantity, through
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 edgeon coplanar configuration is and . The observer is limited to the velocity projected on the k axis. 

Open with DEXTER 
Thus,
can always be set to
in the initial conditions, which leads to
,
and only three parameters are left free
.
They are connected to the three interesting unknowns of the system, namely the mutual inclination I (Eq. (1)), and the two planetary masses ,
as
where is the star mass, the amplitude of the radial velocity variations, the orbital period, 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 , we can control the mutual inclination with 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 ) 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 to a good approximation. We can define two factors and . With and the minimum masses obtained for the edgeon coplanar case (see Sect. 4 and Table 1), we can write
For a given factor k, two values of the inclination are possible: x and (where ). For instance, and give k=2.
Additionally, for a given pair , the accessible mutual inclinations I are limited. Since the inclinations and are angles between 0 and (excluded), prograde coplanar configurations are only possible for and (Eq. (1)). Similarly, retrograde coplanar configurations are only possible for and . This means that in both cases. More generally, , and
For any value of I, two values of are possible: x and x, since appears through its cosine in Eq. (1). The extrema are obtained for (maximum), and (minimum). Finally we notice that one can restrict one of the two inclinations to the line of sight to , which will be the case for .
3.2 Fitting procedure
The influence of , , and 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 timescale 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 , , and are very poorly constrained by the radial velocity data, we cannot fit the data with a model that includes them. Instead, for any chosen (, , ) set, we compute a new best fit using a LevenbergMarquard minimization (Press et al. 1992) and a threebody model, but with eleven free parameters: the center of mass velocity , and for each planet, the semiamplitude of radial velocity K, the period P, eccentricity e, mean anomaly M, and periastron , all given at the initial epoch. Throughout the paper the initial conditions are given at the same initial epoch .
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: and ;
 and 10 fitted: , , , , and .
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)
where , n is the mean motion, and . These equations are averaged to obtain the following firstorder secular perturbations
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 ) of the mean motions n_{1}, obtained over two consecutive time intervals of length T_{1} = T/2. The stability index 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
.
To that end, we used a second stability index D_{2}. Using the same numerical integration, we computed two new determinations of the mean motion n_{2} and
over two consecutive time intervals of length
T_{2} = T_{1}/k, where k
> 1. In the case of quasiperiodic motion, the diffusion should
be close to zero but it is limited by the precision in determining the
frequencies. Since D_{1} is computed over longer time intervals, the frequencies are better determined, and thus, D_{1} should be, on average, smaller than D_{2}.
On the contrary, for chaotic trajectories, the diffusion will
increase on average for longer time intervals. We can then determine an
approximated value
for which
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 for a particular diffusion grid, we look at the distribution of D_{1} < D_{2} trajectories as a function of D_{1}. We actually work with smoothed values and of D_{1} and D_{2} 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 data in 0.5 wide bins, and compute the percentage of 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.
Figure 3: Distribution of D_{1} < D_{2} trajectories from the top panel of Fig. 6. Each integrated trajectory is binned with respect to its diffusion index , after the diffusion grid has been smoothed. For each bin, we computed the proportion of trajectories with a decreasing diffusion index over time: D_{1} < D_{2}. The histogram shows the results for 0.5 wide bins. As expected from Eq. (7), orbits with a high diffusion index D_{1}, and D_{1} < D_{2} are nearly non existent ( ), while we observe the opposite situation for lowdiffusion index orbits ( ). 

Open with DEXTER 
We choose to define as the D_{1} value for which of the trajectories exhibit D_{1} < D_{2}. Graphically, it is the abscissa for which the curve in Fig. 3 crosses y = 0.99. In this example we get .
Figure 4: Percentage of orbits wrongly flagged asstable (false stable) or unstable (false unstable). Once is chosen using a given percentage threshold (see Fig. 3), we have a criterion for stable ( ) and unstable orbits ( ). We compared the results with a diffusion grid computed over a longer time interval (in this case 2 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 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 , the number of false unstables is nearly zero. This is easy to understand, since a higher percentage threshold implies a lower value for , which in turn leads to flagging very few actually unstable orbits, while we might miss several stable ones, and viceversa. To estimate the number of false stables and unstables, we recomputed the diffusion grid on a longer integration time, 2 40 000 years, and used this grid as a reference (Fig. 4).
4 Review of the coplanar edgeon case
4.1 Global dynamics
Following Correia et al. (2005), we study in more detail the dynamics in the neighborhood of the 3body fit obtained in the case of coplanar orbits with (that is, the system seen edgeon). 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 edgeon, that is, with both inclination at and . We let , , and vary. We always keep constant, as it is much better constrained by the radial velocity data. This implies that when we change , the mean anomaly varies accordingly. In the particular case where , this means that the initial mean longitude is kept constant. For each initial conditions, we compute the diffusion index and the square root of the reduced .
Figure 5: Global view of the dynamics of HD 202206 for variations of the semimajor axis and periastrum of the outer planet ( bottom panel) or semimajor axis and eccentricity ( top panel). The step sizes for , , are , , 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 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 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 region: the dark area around AU and . 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 and of initial conditions. The step sizes for , , are , , and 0.004 respectively. The other parameters were kept constant and taken from the fit S1 (Table 1). The level curves give the value computed for each set of initial conditions. The color scale gives the diffusion index . 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 level curves, at the coordinates marked by a cross, inside a highdiffusion (green) area. Several lowdiffusion (blue and dark blue) zones exist for which the orbits are stabilized either by mean motion resonances or by locking of around . Orbits stabilized by the corotation of the apsidal lines are the blue to black zones around (bottom panel). The width (in the direction) increases with the semimajor axis , from 90 degrees at 2 AU, to nearly 360 degrees at 4 AU, since a wider libration of around 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 and ;
 the 5:1 MMR at 2.5 AU with a notable stable island around where the best fit is located;
 the 1/6 MMR at 2.8 AU with a stable island around or high eccentricity.
4.2 Stable fit
We now take a closer look at the 5:1 mean motion resonance island around AU and , 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 , for , and 0.002 for .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 . 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 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: in the structure on the rim, and in the center.
Figure 6: Global view of the dynamics of HD 202206 for variations of the semimajor axis and periastrum ( bottom), and semimajor axis and eccentricity ( top) of the outer planet. The step sizes are respectively 0.0025 AU, , and 0.002 for the eccentricity. The color scale is the stability index (), and the level curves give the 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 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 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 not significantly worse than the best fit. For instance, the orbital solution S2 given in Table 2 is stable and has a of 1.4136 (marked by a filled white circle). The orbital elements are the same as S1, except for which was adjusted from 2.4832 AU to 2.49 AU.
Table 2: Stable orbital parameters S2 for the HD 202206 system for and .
4.3 Resonant and secular dynamics
The orbital solution S2 was integrated over . It remained stable and displays a regular behavior during the whole time. Using frequency analysis on an integration over , we determined its fundamental frequencies (Table 3). Following our notation, and are the mean motions. The secular frequencies are denoted g_{1} and g_{2}. Finally is the frequency associated with the resonance's critical angle . The fundamental secular frequencies g_{1} and g_{2}, related to the periastron of the inner and outer planets correspond to the periods 10^{3} yr and 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: . As a consequence, one of them is superfluous. A new fundamental frequency associated to the resonance 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
(8) 
The variations of versus time are plotted in Fig. 7 and exhibit a libration around . We nonetheless observe that this libration is the modulation of several different terms with similar amplitudes of approximately , but on different time scales. This leads to a libration with an amplitude that can be higher than .
Figure 7: Time variation of the resonant argument for the orbital solution S2 (green line). is in libration around , with a modulation of two terms of period years, and 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
more accurately, we searched for a quasiperiodic decomposition of .
We started with a frequency decomposition using frequency analysis
(Laskar 2005) as
We then decomposed each frequency on the four fundamental frequencies (, g_{1}, g_{2}, )
where , and are integers. Each term of the decomposition follows a D'Alembertlike relationship expressed by
which can be used to simplify the expression of . Indeed, by rearranging the righthand side of Eq. (10), we can write
and using Eq. (11), we get
Table 4: Quasiperiodic decomposition of the resonant angle 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.
Figure 8: Mean motion diffusion of test particles. integration S2 solution with massless particles over 16 000 yr. We computed two determination n and of each particules mean motion over two consecutive time intervals of 8000 yr. The color scale codes the stability index . Initial conditions for the massless particles are . 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 semimajor 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 longterm oscillations with frequency g_{1}  g_{2}. The period corresponding to g_{1}  g_{2}, associated to the angle , is years. It is worth mentioning that the angle 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 shortterm oscillations of frequency , corresponding to a period years. More precisely, the libration of is made of three different kinds of contributions:
 secular terms: , hence of the form ;
 resonant terms: and , of the form ;
 and shortperiod terms: .
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 of each particle mean motion over two consecutive time intervals of 8000 years. We then obtain a stability index 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 semimajor 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 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 semimajor 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 Neptunesized 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 and 3 , depending on the phase. We conclude that a yet undetected planet smaller than half a Jupiter mass could exist at semimajor 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 threebody problem and, more particularly, the planetary problem with a p:p+1 meanmotion resonance. The problem is to find the orbit center of libration starting from a quasiperiodic orbit in the resonance.
In the restricted case, this center of libration is a welldefined 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 quasiperiodic orbits live on 4torus in the phase space, although it can be restricted to 3 degrees of freedom using the angular momentum reduction. Each dimension of the 4torus is associated to one of the four fundamental frequencies (f_{0},f_{1},f_{2},f_{3}) of the orbit. Supposing that f_{0} is the frequency of the resonant mode, the orbit equivalent to the center of libration is living on a 3torus, depending only on (f_{1},f_{2},f_{3}). 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 4torus 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 Quasiperiodic 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 threebody problem (see for instance Henon 1997,1974).
We present here a new numerical method of finding the quasiperiodic
center of libration, using the fact that we can get an accurate
quasiperiodic decomposition of a numerically integrated quasiperiodic
orbit with frequency map analysis (Laskar 2005). Let
be a state vector of such an orbit. Using frequency analysis, we can
obtain a quasiperiodic representation for any component x of the vector
for each coordinate of the vector , with , f = (f_{0}, f_{1}, f_{2}, f_{3}), and . We can separate these sums in two parts: x(t) = u(t) + v(t) where u does not depend on the resonant frequency f_{0}, and v has all the terms depending on f_{0}.
The quasiperiodic orbit described by is precisely lying on a 3torus 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 as a new initial condition, and obtain a new resonant quasiperiodic orbit with
The amplitude of the resonant terms in this new orbit will be smaller. In other words, it lives on 4torus closer to the 3torus of the center of libration. We can then iterate this procedure to suppress all the terms with f_{0}.
Figure 9: Schematic representation of the 3torus of the center of libration. We represent the 4torus of a given resonant orbit by a 2torus (a doughnut) as it is not possible to represent it otherwise. The center of libration 3torus 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 state vector, where and . Each step p of the method is decomposed as follows:
 numerical integration of ;
 determination of (f_{0},f_{1},f_{2},f_{3}) using frequency analysis;
 quasiperiodic decomposition: ;
 new initial conditions: x^{(p+1)}(0) = x^{(p)}(0)  v^{(p)}(0).
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 plane. No resonant terms are found in the first 100 terms of the quasiperiodic decomposition of each variable at the last step (with the exception of ). In fact, in half of the variables, there are no resonant terms left in the first 300 terms (Fig. 11).
Figure 10: Projection of a section of the 4torus of each step's trajectory in the plane. Each torus is approximated using a (truncated) quasiperiodic decomposition of the trajectory. 

Open with DEXTER 
Figure 11: Evolution of the amplitude of the libration mode in the quasiperiodic decomposition of each variable at each step. In the top panel we plot the relative amplitude of the first term depending on f_{0} compared to the first nonconstant 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 3torus it is living on since its
quasiperiodic decomposition gives us a parametrization of the torus.
The three angular variables of the torus appears in the decomposition
(17) 
where , , and are initial phases. With , we can rewrite Eq. (14) to reveal the underlying torus
(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 in the three initial phases , , and .
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 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 ( 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 orbits
In this section we investigate the system behavior when both planets remain in the same orbital plane ( ) and are prograde, but the inclination to the line of sight is lower than :
 ;
 .
In this configuration, the masses of the two companions grow approximately in proportion with when the inclination i diminishes. There are thus only small changes between and , 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 orbits, the ones in the 5:1 mean motion resonance remain stable for inclinations up to . For , 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: and .
Figure 12: Dynamics of a coplanar HD 202206 system for different values of the inclination i. Each panel is a diffusion map in the plane of initial conditions constructed the same way as the top panel of Fig. 5 ( the topleft 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 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 , and for lower inclinations, the low region is slightly shifted towards values that are lower than the resonance island location.
To have a more precise picture, we spanned i values from to with a step size of . For each inclination value, we started by computing a new best fit with the LevenbergMarquard algorithm. The values and stability indexes D_{1} and D_{2} are computed in the plane of initial conditions, around the best fit. The size of each grid is 101 161 dots, with step sizes of 0.0025 AU and 0.002. For each of these maps, we computed (see Sect. 3.4), and detected stable regions and their relative positions to the observations in terms of .
To get a synthetic vision of all this data, we can get an estimation of the lowest stable for each inclination (Fig. 13). We plotthe lowest of 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 orbits is really small from to . Is is actually slightly better between and .
To have a better idea of the correlation between low orbits and stable orbits, we looked at the percentage of stable orbits inside a given level curves (Fig. 14). It gives a synthetic representation of the overlap between low regions and stable regions. It is very clear that inclinations lower than are very unlikely, although possible. It also appears that inclinations between and are the most probable.
If we limit the acceptable inclinations to , we can derive limits for the masses of the inner and outer planets (assuming a coplanar prograde configuration):
 ;
 .
Figure 13: Evolution of the best fit (solid curve) and an estimation of the best stable orbits as a function of the inclination i. For each value of i between and , a new fit is computed, assuming a coplanar configuration. We then look for stable orbits around each of these orbital solutions in a plane of initial conditions. 

Open with DEXTER 
Figure 14: Percentage of stable orbits inside a given level curve against i (see Fig. 13). 

Open with DEXTER 
Figure 15: Differences between the radial velocity of a stable solution with , and (red curve), or (green curve). The starting epoch of the three integrations is . 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
Figure 16: Libration and secular frequencies. For each inclination i, we pick a stable orbit with a low and integrate it over 1 million years. The fundamental frequencies are then determined through frequency analysis. The solid curves represent, for each frequency, a 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 value, and plotted its fundamental frequencies , g_{1}, and g_{2} 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 . 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.
Figure 17: Stable configurations for (i.e. for the minimum mass of planet b). For each value of in , and between and we compute the diffusion index over grids in the plane of initial conditions with step sizes of 0.0025 AU and 0.004, respectively. Each grid is also centered on the minimum computed for each pair . We plot against the square root of the minimum in the top panel, and the proportion of stable orbits inside in the middle panel. The mutual inclination corresponding to each triplet 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 and to vary independently and allow variations in , the longitude of the ascending node of c. To span the possible values for , , and in an efficient manner, we restrict ourselves to two mass ratios (Eq. (3)) for planet b
 ( );
 ( );
 ( );
 ( or );
 ( or ).
6.1 k_{b} = 1
For ( , Fig. 17) we find that significant stable zones inside exist mostly for aligned and antialigned ascending nodes (i.e. and in Fig. 17b). For (red curves) the aligned configurations ( ) are coplanar and prograde ( ), and the antialigned configurations ( ) are coplanar and retrograde. Outside those two particular cases, we find no significant stable zones for . Indeed, while the resonant island is roughly centered on the lowest level curve, it is not stable outside the coplanar configurations. We find that an extended zone of stability exists outside the resonance for , but it lies just outside the 1.5 level curve. Due to symmetries, the situation for mirrors that of .
For (green curves) we again find stable zones for , but not around . However stable regions with potential solutions exist for values up to , corresponding to mutual inclinations between and . This is mostly due to the minimum getting smaller up to (see the green curve in Fig. 17a). While the stable regions, both in the resonance and outside, shrink when augments, the level curves encompass a larger area. Finally, for (blue curves), no significant stable zones are found at low values.
Figure 18: Stable configurations for (i.e. for approximately twice the minimum mass: ). See Fig. 17 caption for more details. 

Open with DEXTER 
To summarize, potential solutions (stable orbit with a low ) for mainly exists for coplanar configurations, inside the 5:1 mean motion resonance. If planet c's mass is kept low ( ), stable regions do exist for non coplanar configurations, but they are located outside the lowest region. A noteworthy exception is the retrograde configuration, where stable orbits are only found close to the coplanar case, which only happens for and .
6.2 k_{b} = 2
When we double the mass of planet b ( , Fig. 18), once again retrograde potential solutions can only occur for coplanar orbits: and (green dotted curve). Concerning prograde orbits, potential stable solutions exist for a mutual inclination of :
 and (red curve);
 and (green solid curve).
6.3 Conclusions
To sum up, we find that the configurations that have a significant stable zone at low values are mostly found when the two lines of nodes are aligned. That is for close to or . In addition, stable orbits with the lowest 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 higher than 1.45. Retrograde configurations seem to be limited to nearly coplanar orbits with antialigned ascending nodes. Other than that, we could not find any clear correlation with mutual inclination.
That we find nonresonant 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 longterm stability point of view. This hypothesis has been reinforced by the recent detection of WASP17b, an ultralow 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: . This means that the companions' masses are most likely not greater than twice their minimum values:
 ;
 .
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 occurs for particular initial values of this angle, providing a stabilizing mechanism outside the meanmotion resonance, not far from the best fits. In most cases other than retrograde coplanar configurations, those orbits in nearcommensurability 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 OC are in the 5:1 mean motion resonance. In fact, the minimum is almost always in the resonance or very close to it, and stable orbits in the resonance can be found with not significantly higher than the best fit. In addition the OC 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 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 edgeon 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 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 semimajor axis lower than 0.12 AU, and that 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.
We acknowledge support from the Swiss National Research Found (FNRS), French PNPCNRS, Fundação para a Ciência e a Tecnologia, Portugal (PTDC/CTEAST/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., FerrazMello, S., & Michtchenko, T. A. 2003, ApJ, 593, 1124 [Google Scholar]
 Beaugé, C., Giuppone, C. A., FerrazMello, 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 ThreeBody 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 [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.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/519/A10  ... measurements^{}
 CORALIE data were not taken during 20072008. 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 corresponds to on the variation in the parameters of the bestfit solution S1.
All Tables
Table 1: Best Newtonian fit S1 for the HD 202206 system assuming and .
Table 2: Stable orbital parameters S2 for the HD 202206 system for and .
Table 3: Fundamental frequencies for S2.
Table 4: Quasiperiodic decomposition of the resonant angle 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
Figure 1: CORALIE radial velocities for HD 202206 superimposed on a 3body Newtonian orbital solution (Table 1). 

Open with DEXTER  
In the text 
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 edgeon coplanar configuration is and . The observer is limited to the velocity projected on the k axis. 

Open with DEXTER  
In the text 
Figure 3: Distribution of D_{1} < D_{2} trajectories from the top panel of Fig. 6. Each integrated trajectory is binned with respect to its diffusion index , after the diffusion grid has been smoothed. For each bin, we computed the proportion of trajectories with a decreasing diffusion index over time: D_{1} < D_{2}. The histogram shows the results for 0.5 wide bins. As expected from Eq. (7), orbits with a high diffusion index D_{1}, and D_{1} < D_{2} are nearly non existent ( ), while we observe the opposite situation for lowdiffusion index orbits ( ). 

Open with DEXTER  
In the text 
Figure 4: Percentage of orbits wrongly flagged asstable (false stable) or unstable (false unstable). Once is chosen using a given percentage threshold (see Fig. 3), we have a criterion for stable ( ) and unstable orbits ( ). We compared the results with a diffusion grid computed over a longer time interval (in this case 2 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 
Figure 5: Global view of the dynamics of HD 202206 for variations of the semimajor axis and periastrum of the outer planet ( bottom panel) or semimajor axis and eccentricity ( top panel). The step sizes for , , are , , 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 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 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 region: the dark area around AU and . 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 
Figure 6: Global view of the dynamics of HD 202206 for variations of the semimajor axis and periastrum ( bottom), and semimajor axis and eccentricity ( top) of the outer planet. The step sizes are respectively 0.0025 AU, , and 0.002 for the eccentricity. The color scale is the stability index (), and the level curves give the 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 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 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 
Figure 7: Time variation of the resonant argument for the orbital solution S2 (green line). is in libration around , with a modulation of two terms of period years, and years, amplitudes of about 35 and 50 degrees, respectively. The black line shows the first secular term contribution. 

Open with DEXTER  
In the text 
Figure 8: Mean motion diffusion of test particles. integration S2 solution with massless particles over 16 000 yr. We computed two determination n and of each particules mean motion over two consecutive time intervals of 8000 yr. The color scale codes the stability index . Initial conditions for the massless particles are . 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 semimajor 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 
Figure 9: Schematic representation of the 3torus of the center of libration. We represent the 4torus of a given resonant orbit by a 2torus (a doughnut) as it is not possible to represent it otherwise. The center of libration 3torus is then represented by the circle in the center of the interior of the doughnut. 

Open with DEXTER  
In the text 
Figure 10: Projection of a section of the 4torus of each step's trajectory in the plane. Each torus is approximated using a (truncated) quasiperiodic decomposition of the trajectory. 

Open with DEXTER  
In the text 
Figure 11: Evolution of the amplitude of the libration mode in the quasiperiodic decomposition of each variable at each step. In the top panel we plot the relative amplitude of the first term depending on f_{0} compared to the first nonconstant 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 
Figure 12: Dynamics of a coplanar HD 202206 system for different values of the inclination i. Each panel is a diffusion map in the plane of initial conditions constructed the same way as the top panel of Fig. 5 ( the topleft 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 
Figure 13: Evolution of the best fit (solid curve) and an estimation of the best stable orbits as a function of the inclination i. For each value of i between and , a new fit is computed, assuming a coplanar configuration. We then look for stable orbits around each of these orbital solutions in a plane of initial conditions. 

Open with DEXTER  
In the text 
Figure 14: Percentage of stable orbits inside a given level curve against i (see Fig. 13). 

Open with DEXTER  
In the text 
Figure 15: Differences between the radial velocity of a stable solution with , and (red curve), or (green curve). The starting epoch of the three integrations is . 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 
Figure 16: Libration and secular frequencies. For each inclination i, we pick a stable orbit with a low and integrate it over 1 million years. The fundamental frequencies are then determined through frequency analysis. The solid curves represent, for each frequency, a law. 

Open with DEXTER  
In the text 
Figure 17: Stable configurations for (i.e. for the minimum mass of planet b). For each value of in , and between and we compute the diffusion index over grids in the plane of initial conditions with step sizes of 0.0025 AU and 0.004, respectively. Each grid is also centered on the minimum computed for each pair . We plot against the square root of the minimum in the top panel, and the proportion of stable orbits inside in the middle panel. The mutual inclination corresponding to each triplet is plotted in the bottom panel. 

Open with DEXTER  
In the text 
Figure 18: Stable configurations for (i.e. for approximately twice the minimum mass: ). See Fig. 17 caption for more details. 

Open with DEXTER  
In the text 
Copyright ESO 2010
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.