Searching for stable orbits in the HD 10180 planetary system

Abstract. A planetary system with at least seven planets has been found around the star HD 10180. However, the traditional Keplerian and n-body fits to the data provide an orbital solution that becomes unstable very quickly, which may quest the reliability of the observations. Here we show that stable orbital configurations can be obtained if general relativity and long-term dissipation raised by tides on the innermost planet are taken into account.


Introduction
A planetary system comprising at least five Neptune-like planets with minimum masses ranging from 12 to 25 M ⊕ has been found orbiting the solar-type star HD 10180 at separations between 0.06 and 1.4 AU (Lovis et al. 2010).The trace of a sixth planet is present at a longer period, probably due to a 65 M ⊕ object.Moreover, another body with a minimum mass as low as 1.4 M ⊕ may be present at 0.02 AU.
The stability of such system is not straightforward, in particular taking into account that the minimum masses of the planets are of the same order as Neptune's mass and the eccentricities are relatively high when compared with the eccentricities of the planets in the Solar System.As a consequence, mutual gravitational interactions between planets in the HD 10180 system cannot be neglected and may give rise to instability.Indeed, if we perform a quick simulation of the evolution of the system directly obtained with a Keplerian or n-body adjustment to the data, we verify that the inner planet's eccentricity grows very fast leading to eccentricities close to 1 in less than 1 Myr (Fig. 1).

The secular planetary equations
Over long times, and in absence of strong mean motion resonances, the variations of the planetary elliptical elements are well described by the secular equations, that is, the equations obtained after averaging over the longitudinal motion of the planets (Laskar 1990).The secular system can even be limited to its first order and linear part, which is usually called the Laplace-Lagrange system (Laskar 1990) which can be written using the classical complex notation where with k = b, c, ..., h.This linear equation is classically solved by diagonalizing the real matrix (L) through the linear transformation to the proper modes Using the initial conditions from a Keplerian fit to the data (Lovis et al. 2010), we have computed analytically the (real) Laplace-Lagrange matrix (L), and derived the (real) matrix (S) of its eigenvectors which gives the relation from the original eccentricity variables z k to the proper modes u k .After the transformation (Eq.3), the linear system (Eq. 1) becomes diagonal where g 1 , . . ., g 7 are the eigenvalues of the Laplace-Lagrange matrix (L).The solution is then trivially given for all k = 1, . . ., 7 by while the secular solution is obtained through (Eq.3).It can be noted that as the matrix (L), and thus (S), only depend on the masses and semi-major axis of the planets, they do not change much for different fits to the data because the periods and masses are well constrained (for a given inclination of the system).

Stability of the short-period planet b
Despite the proximity of planet b to the star, it undergoes strong gravitational perturbations from the remaining bodies, due to secular interactions.Even for a model where the initial eccentricity of the planet b is initially set at zero, its eccentricity shows a rapid increase, which can reach up to 0.9 in only 1 Myr (Fig. 1).The inclusion of general relativity in the model calms down the eccentricity evolution, but still did not prevent the eccentricity of planet b to attain values above 0.4, and planet c to reach values around 0.3, which can then destabilizes the whole system.Most of these variations are well described by the secular system (Eqs.1, 3 and 4).In particular, the effect of general relativity will be to largely increase the first diagonal terms in the Laplace-Lagrange matrix (L), which will increase in z 1 the contribution of the proper mode u 1 and thus decrease the long time oscillations due to the contribution of the modes.
At this stage, one may question the existence of the innermost short-period planet.However, since the mass of this planet should be in the Earth-mass regime, we may assume that it is mainly a rocky planet.As a consequence, due to the proximity of the star, this planet should undergo intense tidal dissipation, that continuously damp its orbital eccentricity.

05001-p.2
Figure 1.: Evolution of the eccentricity of planet HD10180b during 250 kyr for two different models.In red we used a n-body model and the initial eccentricity of planet HD10180b is set at zero, but mutual gravitational perturbations increase its value to nearly 0.9 in 250 kyr.The green curve is the eccentricity of the same planet HD10180b for a model where general relativity is added and the fit is made using the additional constraint (Eq.14) that takes into account tidal dissipation.This last solution is stable for at least 150 Myr.

Tidal contributions
Using a simplified model, the tidal variation of the eccentricity can be written as (Correia 2009): where f 6 (e) = 1 + 45e 2 /14 + 8e 4 + 685e 6 /224 + 255e 8 /448 + 25e 10 /1792 and M ⋆ is the mass of the star, m and R the mass and the radius of the planet, k 2 the second Love number, and Q the dissipation factor.As for the Laplace-Lagrange linear system, we can linearize the tidal contribution from expression (6), and we obtain for each planet k a contribution that is, an additional real contribution on each diagonal term of the Laplace-Lagrange which thus adds an imaginary part iγ k to the diagonal terms of (L).In fact, as apart from the planet b, all planetary masses are relatively large, the dissipation in these planets may be small, and we will uniquely consider the tidal dissipation in the innermost planet b.It should nevertheless be stressed that all the eigenvalues of the matrix (L) will be modified, and will present an imaginary part that will induce an exponential term in the amplitude of the proper modes: Adopting a value similar to Mars k 2 /Q = 0.0015 (which is a minimal estimation for the dissipation), we can see in Figure 2 that the amplitude of the proper mode u 1 will be rapidly damped in a few tens of Myrs, whatever is its original value.Over more than 1 Gyr, the amplitude of the second proper mode u 2 is also most probably damped to a small value.
In order to obtain a realistic solution that is the result of the tidal evolution of the system, it is thus not sufficient to impose that the innermost planets have small eccentricity, as this may only be realized at the origin of time (Fig. 1).It is also necessary that the amplitude of the first two proper modes, and particularly u 1 are set to small values.This will be the way to obtain a solution with moderate variations of the eccentricities on the innermost planets, which will then increase its stability.

Orbital solution with dissipative constraint
As the proper modes of the two innermost planets, u 1 and u 2 are damped after about 1 Gyr, we expect them to be small for the present observations (the age of the star is estimated to be ∼ 4 Gyr).The initial conditions for the HD 10180 planetary system should take into account the tidal damping.By inverting the system (Eq.3)we get which allows us to modify our fitting procedure in order to include a constraint for the tidal damping of the proper modes u 1 and u 2 : For that purpose, we added to the χ 2 minimization, an additional term, corresponding to these proper modes: where R is a positive constant, that is chosen arbitrarily in order to obtain a small value for u 1 and u 2 simultaneously.
Using R = 350 we get u 1 ∼ 0.0017 and u 2 ∼ 0.044 and obtain a final χ 2 = 1.24, which is statistically equivalent to the results obtained without this additional dissipative constraint: R = 0 and χ 2 = 1.22 (Lovis et al. 2010).We believe that this solution is a better and more realistic representation of the true system than the one obtained from the traditional Keplerian or n-body fits.Indeed, with this constraint, the variation of the eccentricity of the two innermost planets are now much smaller (Fig. 1).

Long-term orbital evolution
The estimated age of the HD 10180 system is about 4.3 Gyr (Lovis et al. 2010), indicating that the present planetary system had to survive during such a long time-scale.We tested the stability of the system by performing a numerical integration over 150 Myr of the orbits of the six outermost planets determined with the dissipative constraint from last section, using the symplectic integrator SABAC4 of Laskar & Robutel (2001) with a step size of 10 −3 years, including general relativity.
The orbits remain very stable throughout the simulation.We have also integrated the full 7-planet system over 1 Myr with a step size of 10 −4 years, without any sign of strong instability, although the frequency analysis of the solutions with seven planets shows that these solutions are not as well approximated by quasi-periodic series as the 6-planet solutions.This will have to be analyzed further as like in our solar system, the presence of this innermost planet seems to be critical for long time stability of the system (Laskar & Gastineau 2009).
The fact that we are able to find a stable solution compatible with the observational data can still be considered as a good indicator of the reliability of the determination of the HD 10180 planetary system.As the secular frequencies g k are relatively large, the secular variations of the orbital parameters vary more rapidly than in our Solar System, which should enable them to be detected directly from observations, and hence access the true masses and mutual inclinations of the planets as it was done for the system GJ 876 (Correia et al. 2010).

Conclusions
The HD 10180 star hosts a packed system of at least seven planets which undergo significant mutual gravitational perturbations.The orbital parameters obtained from a Keplerian or n-body fit to the observational data do not provide stable long-term solutions, since the adjustment algorithms tend to overestimate the eccentricities of the planets, which enhances the gravitational interactions.
Since the innermost planet has a small mass and is very close to the host star, it is expected to undergo very intense tidal effects, that will damp its eccentricity close to zero (see Mardling, 2007).However, even if one fixes the eccentricity of the innermost planet at zero during the fitting procedure, we observe that its eccentricity quickly increases to very high and unstable values.It is also useless to compute the direct effect of tides in the n-body model used for adjusting the observational data, since this effect is only visible for long time-scales, and not for the short span of a few years of observations.
In order to obtain a stable orbital solution for the whole system, the effect of tidal friction must indeed be taken into account, but in an indirect way.First, one can solve the secular system and determine which eccentricity modes u k are damped through the age of the star.Then, using the above information, one can use the constraint given by u k ≈ 0 to determine which linear combinations of the eccentricities of all planets should be close to zero.For that purpose we add an additional term to the χ 2 minimization, forcing those modes to be small.
We have applied this methodology to the HD 10180 planetary system with success, since the out-coming orbital solution is very stable over very long time-scales.We conclude that this indirect way of computing the tidal friction is a very powerful way of getting a better determination of the true orbital parameters of the planets.Therefore, we suggest that it should be generally applied to planetary system for which at least one planet undergoes significant tidal interactions with the star.
We stress that even though tidal friction may act only over one planet, the fact that the system is gravitationally coupled will result that all eccentricities are damped.Therefore, the presence of an Earth-mass body at 0.022 AU in the HD 10180 system has important implications for its dynamics and highlights the role of tidal dissipation to guarantee stability.Since this planet is not yet fully confirmed (it has a false alarm probability of 1%), some efforts should still be done on its confirmation.

Figure 2 .:
Figure 2.: Tidal evolution of the amplitude of the proper modes u 1 (red), u 2 (green), u 3 (blue), and u 4 (pink) resulting from the tidal dissipation on planet b with k 2 /Q = 0.0015.