© ESO, 2012

1. Introduction

In several planetary systems, some of the planets are very close to their central star and thus subject to strong tidal interaction. If a planet were alone around its star, a circularization of its orbit would result. However, for multi-planet system, owing to the secular interaction between the planets, the final evolution of the system may be different, and affected by residual eccentricity (Wu & Goldreich 2002; Mardling 2007; Batygin et al. 2009; Mardling 2010; Lovis et al. 2011). This explains why the fitting of a circular orbit to the innermost planets in a system subject to tidal dissipation will not ensure that its eccentricity remains small, because the secular interactions may cause the eccentricities to increase to large values (Lovis et al. 2011). One thus needs to take into account that the observed system is the result of a tidal process (Lovis et al. 2011). Here, we develop in full detail the method that has been used in Lovis et al. (2011) to fit a tidally evolved system. The theory is very general and is compared to the previous results of Mardling (2007, 2010).

2. Model

2.1. Lagrange-Laplace secular equations

In the absence of mean-motion resonances, the long-term evolution of a planetary system is approximated well by the Lagrange-Laplace secular equations, obtained after averaging the motion of the planets over their mean longitudes, and retaining only the terms of first order in eccentricity and inclination. This model is well suited to systems of low eccentricity and inclination such as our Solar System, but even for moderate eccentricity, it usually provides a good approximation to the understanding of the essential features of the long-term evolution of the orbits. A detailed account of the derivation of these equations and their use for planetary systems with two or more planets can be found in Murray & Dermott (1999).

In the secular approximation, the semi-major axes are constant (Laplace 1785). There is thus no exchange of energy among the planetary orbits, and only an exchange of angular momentum that can lead to slow (typical periods of a few thousand years) but significant variations in the eccentricities and inclinations of the planets that are well described by the Lagrange-Laplace first order equations (see Laskar 1990). In this linear approximation, inclinations and eccentricities are decoupled and follow the same kind of evolution. For simplicity, we focus in this short paper only on coplanar systems.

We define n to be the number of planets. Using the classical complex variables, z_k = e_ke^iϖ_k, for k = 1,...,n, where e_k and ϖ_k are, respectively, the eccentricity and the longitude of the periastron of the kth planet (see Fig. 1a), the secular equations, limited to first order in eccentricity read¹(1)and where A is a real matrix whose elements are (Laskar & Robutel 1995) (2)In these expressions, m_k and a_k, are the mass and semi-major axis of the kth planet, while the mean motion n_k is defined by . We assume that the planets are ordered by increasing semi-major axis, while the index 0 represents the star. The functions C₂(α) and C₃(α) are defined by mean of , the usual Laplace coefficients (e.g. Laskar & Robutel 1995), as (3)Using the definition of the Laplace coefficients, one gets (4)As the semi major axes a_i are constant in the secular approximation, the coefficients A_jk of matrix A are constant and the evolution of eccentricities and longitudes of perihelion is provided by the linear equation with constant coefficients (1). If the matrix A were diagonal, Eq. (1) would reduce to dz_k/dt = iA_kkz_k, which has the solution z_k(t) = z_k(0)exp(iA_kkt). The eccentricity would then remain constant with a uniform perihelion precession (Fig. 1b).

Fig. 1 a) The complex variable z = eexp(iϖ) provides the orientation and eccentricity of the orbit in the plane. Alternatively, one can use the rectangular real coordinates k = ecosϖ, h = esinϖ, as in Murray & Dermott (1999). b) If the Lagrange-Laplace matrix A were diagonal, the eccentricity of all planets would be constant with uniform precession. c) As A is not diagonal, one needs to make a linear transformation to transform zk into proper modes uk. The proper modes amplitude are then constant and uk precess uniformly with frequency gk. d) The full solution zj = ejexp(ϖj) are linear combinations of the proper modes uk. In the figure, with two proper modes (i.e. two planets and n = 2). The red curve represents the evolution of the eccentricity and longitude of perihelia of Jupiter over 200 000 years under the perturbation of Saturn.

In general, although the diagonal terms of A are large, the matrix A is not diagonal, and gravitational coupling exists among the planetary orbits. Equation (1) is then classically solved by diagonalizing the matrix A through a linear transformation (5)where S = (S_jk) is an n × n real matrix composed of the eigenvectors of A. With these new variables, the equations of motion become (6)is the diagonal matrix diag(g₁,...,g_n) of the eigenvalues g_k of A. We then have du_k/dt = ig_ku_k, thus (7)As the eigenvalues g_k are real, each proper mode u_k describes a circle in the complex plane at a constant frequency g_k with radius |u_k(0)| (Fig. 1c). The evolution of the planetary eccentricities are then given by (5). Each z_j is thus a linear combination of the proper modes u_k of the form (8)that is on the form (9)where α_jk is real. In Fig. 1d, we represented in red the evolution of z_J = e_Jexp(iϖ_J) for Jupiter in the Solar System. In this case, the solution can essentially be reduced to the first two proper modes (see Laskar 1990), which are represented in the figure. In this case, |α_J2| ≪ |α_J1| and the eccentricity can never reach zero, but for the Earth, in the same decomposition, the leading term amplitude is not larger than the sum of the others, and zero eccentricity can be reached (see Laskar 1990).

2.2. Other effects

The above secular Eqs. (1) describe only the Newtonian interactions between point mass planets. To study exoplanetary systems with short period planets, it is often necessary to add corrections for relativity, oblateness, and tidal friction, at least to the inner planets. The spatial secular equations of motion resulting from all these effects are given in (Lambeck 1980; Eggleton & Kiseleva-Eggleton 2001; Ferraz-Mello et al. 2008) or in Correia et al. (2011) for notations coherent with the present ones. To first order in eccentricity and in the planar case, these corrections modify the diagonal terms of the matrix A in two different ways. There are conservative terms that are purely imaginary, and dissipative ones that are real (and negative). We thus define two new diagonal matrices δA =  ∑ _i = 1,5δA⁽ⁱ⁾ and δB =  ∑ _i = 4,5δB⁽ⁱ⁾ such that the full secular evolution is given by (10)with A_tot = A + δA. The effect of relativity on the kth planet is conservative, and in first order in eccentricity, it leads to (11)The effect of the oblateness of bodies generated by their proper rotation are (12)and (13)for the oblateness of the kth planet and the oblateness of the star, where k_2,k, ω_k, and R_k are, respectively, the second Love number, the proper rotation rate, and the radius of the kth body. Tidal effects have two contributions. With the same approximation, we have (14)where (15)for the tides raised on the kth planet by the star, and (16)where (17)for the tides raised on the star by the kth planet. We consider here the “viscous” approach (Singer 1968; Mignard 1979), where the quality factor of the kth body is Q_k ≡ (n_k(Δt)_k)^-1 and (Δt)_k is a constant time lag.

3. Resolution

3.1. Conservative case

When there is no dissipation (δB = 0) the system (10) is similar to (1) and as above is solved by diagonalizing the matrix A_tot through a linear transformation (18)In the new variables, the equations of motion become (19)is the diagonal matrix diag(g₁,...,g_n) of the eigenvalues g_k of A_tot. We then have²(20)The evolution of the planetary eccentricities are then given by (18). They are linear combinations of the proper modes. The only differences between A and A_tot are their diagonal terms. Those of A_tot are larger or equal to those of A. As a consequence, using A_tot instead of A makes the fundamental frequencies g_k higher and the coupling between the proper modes weaker. The evolution of the eccentricity of each planet is then almost given by one single proper mode, the other modes generate only small oscillations (Fig. 1d).

3.2. General solution

In the full linear secular Eq. (10), the matrix that has to be diagonalized is now iA_tot − δB, where the dissipation part δB comes only from tides. In general, the elements of δB are much smaller than those of the diagonal of A_tot. The correction δB will thus be considered to be a perturbation of the conservative evolution given by A_tot. We define (21)to be the matrix of the linear transformation that diagonalizes the full system. Since δB is a perturbation of A_tot, we hypothesize that the matrix δS₁ is also a perturbation of the matrix S₀. To first order, the inverse of S is (22)and the new diagonal matrix is D = iD₀ − δD₁, with (23)where the bracket is defined by  [δS₁,D₀]  = δS₁D₀ − D₀δS₁. For δD₁ to be actually diagonal, δS₁ is given by (24)As D₀ is diagonal, all terms in the diagonal of  [δS₁,D₀]  vanish. Thus, the diagonal terms of δS₁ do not appear in the computation of δD₁ (23) and they can be set to zero. Let δD₁ = diag(γ₁,...,γ_n). From (23), we then have (25)The coefficients γ_k are real and positive. It turns out that the imaginary part of D is still the one of the conservative case iD₀ (19). The proper frequencies g_k are unaffected by the dissipation δB. However, each proper mode now contains a damping factor γ_k given by (25). The equations of motion in the new variables now read (26)and the solutions are (27)It should be stressed that even if only one planet undergoes tidal dissipation (only (δB)₁₁ is different from 0 for example), because of the linear transformation S₀, all the eigenmodes can be damped (25).

4. Two planet case

In a simpler two planet system where only the first planet undergoes tidal friction, δB = diag(γ,0), the two proper frequencies are given by (28)where T and Δ are the trace and determinant of A_tot. From (25), it can be shown that the two dissipation factors are (29)The sum γ₁ + γ₂ is equal to γ. There is thus always one eigenmode damped on a timescale shorter than 2γ^-1, while the other is damped on a timescale longer than 2γ^-1. In the particular case where A₁₁ = A₂₂, we have γ₁ = γ₂ = γ/2.

Once the first eigenmode is damped, the ratio of the two eccentricities and the difference between the two longitudes of periastron are deduced from (18). We have (30)We note that the matrix δS₁ introduces small corrections to the difference between the longitudes of periastron that are not taken into account in (30).

Fig. 2 Tidal effects on the eccentricity of HD 209458b with one or two companions. a) The companion is a 0.1   MJ planet at 0.4 AU with e2 = 0.4 as in Mardling (2007). The blue curve was obtained without considering the conservative effect of the tides ( in Eq. (14)). The red curve in a), and all the evolutions presented in the subfigures b)–d) take into account this effect. b) The mass of the companion is set to 0.608   MJ to recover the final excentricity of Mardling’s simulation. The red curve is the result of a numerical integration of the full secular equations that are exact in eccentricity. The green one is the analytical solution of the linearized problem. c) Same as b) except for the initial Q-value of HD 209458b, which is set to 15.15 to enable the visualisation of both the damping and the oscillation of the eccentricity. d) Same as c) with an additional 0.1   MJ companion at 1.0 AU with e3 = 0.1.

Fig. 3 Eccentricity e2 of the hypothetical companion of HD 209458b with e1 = 0.01 assuming that the eigenmode with the shortest dissipation timescale is damped (black curves in panel a) and b)). a) The mass of the companion is fixed to m2 = 0.2   MJ. Negative values of e2 correspond to Δϖ = 180deg while positive ones mean Δϖ = 0deg. The dotted line is the eccentricity that the companion would have had 5.5 Gyr ago assuming a dissipation factor computed with (29). b) Same as a) for different masses m2 while the semi-major axis is fixed and set to a2 = 0.4 AU. c) Stellar reflex velocity due to the companion at periastron with the eccentricity of a). d) Idem for the eccentricity of b). In grey regions, the eccentricity of the companion should have been larger than 1 in the past. The configuration appearing in all panels with the same orbital parameters is marked by a fill circle.

5. Application to HD 209458b

Here we compare the results of this paper with those of Mardling (2007) for the example of HD 209458b (Table 1). As in Mardling (2007), we first assume that the non-zero eccentricity of this planet is caused by a m₂ = 0.1   M_J companion at a₂ = 0.4 AU with an eccentricity e₂ = 0.4. For this study, the eccentricities are large and modify the frequencies g_k given by the analytical expression of the matrix A_tot (10). Thus, we chose to compute the matrix A_tot using frequency analysis for a numerical integration of the system without dissipation, exact in eccentricity, and expanded up to the fourth order in the ratio of the semi-major axes (e.g. Mardling & Lin 2002; Laskar & Boué 2010). To first order, the eccentricity variables z₁ and z₂ are linear combinations of two eigenmodes u₁ and u₂ (27) (31)where S is given by (21). With Q₁ = 10⁵ and ω₁ = n₁, the two damping timescales (29) are  Myr and  Gyr. With an age estimate of 5.5 Gyr for this system (Burrows et al. 2007), the first eigenmode should be damped and the modulus of the second should remain almost constant. In consequence, both eccentricity variables should be proportional to u₂. Their modulus should thus be constant and verify (30), or equivalently e₁ ≈ e₂   S₁₂/S₂₂ = 0.0025. In Mardling (2007, Fig. 3), this value is larger, e₁ = 0.012. The difference comes from the tidal deformation of the planet that leads to the coefficient in Eq. (14). This was not taken into account in Mardling (2007), and it accelerates the precession of the periapse of HD 209458b by a factor of 6.6 (see Fig. 2a). In Fig. 2a, the initial Q-value of the planet is set artificially to 100 to enable a direct comparison with Fig. 3 of Mardling (2007). As said by Mardling (2007), and shown in that paper, the Q-value affects the damping timescales but neither the precession frequencies nor the eccentricity amplitudes. However, with a higher precession frequency, the matrix A_tot is closer to a diagonal matrix. The two planets are more weakly coupled and the ratio S₁₂/S₂₂ (31), equal to the final eccentricity ratio e₁/e₂, is smaller.

5.1. Residual eccentricity

One way to recover the final eccentricity of HD 209458b is to increase the mass of the companion up to m₂ = 0.608   M_J (Fig. 2b). Here, our aim is not to explain the large eccentricity of HD 209458b, but simply to illustrate the results of Sect. 3.2.

As the precession of the periastron of the inner planet is faster than in Mardling (2007, Fig. 3), we decreased the initial Q-value to 15.15 to accelerate the damping and to obtain an evolution with the same g₁/γ₁ ratio as in Mardling (2007) (Fig. 2c). As said before, this does not change the final eccentricity, but more clearly illustrates the damping of the first mode with frequency g₁ = 0.14 deg/yr.

After the damping of the first eigenmode, the eccentricities do not oscillate because only one eigenmode with a non-zero amplitude remains. Both eccentricity variables z₁ and z₂ describe a circle in the complex plane with the same frequency g₂. However, if a third planet is added to the system, a new eigenmode appears with a frequency g₃. The eccentricities then oscillate (Fig. 2d). We note that the relative inclination between the planets can also generate another eigenmode and cause the eccentricities to oscillate (Mardling 2010). However, in the linear approach eccentricities and inclinations are decoupled. It is thus necessary to have large eccentricities or inclinations to ensure significant oscillations.

5.2. Possibility of a planet companion

We next attempted to determine which companion parameters can lead to an eccentricity e₁ = 0.01 for HD 209458b. As the system contains two planets, their eccentricities are at most a combination of two eigenmodes. However, since γ^-1 = 46 Myr is less than the age of the system (5.5 Gyr), at least one of the eigenmode is damped. However both eigenmodes cannot have zero amplitude, otherwise the two orbits would be circular. Thus, we assumed that there remains only a single eigenmode with non-zero amplitude, the one with the longer damping timescale. Given a semi-major axis a₂ and a mass m₂, the eccentricity of the companion was then obtained using Eq. (30) to first order. In practice, we numerically integrated the system without dissipation, and found the value of the eccentricity e₂ that cancels the amplitude of the rapidly damped eigenmode. The results are shown with solid curves in Figs. 3a and b. After the current eccentricity e₂ had been given, the initial value (5.5 Gyr ago) was estimated by assuming an exponential decay with a damping factor given by (29) (see the dotted curves Figs. 3a and b). The frequencies g_k and the coefficients A_kk were obtained numerically using frequency analysis. Parameters leading to initial eccentricities larger than 1 were excluded; these correspond to the grey regions in Fig. 3. Although the planets were more weakly coupled than in Mardling (2007), there was still a large range of initial conditions leading to a state compatible with e₁ = 0.01. However, the stellar reflex velocity produced by the companion at periastron (Figs. 3c and d) is above the detectability threshold of about 3 m s^-1. For example, with a₂ = 0.25 AU, and m₂ = 0.05   M_J, the current eccentricity is e₂ = 0.34 and the maximal stellar reflex velocity  m s^-1. It thus seems that the existence of such a planet cannot be assumed to explain the observed eccentricity. Indeed, observations do not strongly constrain the eccentricity of HD 209458b and a circular orbit is not ruled out (Laughlin et al. 2005).

6. Orbit fitting: the HD 10180 case

The analysis of the radial velocity measurement of HD 10180 revealed the potential existence of seven planets in this system (Lovis et al. 2011). The innermost planet, HD 10180b, is a terrestrial planet (m_bsini = 1.35   M_⊕) with a period of  ≈ 1.177 days and a semi-major axis a_b = 0.0223 AU. The planet is thus subject to strong tidal interactions with the central star. During the first orbital fit (Lovis et al. 2011, Table 3), it was thus assumed that its eccentricity had been damped to very small values, and its value was fixed to e_b = 0. Nevertheless, when the system was then numerically integrated over 250 kyr (Fig. 4 (red curve)), owing to secular interactions with the other planets, e_b was found to increase very rapidly to high values, reaching nearly 0.9.

When general relativity (GR) is included in the numerical integration, the main effect is to increase the diagonal terms of the secular matrix (Eq. (11)). As a result, the secular variations in e_b are much smaller (Fig. 4 (green curve)), but still to reach 0.2.

The strategy adopted for the final fit of Lovis et al. (2011) was to include in the fit the constraint that the planetary system that is searched for is the result of the tidal evolution, as described in Sect. 3.2. As the planet has a mass comparable to the Earth, it could be assumed to be terrestrial, with a dissipation factor of the same order of magnitude as (or larger than) Mars (k₂/Q = 0.0015), which is the smallest value among the terrestrial planets in the Solar System. The damping factors e ^− γ_kt can thus be computed using Eq. (25) for all proper modes u_k (Lovis et al. 2011, Table 5). The resulting dampings of the amplitudes of the proper modes u_k are given in Fig. 5.

Fig. 4 Evolution of the eccentricity of planet HD 10180b over 250 kyr starting with eb = 0 at t = 0 (present time) for three different models: In red, the numerical integration is purely Newtonian and do not take into account general relativity (GR). In green, GR is taken into account in the integration. In blue, GR is taken into account and the fit is made with the tidal dissipation constraint (32).

Fig. 5 Tidal evolution of the amplitude of the proper modes |u1| (red), |u2| (green), |u3| (blue), and |u4| (pink) resulting from the tidal dissipation on planet HD 10180b with k2/Q = 0.0015 (Lovis et al. 2011).

From this computation, as the age of the system is estimated to be about 4 Gyr, it can be seen that the first two proper modes amplitudes u₁ and u₂ are certainly reduced to very small values. If the damping factor k₂/Q were ten times smaller, the conclusion would be nearly the same, as the only change in Fig. 5 would be a change in the timescale of the figure, the units being now 10 Gyr instead of 1 Gyr.

To include the above constraint on the tidal damping in the fit, one can then add to the χ² minimization the additional term (32)where R is an empirical constant that needs to be set to a value that will equilibrate the damping constraint with respect to the value of the χ² in the absence of constraint. After various trials, R = 350 was used in Lovis et al. (2011).

The computation of the amplitude of the proper modes |u_k| during the fit is made iteratively. Once a first orbital solution is obtained, the Lagrange-Laplace linear system (Eq. (1)) is computed and thus the matrix S₀ of transformation to proper modes (Eq. (18)). For a given initial condition (z_k) obtained through the fit, the proper modes u_k are computed with (33)and the additional contribution (32) can then be computed in the fitting process. Practically, in an iterative fit that takes into account the Newtonian interactions, the transformation matrix just needs to be computed once, or twice, if one wants to recompute the matrix when the convergence to a final solution is obtained. In Lovis et al. (2011), the final values were u₁ = 0.0017, u₂ = 0.044 for R = 350, with a final , very close to the residuals obtained in the absence of any constraint ().

In this constrained solution, the initial value of e_b is still 0, but the secular change caused by planetary interactions is much smaller (Fig. 4 (blue curve)), ensuring a more stable behavior to the system.

7. Conclusion

We have presented here in full detail the secular theory with tidal dissipation that was outlined in Lovis et al. (2011) for the system HD 10180. The use of Lagrange-Laplace linear theory can very easily include the linear contribution of tidal dissipation and provide an intuitive background for studying multi-planetary systems when one or several planets are close to their central star and subject to tidal damping. Although we have limited ourselves here to the study of the planar case, this formalism can be easily extended to mutually inclined systems.

For the system HD 209458, we could retrieve globally the results of Mardling (2007), although we have found that a companion with mass m₂ = 0.1   M_J, a₂ = 0.4 AU, and e₂ = 0.4 will not lead to e₁ = 0.012, but to a much smaller value of e₁ = 0.0025. This is due to the additional tides generated by the star at the planet (Eq. (14)) in the secular equations (Eq. (10)).

We have examined other configurations that could lead to a final eccentricity of e₁ ≥ 0.01 for HD 209458b, but our conclusions are negative, as we found that a potential companion, massive enough to lead to a final eccentricity e₁ ≥ 0.01 leads to sufficiently large stellar motion that it should already have been detected, assuming a detectability threshold of 3 m s. Our

conclusion is thus that the most probable outcome is that the actual eccentricity of HD 209458b has a value much smaller than 0.01.

Acknowledgments

This work has been supported by PNP-CNRS, by CS of Paris Observatory, by the European Research Council/European Community under the FP7 through a Starting Grant, as well as in the form of grant reference PTDC/CTE-AST/098528/2008, funded by Fundação para a Ciência e a Tecnologia (FCT), Portugal.

References

All Tables

All Figures

Fig. 1 a) The complex variable z = eexp(iϖ) provides the orientation and eccentricity of the orbit in the plane. Alternatively, one can use the rectangular real coordinates k = ecosϖ, h = esinϖ, as in Murray & Dermott (1999). b) If the Lagrange-Laplace matrix A were diagonal, the eccentricity of all planets would be constant with uniform precession. c) As A is not diagonal, one needs to make a linear transformation to transform zk into proper modes uk. The proper modes amplitude are then constant and uk precess uniformly with frequency gk. d) The full solution zj = ejexp(ϖj) are linear combinations of the proper modes uk. In the figure, with two proper modes (i.e. two planets and n = 2). The red curve represents the evolution of the eccentricity and longitude of perihelia of Jupiter over 200 000 years under the perturbation of Saturn.

In the text

Fig. 2 Tidal effects on the eccentricity of HD 209458b with one or two companions. a) The companion is a 0.1   MJ planet at 0.4 AU with e2 = 0.4 as in Mardling (2007). The blue curve was obtained without considering the conservative effect of the tides ( in Eq. (14)). The red curve in a), and all the evolutions presented in the subfigures b)–d) take into account this effect. b) The mass of the companion is set to 0.608   MJ to recover the final excentricity of Mardling’s simulation. The red curve is the result of a numerical integration of the full secular equations that are exact in eccentricity. The green one is the analytical solution of the linearized problem. c) Same as b) except for the initial Q-value of HD 209458b, which is set to 15.15 to enable the visualisation of both the damping and the oscillation of the eccentricity. d) Same as c) with an additional 0.1   MJ companion at 1.0 AU with e3 = 0.1.

In the text

Fig. 3 Eccentricity e2 of the hypothetical companion of HD 209458b with e1 = 0.01 assuming that the eigenmode with the shortest dissipation timescale is damped (black curves in panel a) and b)). a) The mass of the companion is fixed to m2 = 0.2   MJ. Negative values of e2 correspond to Δϖ = 180deg while positive ones mean Δϖ = 0deg. The dotted line is the eccentricity that the companion would have had 5.5 Gyr ago assuming a dissipation factor computed with (29). b) Same as a) for different masses m2 while the semi-major axis is fixed and set to a2 = 0.4 AU. c) Stellar reflex velocity due to the companion at periastron with the eccentricity of a). d) Idem for the eccentricity of b). In grey regions, the eccentricity of the companion should have been larger than 1 in the past. The configuration appearing in all panels with the same orbital parameters is marked by a fill circle.

In the text

Fig. 4 Evolution of the eccentricity of planet HD 10180b over 250 kyr starting with eb = 0 at t = 0 (present time) for three different models: In red, the numerical integration is purely Newtonian and do not take into account general relativity (GR). In green, GR is taken into account in the integration. In blue, GR is taken into account and the fit is made with the tidal dissipation constraint (32).

In the text

	Fig. 5 Tidal evolution of the amplitude of the proper modes |u1| (red), |u2| (green), |u3| (blue), and |u4| (pink) resulting from the tidal dissipation on planet HD 10180b with k2/Q = 0.0015 (Lovis et al. 2011).
In the text