Issue 
A&A
Volume 635, March 2020



Article Number  A37  
Number of page(s)  8  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201936967  
Published online  03 March 2020 
Why do warm Neptunes present nonzero eccentricity?
^{1}
CFisUC, Department of Physics, University of Coimbra,
3004516
Coimbra,
Portugal
email: acor@uc.pt
^{2}
ASD, IMCCE, Observatoire de Paris, PSL Université,
77 Av. DenfertRochereau,
75014
Paris,
France
^{3}
Observatoire de l’Université de Genève,
51 chemin des Maillettes,
1290
Sauverny,
Switzerland
Received:
21
October
2019
Accepted:
17
January
2020
Most Neptunemass planets in closein orbits (orbital periods less than a few days) present nonzero eccentricity, typically around 0.15. This is somehow unexpected, as these planets undergo strong tidal dissipation that should circularize their orbits in a timescale shorter than the age of the system. In this paper we discuss some mechanisms that can oppose to bodily tides, namely, thermal atmospheric tides, evaporation of the atmosphere, and excitation from a distant companion. In the first two cases, the eccentricity can increase consistently, while in the last one, the eccentricity can only be excited for a limited amount of time (that may nevertheless exceed the age of the system). We show the limitations of these different mechanisms and how some of them could, depending on specific properties of the observed planetary systems, account for their presently observed eccentricities.
Key words: celestial mechanics / planetstar interactions / planets and satellites: dynamical evolution and stability / planets and satellites: atmospheres
© ESO 2020
1 Introduction
Planets with orbital periods P_{orb} ≲ 10 day undergo strong tidal interactions with the parent star (e.g., MacDonald 1964). Bodily (or gravitational) tides are raised on the planet by the star because of the gravitational gradient across the planet. Since planets are not perfectly rigid, there is a distortion that gives rise to a tidal bulge. The dissipation of the mechanical energy inside the planet introduces a delay, and, hence, a phase shift between the initial perturbation and the maximal tidal deformation. As a consequence, the star exerts a torque on the tidal bulge which modifies the spin and the orbit of the planet.
The ultimate stage for tidal evolution is the synchronization of the rotation and orbital periods, alignment of the planet spin axis with the normal to the orbit (zero planet obliquity), and the circularization of the orbit (e.g., Hut 1980; Adams & Bloch 2015). Although the spin of closein planets quickly evolves into an equilibrium configuration, the orbital evolution is much slower (e.g., Correia & Laskar 2010). The circularization characteristic timescale is given by Eq. (20) using a constantQ tidal model (1)
where P_{orb} is the orbital period, m is the mass of the planet, M is the mass of the star, a is the semimajor axis of the orbit, R is the average radius of the planet, and k_{2} is the second Love number for potential. Q is the tidal qualityfactor, which measures the ratio of energy dissipated during one period of tidal stress over the peak energy stored in the system during the same period.
For Uranus and Neptune, Q∕k_{2} ~ 10^{5} (Tittemore & Wisdom 1990; Banfield & Murray 1992). Assuming that Neptunemass planets have similar rheologies, we obtain (Eq. (1)) (2)
According to this expression, all Neptunemass planets with P_{orb} < 5 day should circularize in less than 5 Gyr. However, observations of Neptunemass planets such as GJ 436 b put this scenario into question. In a shortperiod orbit with P_{orb} = 2.644 day (Lanotte et al. 2014) around a cool star with M = 0.445 M_{⊙} (Mann et al. 2015), GJ 436 b has a circularization timescale of τ_{circ} ~ 0.2 Gyr (Eq. (2)). This is much shorter than the stellar age of about 4–8 Gyr (Bourrier et al. 2018a), yet the planetary orbit has an eccentricity e = 0.162 ± 0.004 (Lanotte et al. 2014), which is clearly different from zero. Several similar cases have been reported in the literature, such as HD 125612 c (Lo Curto et al. 2010), HATP26 b (Hartman et al. 2011), HATP11 b (Yee et al. 2018), and GJ 3470 b (Kosiarek et al. 2019). In fact, among all Neptunesize planets with P_{orb}< 5 day (Fig. 1), only HD 219828 b (Melo et al. 2007) and HD 47186 b (Bouchy et al. 2009) exhibit a small measured eccentricity (e ≲ 0.05), which is still compatible with a nonzero value.
Figure 1 shows measurements of eccentricity for planets with P_{orb}< 100 day, divided into three groups: Earthsize (R_{⊕} < R < 3 R_{⊕}); Neptunesize (3 R_{⊕}< R < 9 R_{⊕}); and Jupitersize (9 R_{⊕}< R < 25 R_{⊕}). The distribution of Neptunesize planets with P_{orb} < 5 day contrasts with that given for closein Jupiter or Earthsize planets. First, there are no warm Neptune planets with P_{orb}< 3 day, which corresponds to the wellknown “Neptunian desert” (Lecavelier des Etangs 2007; Davis & Wheatley 2009; Szabó & Kiss 2011; Mazeh et al. 2016). We note that the lack of eccentricity values for Earthsize planets is due to the difficulty in measuring them at a high precision level. Second, as noted above, all warm Neptune planets with P_{orb}< 5 day have eccentricities that are consistent with nonzero values despite having damping timescales shorter than 5 Gyr. Finally, there is no apparentcorrelation between orbital period and eccentricity; on average, planets with shorter orbital periods are expected to have smaller eccentricities because they have shorter damping timescales (see Eq. (2)).
We observe that a fraction of Jupiter and Earthsize planets with P_{orb}< 5 day are also observed in eccentric orbits. Nevertheless, the eccentricity values increase with the orbital period, as expected, due to a more efficient tidal damping at smaller distances (Eq. (1)). In the case of Jupitersize planets, the eccentricity leftovers are likely related to their formation process through Lidov–Kozai cycles (e.g., Fabrycky & Tremaine 2007) or planet–planet scattering (e.g., Beaugé & Nesvorný 2012). Indeed, for P_{orb}> 5 day, we observe that the eccentricity of Jupitersize planets can reach extremely high values, while for the lower mass population it remains smaller than about 0.2.
One possibility to explain the observed eccentricities for warm Neptunes is that we are overestimating the tidal dissipation for these planets. If we adopt Q∕k_{2} > 10^{7}, the present nonzero eccentricities can be simply explained as residual values excited during the formation process. However, such high values for Q are not observed for similar mass planets in our Solar System (e.g., Tittemore & Wisdom 1990; Banfield & Murray 1992). In addition, this explanation also fails to explain the “Neptunian desert” and the absence of correlation between orbital period and eccentricity for P_{orb} < 5 day (Sect. 2). Another possibility to explain the anomalous eccentricity distribution is that additional mechanisms pump the eccentricity or delay the circularization of the orbit. In this paper, we discuss the most likely mechanisms that act in opposition to gravitational tides, namely, thermal atmospheric tides (Sect. 3), evaporation of the atmosphere (Sect. 4), and excitation from a distant companion (Sect. 5). We discuss their limits of applicability in Sect. 6.
Fig. 1 Distribution of eccentricities as a function of orbital period. Only eccentricities measured with uncertainties smaller than 0.1 are shown. Planets are arbitrarily separated in three categories by radius: Jupitersize (red, top panel), Neptunesize (green, middle panel), and Earthsize (blue, bottom panel). In the middle panel, black squares highlight three iconic planets representative of the warm Neptunes population. 
2 High initial eccentricity
The damping timescale strongly depends on the orbital period, (Eq. (2)). Then, if the observed eccentric planet was initially on a wider and more eccentric orbit, the initial damping timescale is much longer than the τ_{circ} obtained with the current orbital period. In this section, we estimate this effect.
We assume that tides in the planet circularize its orbit following the law (see Eq. (20)) (3)
Adopting the constantQ model^{1} (Eq. (1)), we have α = 13∕2. We also assume that the tidal effect in the planet does not affect its orbital angular momentum significantly. We thus have (5)
where a_{0}, e_{0}, and τ_{0} are the current (at the time of observation) values of the parameters. In the case of small eccentricity, this equation can be approximated by (7)
which is equivalent to assuming a constant circularization timescale τ_{circ}. The evolution of the eccentricity in this regime is the classical exponential decay, (8)
On the contrary, if the eccentricity is close to one, Eq. (6) can be approximated by (9)
with x = 1 − e^{2}, and we find (10)
Therefore, asymptotically (for e → 1), we have t →−∞. This means that we can make the tidal circularization process arbitrarily long by taking an initial eccentricity very close to one, together with a very large initial semimajor axis.
In the general case, the evolution of the eccentricity follows (11)
which can be rewritten as (12)
is the hypergeometric function. The above expression allows us to estimate the time needed to reach the currently observed orbit (a_{0}, e_{0}) as a function of the initial eccentricity e. Figure 2 shows an example of the evolution of the semimajor axis and eccentricity as a function of time. We assume the observed eccentricity to be e_{0} = 0.1, typical of the values measured for warm Neptunes (Fig. 1). As expected, for small eccentricities (e ≲ 0.4), we observe an exponential decay of the eccentricity (Eq. (8)), while at high eccentricities (e ≳ 0.6) the evolution is much slower (Eq. (10)). This counterintuitive behaviour is due to the fact that when the eccentricity was large, the semimajor axis was also large and, thus, tidal effects were much less efficient.
In Fig. 2, we see that if the initial eccentricity was around 0.8, and the initial semimajor axis around 2.8 times the observed one, the planet might have spent around 40 τ_{0} to reach the observed configuration (where τ_{0} is the current instantaneous damping timescale). For a warm Neptune like GJ 436 b, this corresponds to about 8 Gyr. Therefore, any mechanism (scattering, Kozai, tidal pumping, etc.) that could have excited the eccentricity to a high level during the formation of the planet (or at a lower level more recently) might be responsible for the currently remaining eccentricity.
A possible limitation of this scenario is that we have to observe the system at a specific time in its evolution. Indeed, the eccentricity spends a long time above 0.6, and a long time close to 0, but a short time inbetween (see Fig. 2). The probability to observe the system at this specific time (i.e. during the phase of exponential decay) is about τ_{0} ∕T, where T is the star age. Therefore, if the observed population of warm Neptunes with eccentricity of ~0.1 formed this way (Fig. 1), one would expect to also observe sizeable populations of Neptunemass planets in the high eccentricity phase (on wider orbits) and in the circular phase (on closerin orbits).
Fig. 2 Evolution (past and future) of the eccentricity (top) and semimajor axis (bottom) of a planet undergoing tidal dissipation with a power law ė ∝−ea^{−13∕2} (i.e. α = 13∕2). The timescale τ_{0} is the currently observed instantaneous damping timescale (in this example at a = a_{0}, e = e_{0} =0.1). However,the shape of the curve is always the same and does not depend on e_{0}. 
3 Thermal atmospheric tides
The differential absorption of the stellar heat by the atmosphere of a planet gives rise to local variations of temperature and, consequently, to pressure gradients. The mass of the atmosphere is then permanently redistributed, adjusting for an equilibrium position. More precisely, the particles of the atmosphere move from the high temperature zone (at the substellar point) to the low temperature areas. Observations on Earth show that the pressure redistribution is essentially asuperposition of two pressure waves (Chapman & Lindzen 1970; AuclairDesrotour et al. 2017): a daily (or diurnal) tide of small amplitude (the pressure is minimal at the subsolar point and maximal at the antipode) and a strong halfdaily (semidiurnal) tide (the pressure is minimal at the subsolar point and at the antipode).
As for bodily tides, there is a delay between the response of the atmosphere and the thermal excitation. The phase shifted elongation of the shell induces a torque which modifies the rotational dynamics of the body as well as its orbit. For atmospheric tides, the dominating semidiurnal pressure wave usually leads the perturbation. As a consequence, the thermal tide moves the planet rotation away from the synchronous equilibrium, and determines possible new states of equilibrium balanced by bodily tides (Gold & Soter 1969; Correia & Laskar 2001).
The combined effect of bodily and thermal tides tend to align the equator of the planet with the orbital plane (Correia et al. 2003). Therefore, for simplicity, we let the obliquity of the planet to be zero. In this case, the secular evolution of the rotation rate, Ω, and semimajor axis, a, due to a given kind of tidal effect (bodily or thermal) can be written for small eccentricity, e, as (Correia & Laskar 2003, 2010): (14) (15)
where terms in e^{2} have been neglected. C is the principal moment of inertia, L = βna^{2} is the orbital angular momentum, n is the mean motion, and β = mM∕(m + M) ≈ m is the reduced mass of the system. K_{τ} is a constant factor related to the strength of the tide, and b_{τ}(x) is an odd function related to the dissipation within the planet. For bodily (or gravitational) tides we denote τ ≡ g, while for thermal (atmospheric) tides we denote τ ≡ a. For the first kind we have (16)
while for the second kind (17)
where G is the gravitational constant. k_{2} (x) > 0 is the second Love number, p_{2} (x) > 0 is the amplitude of the pressure variations at the ground, and 0 ≤ δ_{τ}(x) < 90° is the phase shift. The minus sign in the expression of b_{a}(x) already accounts for the fact that thermal tides lead the perturbation (Correia et al. 2003; AuclairDesrotour et al. 2017). The phase shift can be related with the time delay Δt_{τ} between the perturbation and the maximal amplitude of the tide through δ_{τ} (x) = xΔt_{τ}(x). We thus have b_{τ}(0) = 0, and for x > 0, b_{g} (x) > 0 and b_{a} (x) < 0. The phase shift can also be related to the quality factor through (Efroimsky 2012, Eq. (141)).
For the eccentricity evolution, bodily and thermal tides give different evolutions^{2}. For the first we have (18)
If we consider bodily tides alone, from Eq. (14) we see that the equilibrium configuration is obtained when for Ω = n, that is, for synchronous rotation. In addition, this equilibrium is stable because b_{g} (x) is an odd function for which b_{g}(x) > 0 when x > 0. Replacing Ω = n in Eqs. (15) and (18), we get ȧ = 0, and (20)
which leads to an exponential damping of the eccentricity, as expected. For thermal tides alone, we can derive similar conclusions, but Ω = n becomes an unstable equilibrium and (Eq. (19)) (21)
since b_{a}(x) < 0 when x > 0. Hence, we conclude that thermal tides may be able to counterbalance the damping effect of bodily tides on the eccentricity.
When bodily and thermal tides are considered together, the complete evolution of the rotation rate is given by (Eq. (14)) (22)
New equilibrium solutions for the rotation rate () are then obtained for (Correia & Laskar 2001) (23)
where f(x) = b_{a}(x)∕b_{g}(x). This simple condition explains why the spin of the planet Venus is not synchronous. Let us express the stable solutions of Eq. (23) as Ω = n + Ω_{0}∕2. For the stable states we then have f(Ω_{0}) = −K_{g}∕K_{a}. Thus, for the semimajor axis secular evolution we get (Eq. (15)) (24)
and for the eccentricity (Eqs. (18) and (19)) (25)
The eccentricity evolution depends on the balance between the bodily and thermal tides harmonic functions b_{τ} (x), that is, it depends on the tidal models that we adopt. As the rheology of planets is poorly known, these functions are usually unknown. However, for closein planets it is expected that bodily tides dominate thermal tides (Correia et al. 2008; Cunha et al. 2015; Leconte et al. 2015) and, thus, Ω_{0}≪ n. We can then simplify the previous expression as (26)
that is, the eccentricity increases whenever (27)
We immediately see that the eccentricity always decreases if we adopt for both tides a constant (b_{τ} = cte) or a linear model (b_{τ} ∝ x). On the other hand, if we adopt a linear model for thermal tides and a constant model for bodily tides, the eccentricity always increases, since we get Ω_{0}∕n < 6∕7, which is consistent with Ω_{0}≪ n. The constant and linear models are widely used, because they are simple, but they are not very realistic. A more correct description of the rheology involves the use of viscoelastic models (e.g., Henning et al. 2009). In particular, thermal tides are believed to follow a Maxwell rheology (Leconte et al. 2015; AuclairDesrotour et al. 2017), for which , where τ is a viscous relaxation time. Assuming also a Maxwell model for bodily tides, with (Remus et al. 2012; Correia et al. 2014), the eccentricity increases for (Eq. (27)) (28)
The value of τ_{g} can be estimated from Q ~ τ_{g}n + 1∕(τ_{g}n), which has two solutions. With Q ~ 10^{4}, for orbital periods P_{orb} ~ 5 day, we get τ_{g} ~ 10^{4} day or τ_{g} ~ 10^{−4} day. At present, there are no estimations of τ_{a} for Neptunelike planets. However, for an Earthlike planet, we get τ_{a} ~ 10 day (Leconte et al. 2015). Assuming an identical value for warm Neptunes, we see that condition in Eq. (29) is widely assured for τ_{g} ~ 10^{4} day. We then conclude that although there is a limited set of values (τ_{g}, τ_{a}) that allow an increase in the eccentricity due to the combined action of bodily and thermal tides, the rheology of warm Neptunes may be such that condition in Eq. (29) can be met.
4 Evaporation of the atmosphere
Interestingly, it has been observed that the atmospheres of some of the warm Neptunes with nonzero eccentricity are also undergoing strong evaporation. This is the case for the planets GJ 436 b (Ehrenreich et al. 2015; Lavie et al. 2017; dos Santos et al. 2019) and GJ 3470 b (Bourrier et al. 2018b). Therefore, a fraction of the mass of these planets is being lost through this process. Isotropic planetary evaporation generates almost no orbital evolution. However, the hottest region in a planet atmosphere could shift with respect to the substellar point, usually eastward, depending on the strength of the imposed stellar heating and other factors such as the radiative timescale and drag time constant (e.g., Showman & Guillot 2002; Showman & Polvani 2011). This displacement could be associated with anisotropic massloss, and a modification in the orbit of the planet (Boué et al. 2012; Teyssandier et al. 2015).
We let θ be the angle between the star and the direction of maximal mass ejection, φ be the solid angle aperture of the stream, and v_{e} be the ejection speed of atmospheric particles with respect to the planet’s barycenter (see Fig. 1 in Boué et al. 2012). The secular evolution of the semimajor axis and eccentricity due to the massloss can be written as (Boué et al. 2012): (30) (31)
where terms in e^{2} have been neglected, and (32)
is a dimensionless parameter related to the efficiency of the evaporation. We have γ = 0 for substellar (θ = 0, π) or isotropic evaporation (φ = π). We also note that the orbital speed v_{orb}≈ na. Since we assume that the planet is losing mass (ṁ < 0), this leads to an increase in both semimajor axis and eccentricity.
The secular evolution of the eccentricity due to the combined effect of bodily tides (Eq. (20)) and evaporation (Eq. (31)) is then (33)
For GJ 436 b and GJ 3470 b we need − γ ṁ∕m > 10 Gyr^{−1} and 3.3 Gyr^{−1}, respectively (Eq. (2)). The two planets have similar orbital velocities, v_{orb} ≈ 120 km s^{−1}. Theoretical simulations of escaping outflows yield typical escape velocities v_{e} between 1− 10 km s^{−1} (e.g., Salz et al. 2016). Simulations of GJ 436 b exosphere based on Lymanα transit observations suggest that v_{e} could reach up to 50− 60 km s^{−1}. Adopting the upper limit, γ still remains lower than about 1. The present mass loss rates from GJ 436 b and GJ 3470 b atmospheres are 2.2 × 10^{10} g s^{−1} (Bourrier et al. 2016) and 8.5 × 10^{10} g s^{−1} (Bourrier et al. 2018b), respectively. These values are also upper limits, assuming the planets’ atmospheres are in the energylimited regime and all input stellar energy is used for the escape. Even with maximum values, − γ ṁ∕m would thus reach 0.0025 Gyr^{−1} for GJ 436 b and 0.016 Gyr^{−1} for GJ 3470 b, far too low for anisotropic escape to balance the damping effect from tides. Indeed, in the optimistic case for which γ ~ 1, the planet would have to lose roughly half of its mass every 100 Myr to balance the damping effect from tides; after a few Myr the planet would have completely lost its atmosphere. We then conclude that, unless there is an unknown mechanism that produces an ejection speed of v_{e} > 10^{3} km s^{−1}, the orbital effect of atmospheric evaporation is not sufficient to explain the eccentricities observed for warm Neptunes.
5 Excitation from a distant companion
Most warm Neptune eccentric planets were discovered using radial velocity combined with transit observations. At present, these are the two most successful techniques for detecting planets, but they share a caveat: they are unable to spot planets whose orbital plane lie near the observer’s plane of the sky. In the case of radial velocity, the amplitude of the signal is proportional to the sinus of the inclination of the orbital plane to the plane of the sky, while for transits we can only detect planets passing in front of the star. Therefore, we cannot rule out that a companion is present in the system, provided that it lies in a inclined or very distant orbit with respect to the warm Neptune. Indeed, some warm Neptune planets have been observed in systems with more massive distant companions (HATP11, HD 219828, HD 47186 and HD 125612). We note, however, that even in these positive cases, we also lack knowledge for their true masses and inclinations, since they were all spotted by radialvelocity alone. In other cases, such as GJ 436, GJ 3470 and HATP26, no companion has been identified so far.
In planetary dynamics, still undetected companions are often suggested to solve unexplained observations. For instance, a ninth planet was proposed to explain the observed clustering in the perihelia of the distant Kuiper Belt objects (Trujillo & Sheppard 2014; Batygin & Brown 2016). A great success was the discovery of the planet Neptune by Le Verrier (1846), that probably triggered all the other undetected companion scenarios. A great failure was the prediction of planet Vulcan to explain the perihelion advance of Mercury, also by Le Verrier (1859), which was later explained in the frame of general relativity (Einstein 1915). Therefore, the undetected companion scenario should always be considered with a degree of caution.
When the orbit of a planet is excited by an outer companion in a inclined orbit, its eccentricity undergoes some perturbations. Assuming for simplicity that the perturber is on a circular orbit with radius a_{p}, mutual inclination I, and its mass m_{p} ≫ m, the dominating (quadrupolar) contribution for the eccentricity and pericenter are (e.g., Correia et al. 2013): (35) (36)
is due to the perturber, while (38)
are due to general relativity and planet oblateness, respectively. The mutual inclination can also change over time, but it can be related with the eccentricity using the conservation of the total orbital angular momentum: (39)
also known by the “Kozai constant”.
5.1 Lidov–Kozai cycles
When a massive companion in a inclined orbit is present, there is an alternative mechanism to the migration in the accretion disk for the formation of closein planets: high eccentricity migration through Lidov–Kozai cycles combined with tidal friction (e.g., Wu & Murray 2003; Correia et al. 2011; Anderson et al. 2016). This mechanism has been proposed byBeust et al. (2012) to explain the present eccentricity of GJ 436 b and it is strengthened by the measurement of the planet misaligned orbit (Bourrier et al. 2018a), which is also a natural consequence of the Kozai cycles.
For planets far from the star, the gravitational perturbations from the companion dominate over the general relativity and the oblateness terms (ν ≫ ν_{g}, ν_{r}) in the precession of the periapsis (Eq. (36)). As a consequence, it is possible to find stable equilibria points for the eccentricity (ė = 0, ), known as Kozai equilibria (Lidov 1962; Kozai 1962). These points are centered at ω = ±π∕2 and . They correspond to a resonance between the precession of the longitude of the node and the longitude of the pericenter of the inner orbit. Therefore, there is a libration zone associated to this resonance (for ), a separatrix, and a nonresonant circulation zone.
Higheccentricity migration through Lidov–Kozai cycles is very efficient, provided that the planet is close to the separatrix of the Lidov–Kozai equilibria. Indeed, near the separatrix, the mutual inclination undergoes large variations, that, in turn, induce large eccentricity oscillations (Eq. (39)). Even at large distances, the outer companion can significantly perturb the inner orbit as long as the initial mutual inclination is . If the inner orbit is initially circular, the maximum eccentricity achieved in a Lidov–Kozai cycle is (e.g., Lidov & Ziglin 1976). When the eccentricity reaches very high values tidal dissipation is enhanced and the semimajor axis of the inner orbit decreases.
Lidov–Kozai cycles persist as long as the perturbation from the outer companion is the dominant cause of periapse precession in the planetary orbit. However, small additional sources of periapse precession, such as general relativity (ν_{g}) or oblateness (ν_{r}), can compensate the gravitational precession (ν) and suppress the eccentricity oscillations (Eq. (35)). The inner planet is then left in a closein orbit with large eccentricity that is then progressively damped by tides (see Sect. 2).
The presence of the Lidov–Kozai cycle from an undetected companion then act as possible mechanism of formation of the warm Neptune that simultaneously accounts for the possible observed eccentricity, as a remnant of this formation mechanism. Although the circularization timescale for GJ 436 b is only about 0.2 Gyr (Eq. (2)), there is a range of parameter values for its companion for which the Lidov–Kozai cycles could delay the deliver of the planet close to the star; that would explain why we still observe e ≈ 0.16 at present (Bourrier et al. 2018a).
5.2 Spindriven eccentricity pumping
This effect has been first described for coplanar orbits (Correia et al. 2012; Greenberg et al. 2013), but it is also much more efficient for mutually inclined orbits (Correia et al. 2013).
When the general relativity or the oblateness terms dominate over the gravitational perturbations (ν_{g} or ν_{r} ≫ ν), the Lidov–Kozai cycles no longer work. This is the case for all currently observed warm Neptunes. Thus, at first order, the precession of the periastron is constant, (Eq. (36)), and the eccentricity shows small oscillations around a mean value, e = e_{0} + δe, with (40)
where e_{0}, Δe and ω_{0} are constant. The evolution of the rotation rate (Eq. (14)) also depends on the eccentricity^{3}, so the rotation rate experiences a variation Ω = Ω_{0} + δΩ, with (41)
where Ω_{0}, ΔΩ and ϕ are constant. The rotation rate thus presents an oscillation identical to the eccentricity (Eq. (40)), but delayed by an angle ϕ, such that (Correia et al. 2013) (42)
The damping timescale for the spin can be seen as τ_{spin} ≈ k^{−1}, which is usually much shorter than the circularization timescale, τ_{circ} ≈ LQ∕k_{2}K_{g} (Eq. (1)), since L ≫ Cn. At second order, the precession of the periastron depends on the eccentricity (Eq. (36)), but also on the rotation rate (term in ν_{r}). Replacing Eqs. (40) and (41) in the expression of the periastron (Eq. (36)), integrating and replacing the result again in Eq. (40), we obtain a secular contribution for the eccentricity (see Appendix A in Correia et al. 2013) (43)
Thus, the secular term vanishes when ϕ = 0 or π∕2, that is, for strong dissipation (k ≫ g), where δė ~ ν_{r}g∕k, or for weak dissipation (k ≪ g), where δė ~ ν_{r}k∕g, respectively. The effect on the eccentricity is then maximized when ϕ = π∕4, which occurs for k ~ g, for which δė ~ ν_{r}. We conclude that, in order to observe the spindriven eccentricity pumping effect, the damping timescale of the spin (k^{−1}) should be of the order of the period of eccentricity oscillations (g^{−1}).
The secular evolution that results from expression (43) must be added to the usual orbital damping (Eq. (20)). As long as the spindriven term counterbalances the damping effect, the eccentricity can increase to high values. However, when we consider the full nonlinearized problem, the secular increase (Eq. (43)) cannot last indefinitely because when the eccentricity reaches high values, the damping effect is also enhanced. As a consequence, the spindriven pumping of the eccentricity is never permanent, although it can last through the age of the system.
6 Discussion
The distribution of eccentricity measured for warm Neptunes with P_{orb} < 5 day contrasts with that given for smaller and larger closein planets. In particular, most warm Neptunes present nonzero eccentricity, a surprising feature considering that bodily tides should have circularized their orbits in a timescale shorter than the age of the systems. In this paper, we investigate several mechanisms that could counterbalance the bodily tides, namely, thermal atmospheric tides, evaporation of the atmosphere, and excitation from a distant companion.
The combined action of bodily and thermal tides can be simplified using their relaxation times, (τ_{g}, τ_{a}), respectively. A dynamical analysis shows that the eccentricity is allowed to increase under some conditions for these two parameters, in particular, for τ_{a} < τ_{g} (Eq. (29)). Observational constraints on the rheology and atmospheric composition of warm Neptunes are required to assess the fraction of these planets maintained on eccentric orbits because of thermal tides.
Spectacular evaporation has been observed for the warm Neptunes GJ 436 b (Bourrier et al. 2018a) and GJ 3470 b (Bourrier et al. 2018b). This class of planets is expected to be particularly sensitive to atmospheric escape (Owen & Wu 2017). If anisotropic, their evaporation would always contribute to increase the orbital eccentricity. Our estimates, however, show that the velocity and mass loss rate of the gas escaping from warm Neptunes (as derived from observations and predicted theoretically) are not sufficient to counterbalance the damping effect from bodily tides.
Distant, more massive planets have been found to accompany some warm Neptunes. Those that appear alone may still have undetected companions. In that case, mutual gravitational interactions between the two planets can temporarily excite the eccentricity of the inner warm Neptune, thus delaying the circularization of the orbit up to several Gyrs. This excitation can be induced via Lidov–Kozai cycles for companions with mutual inclinations , or via spindriven eccentricity pumping when the damping timescale of the spin (k^{−1}) is of the order of the period of eccentricity oscillations (g^{−1}). Because they excite the eccentricity for a limited amount of time, these mechanisms alone cannot explain the distribution of eccentricity observed for warm Neptunes, in particular the absence of closein Neptunes on circular orbits.
Although all the proposed scenarios have some limitations, we cannot rule out that more than one of these mechanisms is simultaneously at work. It is possible that the hot expanding atmosphere of an evaporating Neptune creates an additional thermal tidal torque on the planet that would help to increase its eccentricity. Alternatively, a distant companion can excite the eccentricity, and as the Neptunesize planet migrates closer to the star, it would evaporate increasingly and erode into a smaller planet before it can fully circularize, thus explaining why no Neptunemass planets are found on circular orbits close to their star. This combination of higheccentricity migration and evaporation as the origin of eccentric warm Neptunes was proposed by Bourriere t al. (2018a), based on the study of GJ 436b, and is now gaining interest to explain the structure of the “Neptune desert” (e.g., Owen & Lai 2018). In this scenario, the warm Neptunes on eccentric orbits are either still undergoing their migration and will erode as they enter the desert or they have already reached a stable orbit far enough from the star to be safe from evaporation. The orbital properties and evaporation status of the warm Neptunes GJ 436 b (Bourrier et al. 2018a) and GJ 3470 b (Bourrier et al. 2018b) are consistent with this scenario.
In conclusion, thermal atmospheric tides or excitation from a distant companion combined with evaporation could explain the distribution of eccentricities observed for warm Neptunes. These scenarios require specific conditions for the warm Neptunes, their host star, and possible companions. At present we cannot conclude whether one of these mechanisms is dominant, whether they act together, or whether each influence a fraction of the observed warm Neptunes.
Further theoretical modeling of these mechanisms should be attempted to better understand their impact, complemented by indepth characterizations of the orbital and atmospheric properties of known warm Neptunes, along with a search for additional Neptunesize planets and their putative companions. Studying the distribution of eccentricities as a function of the tidal quality factor could also allow us to better assess the role of tidal circularization for each class of planets, though this requires a larger sample of highprecision mass and radius measurements, and a finer understanding of planetary internal structures.
Acknowledgements
A.C. acknowledges support by CFisUC strategic project (UID/FIS/04564/2019), ENGAGE SKA (POCI010145FEDER022217), and PHOBOS (POCI010145FEDER029932), funded by COMPETE 2020 and FCT, Portugal. V.B. acknowledges support by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Four Aces grant agreement No 724427). J.B.D. acknowledges support by the Swiss National Science Foundation (SNSF). This work has, in part, been carried out within the framework of the National Centre for Competencein Research PlanetS supported by SNSF.
References
 Adams, F. C., & Bloch, A. M. 2015, MNRAS, 446, 3676 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, K. R., Storch, N. I., & Lai, D. 2016, MNRAS, 456, 3671 [NASA ADS] [CrossRef] [Google Scholar]
 AuclairDesrotour, P., Laskar, J., & Mathis, S. 2017, A&A, 603, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Banfield, D., & Murray, N. 1992, Icarus, 99, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Brown, M. E. 2016, AJ, 151, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Beaugé, C., & Nesvorný, D. 2012, ApJ, 751, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Beust, H., Bonfils, X., Montagnier, G., Delfosse, X., & Forveille, T. 2012, A&A, 545, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bouchy, F., Mayor, M., Lovis, C., et al. 2009, A&A, 496, 527 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boué, G., Figueira, P., Correia, A. C. M., & Santos, N. C. 2012, A&A, 537, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., Tanaka, Y. A., & Vidotto, A. A. 2016, A&A, 591, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., Lovis, C., Beust, H., et al. 2018a, Nature, 553, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., et al. 2018b, A&A, 620, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chapman, S., & Lindzen, R. 1970, Atmospheric Tides, Thermal and gravitational (Dordrecht: Springer) [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2003, J. Geophys. Res. (Planets), 108, 5123 [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2010, Tidal Evolution of Exoplanets, ed. S. Seager (Tucson: University of Arizona Press), 239 [Google Scholar]
 Correia, A. C. M., Laskar, J., & Néron de Surgy, O. 2003, Icarus, 163, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Udry, S., Mayor, M., et al. 2008, A&A, 479, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dyn. Astron., 111, 105 [Google Scholar]
 Correia, A. C. M., Boué, G., & Laskar, J. 2012, ApJ, 744, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Boué, G., Laskar, J., & Morais, M. H. M. 2013, A&A, 553, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cunha, D., Correia, A. C. M., & Laskar, J. 2015, Int. J. Astrobiol., 14, 233 [CrossRef] [Google Scholar]
 Davis, T. A., & Wheatley, P. J. 2009, MNRAS, 396, 1012 [NASA ADS] [CrossRef] [Google Scholar]
 dos Santos, L. A., Ehrenreich, D., Bourrier, V., et al. 2019, A&A, 629, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Efroimsky, M. 2012, Celest. Mech. Dyn. Astron., 112, 283 [NASA ADS] [CrossRef] [Google Scholar]
 Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459 [Google Scholar]
 Einstein, A. 1915, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), Seite 844 [Google Scholar]
 Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298 [NASA ADS] [CrossRef] [Google Scholar]
 Gold, T., & Soter, S. 1969, Icarus, 11, 356 [NASA ADS] [CrossRef] [Google Scholar]
 Greenberg, R., Van Laerhoven, C., & Barnes, R. 2013, Celest. Mech. Dyn. Astron., 117, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Hartman, J. D., Bakos, G. Á., Kipping, D. M., et al. 2011, ApJ, 728, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Henning, W. G., O’Connell, R. J., & Sasselov, D. D. 2009, ApJ, 707, 1000 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1980, A&A, 92, 167 [NASA ADS] [Google Scholar]
 Kosiarek, M. R., Crossfield, I. J. M., HardegreeUllman, K. K., et al. 2019, AJ, 157, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
 Lanotte, A. A., Gillon, M., Demory, B.O., et al. 2014, A&A, 572, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lavie, B., Ehrenreich, D., Bourrier, V., et al. 2017, A&A, 605, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Le Verrier, U. J. 1846, Astron. Nachr., 25, 65 [NASA ADS] [Google Scholar]
 Le Verrier, U. J. 1859, Annales de l’Observatoire de Paris, 5 [Google Scholar]
 Lecavelier des Etangs, A. 2007, A&A, 461, 1185 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Lidov, M. L. 1962, Planet. Space Sci., 9, 719 [NASA ADS] [CrossRef] [Google Scholar]
 Lidov, M. L., & Ziglin, S. L. 1976, Celest. Mech., 13, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Lo Curto, G., Mayor, M., Benz, W., et al. 2010, A&A, 512, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MacDonald, G. J. F. 1964, Rev. Geophys., 2, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015, ApJ, 804, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Melo, C., Santos, N. C., Gieren, W., et al. 2007, A&A, 467, 721 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owen, J. E., & Lai, D. 2018, MNRAS, 479, 5012 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Remus, F., Mathis, S., Zahn, J.P., & Lainey, V. 2012, A&A, 541, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016, A&A, 586, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Showman, A. P., & Guillot, T. 2002, A&A, 385, 166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Szabó, G. M., & Kiss, L. L. 2011, ApJ, 727, L44 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssandier, J., Owen, J. E., Adams, F. C., & Quillen, A. C. 2015, MNRAS, 452, 1743 [NASA ADS] [CrossRef] [Google Scholar]
 Tittemore, W.C., & Wisdom, J. 1990, Icarus, 85, 394 [NASA ADS] [CrossRef] [Google Scholar]
 Trujillo, C. A., & Sheppard, S. S. 2014, Nature, 507, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y., & Murray, N. 2003, ApJ, 589, 605 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yee, S. W., Petigura, E. A., Fulton, B. J., et al. 2018, AJ, 155, 255 [NASA ADS] [CrossRef] [Google Scholar]
For closein planets, the orbital mean motion does not change much during the circularization process (see Fig. 2), so the constantQ model is a good approximation. Alternative tidal models provide different α, which does not modify the final conclusions of this section.
The bodily tides potential is proportional to r^{−6}, while the thermal tides potential is proportional to r^{−5}, where r is the relative distance between the planet and the star (see Correia et al. 2003).
In Eq. (14) we neglected terms in e^{2}, but expanding to a higher order we get (Correia et al. 2008): .
All Figures
Fig. 1 Distribution of eccentricities as a function of orbital period. Only eccentricities measured with uncertainties smaller than 0.1 are shown. Planets are arbitrarily separated in three categories by radius: Jupitersize (red, top panel), Neptunesize (green, middle panel), and Earthsize (blue, bottom panel). In the middle panel, black squares highlight three iconic planets representative of the warm Neptunes population. 

In the text 
Fig. 2 Evolution (past and future) of the eccentricity (top) and semimajor axis (bottom) of a planet undergoing tidal dissipation with a power law ė ∝−ea^{−13∕2} (i.e. α = 13∕2). The timescale τ_{0} is the currently observed instantaneous damping timescale (in this example at a = a_{0}, e = e_{0} =0.1). However,the shape of the curve is always the same and does not depend on e_{0}. 

In the text 
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.