PUMPING THE ECCENTRICITY OF EXOPLANETS BY TIDAL EFFECT

Close-in exoplanets, as for Mercury, Venus, and the majority of the natural satellites in the solar system, are supposed to undergo significant tidal interactions, resulting in their spins and orbits being slowly modified. The ultimate stage for tidal evolution is the synchronization of the spin and the circularization of the orbit (e.g., Correia 2009). The spin evolves in a shorter timescale than the orbit, so long-term studies on the tidal evolution of exoplanets usually assume that their rotation is synchronously locked, and therefore limit the evolution to the orbits. However, these two kinds of evolution cannot be dissociated because the total angular momentum must be conserved. Synchronous rotation can only occur when the eccentricity is very close to zero. Otherwise, the rotation rate tends to be locked with the orbital speed at the periapsis, because tidal effects are stronger when the two bodies are closer to each other. In addition, in the presence of a companion body, the eccentricity undergoes oscillations (e.g., Mardling 2007), and the rotation rate of the planet shows variations that follow the eccentricity (Correia & Laskar 2004). As a consequence, some unexpected behaviors can be observed, such as a secular increase of the eccentricity. In this Letter, we provide a simple averaged model for the orbital and spin evolution of an exoplanet with a companion (Section 2) and apply it to the HD 117618 planetary system (Section 3). We then give an explanation for the eccentricity pumping (Section 4) and derive some conclusions (Section 5).

We consider here a system consisting of a central star of mass m₀, an inner planet of mass m₁, and an outer companion of mass m₂. We use Jacobi canonical coordinates, with $\mathbf {r_1}$ being the position of m₁ relative to m₀, and $\mathbf {r_2}$ the position of m₂ relative to the center of mass of m₁ and m₀. We further assume that $|\mathbf {r_1}| \ll |\mathbf {r_2}|$, that the system is coplanar, and that the obliquity of the planet is zero. The inner planet is considered an oblate ellipsoid with gravity field coefficient given by J₂, rotating about the axis of maximal inertia, with rotation rate ω, such that (e.g., Lambeck 1988)

$$\begin{equation} J_2 = k_2 \frac{\omega _{}^2 R_{}^3}{3 G m_1}, \end{equation} \tag{ 1 }$$

where G is the gravitational constant, R is the radius of the planet, and k₂ is the second Love number for potential.

Since we are interested in the secular behavior of the system, we average the motion equations over the mean anomalies of both orbits. The averaged potential, quadrupole level for the spin (e.g., Correia & Laskar 2010a), octopole level for the orbits (e.g., Lee & Peale 2003; Laskar & Boué 2010), and with general relativity corrections (e.g., Touma et al. 2009), is given by

$$\begin{eqnarray} U & =& {-}C_0 \left(1-e_1^2\right)^{-1/2} - C_1 \left(1-e_1^2\right)^{-3/2} \cr & & - C_2 \frac{\left(1 + \frac{3}{2} e_1^2\right)}{\left(1-e_2^2\right)^{3/2}} + C_3 \frac{e_1e_2\left(1 + \frac{3}{4} e_1^2\right)}{\left(1-e_2^2\right)^{5/2}} \cos \varpi , \end{eqnarray} \tag{ 2 }$$

with

$$\begin{equation} C_0 = \frac{3 \beta _1G^2 (m_0+ m_1)^2}{a_1^2 c^2}, \quad C_1 = \frac{G m_0m_1J_2 R_{}^2}{2 a_1^3} , \end{equation} \tag{ 3 }$$

$$\begin{equation} C_2 = \frac{G \beta _1m_2a_1^2}{4 a_2^3 }, \quad C_3 = \frac{15 G \beta _1m_2a_1^3}{16 a_2^4} \frac{(m_0-m_1)}{m_0+ m_1} , \end{equation} \tag{ 4 }$$

where a_i is the semimajor axis (that can also be expressed using the mean motion n_i), e_i is the eccentricity, and ϖ = ϖ₁ − ϖ₂ is the difference between the longitudes of the periastron, ϖ_i. We also have β₁ = m₀m₁/(m₀ + m₁) and β₂ = (m₀ + m₁)m₂/(m₀ + m₁ + m₂).

The contributions to the orbits are easily obtained using the Lagrange planetary equations (e.g., Murray & Dermott 1999):

$$\begin{equation} \dot{e}_i= \frac{\sqrt{1-e_i^2}}{\beta _in_ia_i^2 e_i} \frac{\partial U}{\partial \varpi _i}, \quad \dot{\varpi }_i= - \frac{\sqrt{1-e_i^2}}{\beta _in_ia_i^2 e_i} \frac{\partial U}{\partial e_i} . \end{equation} \tag{ 5 }$$

Thus,

$$\begin{equation} \dot{e}_1= - \nu _{31} \frac{e_2\left(1 + 3/4 e_1^2\right) \sqrt{1-e_1^2}}{\left(1-e_2^2\right)^{5/2}} \sin \varpi , \end{equation} \tag{ 6 }$$

$$\begin{equation} \dot{e}_2= \nu _{32} \frac{e_1\left(1 + 3\big/4 e_1^2\right)}{\left(1-e_2^2\right)^{2}} \sin \varpi , \end{equation} \tag{ 7 }$$

and

$$\begin{eqnarray} \dot{\varpi }&=& \frac{\nu _0}{ \left(1-e_1^2\right)} + \frac{\nu _1 \, x^2}{\big(1-e_1^2\big)^2}\nonumber\\ &&+ \nu _{21} \frac{\sqrt{1-e_1^2}}{ \left(1-e_2^2\right)^{3/2}} - \nu _{22} \frac{ \left(1+\frac{3}{2} e_1^2\right)}{ (1-e_2^2)^2}\nonumber\\ &&- \nu _{31} \frac{e_2\sqrt{1-e_1^2} \left(1+\frac{9}{4} e_1^2\right)}{e_1\big(1-e_2^2\big)^{5/2}} \cos \varpi\nonumber\\ &&+ \nu _{32} \frac{e_1\left(1+\frac{3}{4} e_1^2\right) \big(1+4 e_2^2\big)}{e_2\left(1-e_2^2\right)^3} \cos \varpi, \end{eqnarray} \tag{ 8 }$$

where x = ω/n₁, and the constant frequencies

$$\begin{equation} \nu _0= 3 n_1\!\left(\!\frac{n_1a_1}{c}\!\right)^{\!\!2} , \end{equation} \tag{ 9 }$$

$$\begin{equation} \nu _1 = n_1\frac{k_2}{2} \frac{m_0+ m_1}{m_1} \!\left(\!\frac{R_{}}{a_1}\!\right)^{\!\!5} , \end{equation} \tag{ 10 }$$

$$\begin{equation} \nu _{21} = n_1\frac{3}{4} \frac{m_2}{m_0+ m_1} \!\left(\!\frac{a_1}{a_2}\!\right)^{\!\!3} , \end{equation} \tag{ 11 }$$

$$\begin{equation} \nu _{22} = n_2\frac{3}{4} \frac{m_0m_1}{(m_0+m_1)^2} \!\left(\!\frac{a_1}{a_2}\!\right)^{\!\!2} , \end{equation} \tag{ 12 }$$

$$\begin{equation} \nu _{31} = n_1\frac{15}{16} \frac{m_2}{m_0+m_1} \frac{m_0-m_1}{m_0+m_1} \!\left(\!\frac{a_1}{a_2}\!\right)^{\!\!4} , \end{equation} \tag{ 13 }$$

$$\begin{equation} \nu _{32} = n_2\frac{15}{16} \frac{m_0m_1}{(m_0+m_1)^2} \frac{m_0-m_1}{m_0+m_1} \!\left(\!\frac{a_1}{a_2}\!\right)^{\!\!3} . \end{equation} \tag{ 14 }$$

Note that the variations in e₁ and e₂ (Equations (6) and (7)) are related by the conservation of the total angular momentum (after dividing by m₁n₂a²₂):

$$\begin{equation} \xi \frac{\omega _{}}{n_2} \!\left(\!\frac{R_{}}{a_2}\!\right)^{\!\!2} + \frac{\beta _1}{m_1} \frac{n_1}{n_2} \!\left(\!\frac{a_1}{a_2}\!\right)^{\!\!2} \!\!\sqrt{1-e_1^2} + \frac{\beta _2}{m_1} \!\sqrt{1-e_2^2} = Cte, \end{equation} \tag{ 15 }$$

where ξ is a structure coefficient. The conservative system (Equation (2)) can thus be reduced to one degree of freedom.

In our model, we additionally consider tidal dissipation raised by the central star on the inner planet. The dissipation of the mechanical energy of tides in the planet's interior is responsible for a time delay Δt between the initial perturbation and the maximal deformation. As the rheology of planets is badly known, the exact dependence of Δt on the tidal frequency is unknown. Several models exist (for a review see Correia et al. 2003; Efroimsky & Williams 2009), but for simplicity we adopt here a model with constant Δt, which can be made linear (Singer 1968; Mignard 1979). The contributions to the equations of motion are given by (e.g., Correia 2009)

$$\begin{equation} \frac{\dot{\omega }_{}}{n_1} = - K_{}\left(f_1(e_1) x_{}- f_2(e_1) \right) , \end{equation} \tag{ 16 }$$

$$\begin{equation} \frac{\dot{a}_1}{a_1} = 2 K^{\prime }(f_2(e_1) x - f_3(e_1)) , \end{equation} \tag{ 17 }$$

$$\begin{equation} \dot{e}_1= 9 K^{\prime }_{}\left(\frac{11}{18} f_4(e_1) x - f_5(e_1) \right) e_1, \end{equation} \tag{ 18 }$$

where

$$\begin{equation} K_{}= n_1\frac{3 k_2}{\xi Q} \frac{m_0\beta _1}{m_1^2} \!\left(\!\frac{R_{}}{a_1}\!\right)^{\!\!3} , \quad K^{\prime }_{}= \frac{K_{}}{1/\xi } \frac{m_1}{\beta _1} \!\left(\!\frac{R_{}}{a_1}\!\right)^{\!\!2} , \end{equation} \tag{ 19 }$$

Q⁻¹ = n₁Δt, and f₁(e) = (1 + 3e² + 3e⁴/8)/(1 − e²)^9/2, f₂(e) = (1 + 15e²/2 + 45e⁴/8 + 5e⁶/16)/(1 − e²)⁶, f₃(e) = (1 + 31e²/2 + 255e⁴/8 + 185e⁶/16 + 25e⁸/64)/(1 − e²)^15/2, f₄(e) = (1 + 3e²/2 + e⁴/8)/(1 − e²)⁵, f₅(e) = (1 + 15e²/4 + 15e⁴/8 + 5e⁶/64)/(1 − e²)^13/2.

We neglect the effect of tides over the longitude of the periastron as well as the flattening of the central star. Their effect is only to add a small supplementary frequency to ϖ₁, similar to the contributions from general relativity (for a complete model see Correia et al. 2011).

Under the effect of tides alone, the equilibrium rotation rate, obtained when $\dot{\omega }_{}= 0$, is attained for (Equation (16))

$$\begin{equation} \frac{\omega _{}}{n_1} = \frac{f_2(e_1)}{f_1(e_1)} = 1 + 6 e_1^2 + {\cal O}\left(e_1^4\right) . \end{equation} \tag{ 20 }$$

Usually K' ≪ K, so tidal effects modify the rotation rate much faster than the orbit. It is thus tempting to replace the equilibrium rotation in expressions (17) and (18). With this simplification, one always obtains negative contributions for $\dot{a}_1$ and $\dot{e}_1$ (Correia 2009),

$$\begin{equation} \frac{\dot{a}_1}{a_1} = - 7 K_{}^{\prime } \, f_6(e_1) e_1^2 , \end{equation} \tag{ 21 }$$

$$\begin{equation} \dot{e}_1= - \frac{7}{2} K_{}^{\prime } f_6(e_1) \left(1-e_1^2\right) e_1, \end{equation} \tag{ 22 }$$

with f₆(e) = (1 + 45e²/14 + 8e⁴ + 685e⁶/224 + 255e⁸/448 + 25e¹⁰/1792)(1 − e²)^−15/2/(1 + 3e² + 3e⁴/8). Thus, the semimajor axis and the eccentricity can only decrease until the orbit of the planet becomes circular (Figure 1(a)). However, planet–planet interactions can produce eccentricity oscillations with a period shorter than, or comparable to, the damping timescale of the spin. In that case, expression (20) is not satisfied and multi-planetary systems may show non-intuitive eccentricity evolutions, such as eccentricity pumping of the inner orbit (e₁ increases while e₂ decreases).

Zoom In Zoom Out Reset image size

As an illustration of the eccentricity pumping, we apply our model to the HD 117618 system. This Sun-like star (m₀ ≈ M_☉) has been reported to host a single Saturn-like planet on a eccentric orbit (Butler et al. 2006). The residuals of the best-fit solution to the observational data are 5.5 m s⁻¹, so we assume that any additional companion with a Doppler shift semi-amplitude smaller than this value is presently undetected, that is, any planet with m₂ < 0.2 M_J and a₂ > 1.4 AU.

In our simulations we adopt for the observed planet the same geophysical parameters as for Saturn, R = R_Sat, k₂ = 1/2, ξ = 1/5, and a dissipation time lag Δt = 200 s (which is equivalent to Q ≈ 3 × 10⁴). Since the semimajor axis of the planet undergoes tidal dissipation, its value was certainly larger when the system formed. We then adopt a₁ = 0.25 AU as the initial value for all simulations. The initial eccentricity is chosen such that e₁ ≈ 0.4 when a₁ = 0.175, the present observed values (Table 1). We further assume initial ϖ = 180°, and 2π/ω = 50 day, which quickly evolves near the equilibrium rotation (Equation (20)).

In the absence of a companion, the eccentricity and the semimajor axis are damped following an exponential decay (Correia & Laskar 2010b), and the present configuration is attained after 1 Gyr (Figure 1(a)). At present the observed eccentricity would be around 0.1, and we still needed to explain the high initial value near 0.7.

We now add a companion to the system with m₂ = 0.2 M_J, a₂ = 1.8 AU, and e₂ = 0.4, and set e₁ = 0.3. At first, we only consider dissipation in the spin (Equation (16)) and neglect its effect on the orbit (Equations (17) and (18)) in order to highlight the eccentricity pumping (Figure 1(b)). We then clearly observe this effect, with the eccentricity of the inner planet rising up to 0.7. We also observe that the eccentricity of the outer planet is simultaneously damped because of the conservation of the total angular momentum (Equation (15)).

Orbital and spin evolution cannot be dissociated, so we then integrate the full system (Figure 1(c)). We observe that the initial behavior of the system is identical to the situation without dissipation on the orbit (Figure 1(b)). However, as the eccentricity increases, the inner planet comes closer to the star at periastron, and tidal effects on the orbit become stronger. As a consequence, the semimajor axis decreases and the damping effect on the eccentricity (Equation (18)) overrides the pumping drift. The system ultimately evolves into a circular orbit. The present configuration is attained around 4 Gyr of evolution, which is compatible with the present estimated age of the star. The pumping effect is then responsible for a delay in the final evolution of planetary systems and may explain the high values observed for some of them (Table 1).

Finally, we repeat the integration of the full system, but with a smaller-mass companion m₂ = 0.15 M_J at a₂ = 1.4 AU. The companion eccentricity is still e₂ = 0.4, but the inner planet now begins with e₁ = 0.15 (Figure 1(d)). The initial evolution is still similar to the previous simulation (Figure 1(c)), except that the eccentricity oscillations of both planets are higher because the orbits are closer. Since the companion mass is smaller, its eccentricity also decreases more than before and reaches zero around 3 Gyr. At this stage, the angle ϖ stops circulating and begins librating around 180°. The two orbits are then tightly coupled and evolve together, showing an identical behavior to close-in planets with an orbital period of a few days (e.g., Mardling 2007; Laskar et al. 2011). As a consequence, the evolution timescale is much longer, allowing the inner planet to maintain high eccentricity for longer periods of time (Figure 1(d)). The present eccentricity is only observed after 6 Gyr.

In order to understand the unexpected behavior of the eccentricity during the initial stages of the evolution, we can perform some simplifications in the equations of motion without loss of generality (Section 2). We can neglect tidal effects on orbital quantities (Equations (17) and (18)), which is justified since K' ≪ K (Equation (19)). The only contribution of tides is then on the rotation rate (Equation (16)). The semimajor axis and the mean motion are thus constant, and the eccentricity only varies due to the gravitational perturbations (Equation (6)). In addition, we linearize the set of equations of motion in the vicinity of the averaged values of x, e₁, and e₂. Let x = x₀ + δx, where x₀ is the solution of Equation (20), e₁ = e₁₀ + δe₁, and e₂ = e₂₀ + δe₂. In the following, δe₂ is expressed as a function of δx and δe₁ using the conservation of the angular momentum (Equation (15)). Then, the set of equations of motion (6), (8), and (16) reduces to

$$\begin{equation} \delta \dot{e}_1= - A \sin \varpi , \end{equation} \tag{ 23 }$$

$$\begin{equation} \dot{\varpi }= g + g_x \delta x + g_e \delta e_1, \end{equation} \tag{ 24 }$$

$$\begin{equation} \delta \dot{x} = -\nu _x \delta x + \nu _e \delta e_1, \end{equation} \tag{ 25 }$$

with

$$\begin{equation} A = \nu _{31} \frac{e_{20} \big(1 + 3/4 e_{10}^2\big) \sqrt{1-e_{10}^2}}{\left(1-e_{20}^2\right)^{5/2}} , \end{equation} \tag{ 26 }$$

$$\begin{eqnarray} g &=& \frac{\nu _0}{\left(1-e_{10}^2\right)} + \frac{\nu _1 x_0^2}{\left(1-e_{10}^2\right)^2} \cr && + \nu _{21} \frac{\sqrt{1-e_{10}^2}}{\left(1-e_{20}^2\right)^{3/2}} - \nu _{22} \frac{\left(1+3 e_{10}^2\big/2\right)}{\left(1-e_{20}^2\right)^2} , \end{eqnarray} \tag{ 27 }$$

$$\begin{equation} g_x = \nu _1 \frac{2 x_0}{\left(1-e_{10}^2\right)^2} , \end{equation} \tag{ 28 }$$

$$\begin{eqnarray} g_e &=& \nu _0 \frac{2 e_{10}}{\left(1-e_{10}^2\right)^2} + \nu _1 \frac{4 x_0^2 e_{10}}{\left(1-e_{10}^2\right)^3} - \nu _{22} \frac{3 e_{10}}{\left(1-e_{20}^2\right)^2} \cr && - \nu _{21} \frac{e_{10}}{\sqrt{1-e_{10}^2} \left(1-e_{20}^2\right)^{3/2}} , \end{eqnarray} \tag{ 29 }$$

$$\begin{equation} \nu _x = K_{}f_1(e_{10}) , \end{equation} \tag{ 30 }$$

$$\begin{equation} \nu _e = -K_{}(f^{\prime }_1 (e_{10}) x_0 - f^{\prime }_2(e_{10})) , \end{equation} \tag{ 31 }$$

where f'₁(e) = 15(e + 3e³/2 + e⁵/8)/(1 − e²)^11/2 and f'₂(e) = 3(9e + 65e³/2 + 125e⁵/8 + 5e⁷/8)/(1 − e²)⁷. We neglected the octopole terms since ν_3i ≪ ν_2i, the contributions from δe₂, and assumed that e_i0 ≠ 0.

At first order, the precession of the periastron is constant, $\dot{\varpi }\simeq g$, and the eccentricity is simply given from expression (23) as

$$\begin{equation} \delta e_1= \Delta e\cos (g t + \varpi _0) , \end{equation} \tag{ 32 }$$

where Δe = A/g. That is, the eccentricity e₁ presents periodic variations around an equilibrium value e₁₀, with amplitude Δe and frequency g. Since g_xδx, g_eδe₁ ≪ g, the above solution for the eccentricity can be adopted as the zeroth-order solution of the system of Equations (23)–(25). With this approximation, the equation of motion (25) of δx becomes that of a driven harmonic oscillator whose steady state solution is

$$\begin{equation} \delta x = \Delta x\cos (g t + \varpi _0 - \phi) , \end{equation} \tag{ 33 }$$

with $\Delta x= \nu _e \Delta e/ \sqrt{\nu _x^2 + g^2 }$ and $\sin \phi = g/\sqrt{\nu _x^2 + g^2}$. The rotation rate thus presents an oscillation identical to the eccentricity (Equation (32)), but with smaller amplitude and delayed by an angle ϕ (see Correia 2011). Using the above expression in Equation (24) and integrating gives for the periastron:

$$\begin{equation} \varpi = g t + \varpi _0 + \frac{g_x}{g} \Delta x\sin (g t + \varpi _0 - \phi) + \frac{g_e}{g} \Delta e\sin (g t + \varpi _0) . \end{equation} \tag{ 34 }$$

Finally, substituting into Equation (23) and using the approximation g_xΔx, g_eΔe ≪ g gives

$$\begin{eqnarray} \delta \dot{e}_1& \approx & {-}A \sin (g t + \varpi _0) - \frac{g_e A}{g} \Delta e\sin (g t + \varpi _0) \cos (g t + \varpi _0) \cr & & - \frac{g_x A}{g} \Delta x\sin (g t + \varpi _0 - \phi) \cos (g t + \varpi _0) , \end{eqnarray} \tag{ 35 }$$

or, combining the two products of periodic functions,

$$\begin{eqnarray} \delta \dot{e}_1& = & {-}A \sin (g t + \varpi _0) - \frac{g_x A}{2 g} \Delta x\sin (2 g t + 2 \varpi _0 - \phi) \cr & & - \frac{g_e A}{2 g} \Delta e\sin (2 g t + 2 \varpi _0) + \frac{g_x A}{2 g} \Delta x\sin \phi . \end{eqnarray} \tag{ 36 }$$

The two middle terms in the above equation can be neglected since they are periodic and have a very small amplitude (g_xΔx, g_eΔe ≪ g). However, the last term in sin ϕ is constant and it adds a small drift to the eccentricity,

$$\begin{equation} \langle e_1\rangle = \frac{\nu _e g_x A^2}{2 g \left(\nu _x^2 + g^2\right)} t . \end{equation} \tag{ 37 }$$

The drift is maximized for g ∼ ν_x, and it vanishes for weak dissipation (Δx → 0), and also for strong dissipation (ϕ → 0). Note that the phase lag ϕ between the eccentricity (Equation (32)) and the rotation variations (Equation (33)) is essential to get a drift on the eccentricity. That is why the eccentricity pumping was never observed in previous studies that did not take into account the spin evolution.

The major difference when we consider the full non-linearized problem is that the drift (Equation (37)) cannot grow indefinitely. Indeed, when the eccentricity reaches high values, the drift vanishes (Figure 1(b)). Moreover, tidal effects are also enhanced for high eccentricities and counterbalance the drift (Equation (18)). As a consequence, the drift on the eccentricity is never permanent, although it can last for the age of the system (Figures 1(c) and (d)).

In order to observe the pumping effect, the eccentricity should not be damped, while the damping timescale of the spin of the planet should be of the order of the period of eccentricity oscillations. This is valid for gaseous planets roughly within 0.1 < a₁ < 0.3. About half of the planets in this range present eccentricities higher than 0.3 (Table 1). This is somehow unexpected, since the timescale for damping the eccentricity (τ ∼ 1/K') is shorter than the age of those systems (Table 1, Figure 1(a)). Thus, unless the initial eccentricity of those planets was extremely high, one may suspect the existence of undetected companions that help the eccentricity to maintain present values.

In this Letter we have shown that, under some particular initial conditions, orbital and spin evolution cannot be dissociated. Indeed, counterintuitive behaviors can be observed, such as the secular augmentation of the eccentricity. This effect can last over long timescales and may explain the high eccentricities observed for moderate close-in planets.

The variation in the flattening of the inner planet due to rotation (Equation (1)) is a key element to pump the eccentricity by means of g_x (Equation (37)). We have considered an instantaneous response to the rotation in J₂, but a time delay until the maximum deformation is reached is to be expected. This may increase the phase lag between the eccentricity forcing (Equation (32)) and the precession angle oscillations (Equation (34)). The result is that the drift effect on the eccentricity can be even more pronounced than the one presented here.

We acknowledge support from the PNP-CNRS, the FCT (grant PTDC/CTE-AST/098528/2008), and the European Research Council.