MercuryT: A new code to study tidally evolving multiplanet systems. Applications to Kepler62^{⋆}
^{1}
Univ. Bordeaux, LAB, UMR 5804,
33270
Floirac,
France
^{2}
CNRS, LAB, UMR 5804, 33270
Floirac,
France
^{3}
Now at: NaXys, Department of Mathematics, University of
Namur, 8 Rempart de la
Vierge, 5000
Namur,
Belgium
email:
emeline.bolmont@unamur.be
^{4} Canadian Institute for Theoretical Astrophysics, 60st St
George Street, University of Toronto, Toronto, ON, M5S3H8, Canada
^{5}
Center for Planetary Sciences, Department of Physical &
Environmental Sciences, University of Toronto Scarborough,
Toronto, ON, M1C
1A4, Canada
^{6}
CIDMA, Departamento de Física, Universidade de
Aveiro, Campus de
Santiago, 3810193
Aveiro,
Portugal
^{7}
ASD, IMCCECNRS UMR 8028, Observatoire de Paris,
UPMC, 77 Av. DenfertRochereau,
75014
Paris,
France
Received: 16 February 2015
Accepted: 14 July 2015
A large proportion of observed planetary systems contain several planets in a compact orbital configuration, and often harbor at least one closein object. These systems are then most likely tidally evolving. We investigate how the effects of planetplanet interactions influence the tidal evolution of planets. We introduce for that purpose a new opensource addition to the MercuryNbody code, MercuryT, which takes into account tides, general relativity and the effect of rotationinduced flattening in order to simulate the dynamical and tidal evolution of multiplanet systems. It uses a standard equilibrium tidal model, the constant time lag model. Besides, the evolution of the radius of several host bodies has been implemented (brown dwarfs, Mdwarfs of mass 0.1 M_{⊙}, Sunlike stars, Jupiter). We validate the new code by comparing its output for oneplanet systems to the secular equations results. We find that this code does respect the conservation of total angular momentum. We applied this new tool to the planetary system Kepler62. We find that tides influence the stability of the system in some cases. We also show that while the four inner planets of the systems are likely to have slow rotation rates and small obliquities, the fifth planet could have a fast rotation rate and a high obliquity. This means that the two habitable zone planets of this system, Kepler62e ad f are likely to have very different climate features, and this of course would influence their potential at hosting surface liquid water.
Key words: planets and satellites: dynamical evolution and stability / planets and satellites: terrestrial planets / planets and satellites: individual: Kepler 62 / planetstar interactions
The code is only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/583/A116
© ESO, 2015
1. Introduction
More than 1400 exoplanets have now been detected and about 20% of them are part of multiplanet systems^{1}. Many of these systems are compact and host closein planets for which tides have an influence. In particular, tides can have an effect on the eccentricities of planets, and also on their rotation periods and their obliquities, which are important parameters for any climate studies. Besides tides can influence the stability of multiplanet systems, by their effect on the eccentricities but also by their effect on precession rates.
We present here a new code, MercuryT^{2}, which is based on the Nbody code Mercury (Chambers 1999). It allows us to calculate the evolution of semimajor axis, eccentricity, inclination, rotation period and obliquity of the planets as well as the rotation period evolution of the host body. This code is flexible in so far as it allows to compute the tidal evolution of systems orbiting any nonevolving object (provided we know its mass, radius, dissipation factor and rotation period), but also evolving brown dwarfs (BDs), an evolving Mdwarf of 0.1 M_{⊙}, an evolving Sunlike star, and an evolving Jupiter.
The dynamics of multiplanet systems with tidal dissipation have been studied (evolution of the orbit in Wu & Goldreich 2002; Mardling 2007; Batygin et al. 2009; Mardling 2010; Laskar et al. 2012 and also of the spin in Wu & Murray 2003; Fabrycky & Tremaine 2007; Naoz et al. 2011; Correia et al. 2012), but most of these studies use averaged methods, they do not study the influence of an evolving host body radius, and they often consider only coplanar systems. We introduce here a tool which allows more complete studies. Indeed, the tidal equations used in this code are not averaged equations so it is possible to study phenomena such as resonance crossing or capture. Contrary to other codes using semiaveraged (Mardling & Lin 2002) or nonaveraged equations (Touma & Wisdom 1998; Mardling & Lin 2002; Laskar et al. 2004; Fienga et al. 2008; Beaugé & Nesvorný 2012; Correia & Robutel 2013; Makarov & Berghea 2014), our code will be freely accessible and open source.
After describing the tidal model used here, we seek to validate the code by comparing one planet evolutions around BDs with evolutions computed with a secular code (as used in Bolmont et al. 2011, 2012). We use systems around BDs to test systems for which tides are really strong and lead to important orbital changes. We then offer a glimpse of possible investigations possible with our code taking the example of the dynamical evolution of the Kepler62 system (Borucki et al. 2013).
2. Model description
The major addition to Mercury provided in MercuryT is the addition of the tidal forces and torques. But we also added the effect of general relativity and rotationinduced deformation. We explain in the following sections how these effects were incorporated in the code.
We also give the planets and star/BD/Jupiter parameters which are implemented in the code.
2.1. Tidal model
To compute the tidal interactions we used the tidal force as expressed in Mignard (1979), Hut (1981), Eggleton et al. (1998) and Leconte et al. (2010) for the constant time lag model. This model consists in assuming that the bodies considered are made of a weakly viscous fluid (Alexander 1973).
We added this force in the Nbody code Mercury (Chambers 1999). We consider the tidal forces between the star and the planets but we neglect the tidal interaction between planets. We consider here a population of N planets orbiting a star.
As in Hut (1981), we stop the development at the quadrupole order. At this order, we can use the point mass description of the tidal bulges. Star and planets are deformed. Due to the presence of planet j, the star of mass M_{⋆} is deformed and can be decomposed in a central mass M_{⋆} − 2μ_{⋆}, and 2 bulges of mass μ_{⋆}. As in Hut (1981), each bulge is located at a radius R_{⋆} from the center of the star and they are diametrically opposed. Figure 1 shows the geometrical setting of the problem. The central mass of the star is labeled by S, and the bulges by S′ and S′′. The mass of a bulge depends on the time lag and is given by: (1)where r_{j} is the distance between the star and planet j at time t − τ_{⋆}. R_{⋆} is the radius of the star, k_{2,⋆} its potential Love number of degree 2 and τ_{⋆} its constant time lag.
Due to the presence of the star, the planet j is deformed and can be decomposed in a central mass M_{pj} − 2μ_{pj}, and 2 bulges of mass μ_{pj}. The central mass of the planet j is labeled by P_{j}, and the bulges by P and P. The bulges mass is given by: (2)where R_{pj} is the radius of planet j, k_{2,pj} its potential Love number of degree 2 and τ_{pj} its time lag. To the lowest order in τ_{pj}: (3)Up at the third order in R_{pj}/r_{j} and R_{⋆}/r_{j} the forces exerted by the primary on the secondary are the following gravitational forces: f_{S → Pj}, , , f_{S′ → Pj} and f_{S′′ → Pj}, where the latter expression is given by: (4)where P_{j}S′′ is the vector , defined in Fig. 1.
Fig. 1 Twodimensional schematic representing the two deformed bodies. The star is divided into three masses: the central one of mass M_{⋆} − 2μ_{⋆} at S, and two bulges of mass μ_{⋆} at S′ and S′′. Planet j is divided into three masses: the central one of mass M_{pj} − 2μ_{pj} at P_{j}, and two bulges of mass μ_{pj} at and . Ω_{⋆} is the star rotation vector (its norm is Ω_{⋆}, the star rotation frequency), Ω_{pj} is planet j rotation vector (its norm is Ω_{pj}, the planet rotation frequency), and is a vector collinear with the orbital angular momentum of planet j (its norm is equal to the derivative of the true anomaly). e_{rj} is the radial vector. 

Open with DEXTER 
Let us define F_{tr} (for tides radial), P_{to,⋆} and P_{to,pj} (for tides orthoradial): (5)F_{tr} has the dimension of a force (M.L.T^{2}), while P_{to,pj} and P_{to,⋆} have a dimension of a momentum (M.L.T^{1}).
Then the total resulting force due to the tides acting on planet j is: (6)where Ω_{⋆} is the star rotation vector, Ω_{pj} is planet j rotation vector and v_{j} = ṙ_{j} is the velocity of planet j. e_{rj} is the unit vector defined as: SP_{j}/ SP_{j}. is a vector collinear with the orbital angular momentum of planet j (L_{horb} defined hereafter), whose norm is equal to the derivative of the true anomaly. The term can be rewritten as follows: (7)What modifies the spin of the star is the following torque contribution . To calculate the torque on the star, , we consider the planet as a point mass, meaning that we neglect the gravitational interaction between the bulges of the planet and the bulges of the star. Likewise, the torque contribution which modifies the spin of planet j, , is . So the torque exerted by planet j on the star is given by: (8)and the torque exerted by the star on the planet j is: (9)With this description of the phenomenon, we consider that each planet creates an independent tidal bulge on the star, and that the bulge created by planet j doesn’t affect planet i ≠ j.
2.2. General relativity
We added to Mercury the force due to general relativity as given in Kidder (1995), Mardling & Lin (2002). This force corresponds to the orbital acceleration due to the postNewtonian potential and is given by: (10)F_{GRr} and F_{GRo} are given by: (11)where v_{j} is the norm of the velocity v_{j} of the planet and c is the speed of light and η is given by: (12)
2.3. Rotational deformation
The equilibrium figure of a viscous body in rotation is an triaxial ellipsoid symmetric with respect to the rotation axis (Murray & Dermott 1999). The rotational deformation is quantified by the parameter J_{2}, defined as follows for planet j: (13)and for the star: (14)where k_{2f,pj} is the fluid Love number of planet j and k_{2f,⋆} that of the star. The fluid Love number is here defined as the potential Love number for a perfectly fluid planet (see for example Fig. 2 of Correia & Rodríguez 2013, for the Earth’s potential Love number and fluid Love number). Our code allows the user to choose different values for the fluid Love number k_{2f,p} and the potential Love number k_{2,p}.
The total resulting force due to the rotational deformation of star and planet j on planet j is (Murray & Dermott 1999; Correia et al. 2011): (15)where C_{⋆} and C_{pj} are defined as follows: (16)The torque exerted by planet j on the star is given by: (17)and the torque exerted by the star on the planet is: (18)These torques are responsible for the precession of the orbit normal. This precession has an influence on the mean eccentricity of the planets and also on their mean obliquity.
2.4. Summing all the effects
2.4.1. Corrective acceleration
In order to compute the evolution of the planetary orbits, we have to take into account all these effects. The orbital part of the acceleration is handled by Mercury, so we give here the expression of the corrective acceleration of planet j in the astrocentric coordinates: (19)where , and are defined in the previous sections.
2.4.2. Spin equations
Let us consider first a system with one planet. We make the hypothesis that we can decouple the torque equation given by the conservation of total angular momentum L: (20)where I_{⋆} is the principal moment of inertia of the star, and I_{p} that of the planet. L_{orb} is the orbital angular momentum: (21)where r^{∗} and v^{∗} are the position and velocity of the planet in the reference frame of the center of mass of the system. and are the position and velocity of the star in the reference frame of the center of mass of the system.
In astrocentric coordinates, r and v: (22)so, (23)where N_{⋆ → p} is the total torque exerted on the planet and N_{p → ⋆} is the total torque exerted on the star.
Assuming that the spins of the star and the planet evolve solely as a result of tidal and rotational flattening torques, this means that we can decouple Eq. (20) to obtain the following spin equations: (24)If there are more than one planet in the system, the orbital angular momentum involves planet j– planet i≠j cross terms (e.g., in r_{j} × v_{i}, r_{i} × v_{j}) that we neglect for the spin calculation. The equations governing the evolution of the spin of the star and the spin of planet j are therefore: (25)where , , and are defined in the previous sections.
Our MercuryT code provides the x, y and z components of the spin of the planets, of their orbital angular momentum, and of the spin of the star. Then the obliquity ϵ_{pj} and the inclination i_{j} of planet j are obtained by calculating: (26)where L_{orbj} is a vector normal to the orbit of planet j.
The spin of planet j is the norm of the vector Ω_{pj} and the spin of the star is the norm of the vector Ω_{⋆}.
2.5. Integration of the spin
In this work, we use the Mercury code’s hybrid routine, which relies on a hamiltonian description of the problem. The hamiltonian is divided in three parts that are integrated consecutively (Chambers 1999). Mercury allows the user to add other forces on the planets via a routine. In the hamiltonian description, these extra forces are treated as a perturbation to the Keplerian potential. The user should keep in mind that this violates the symplectic properties of the integrator.
In this user routine, we added the tidal forces, rotationinduced flattening resulting forces, GR forces. On the one hand, we compute the resulting acceleration on the planets, and this acceleration is used by the bulk of Mercury to compute the evolution of the system. On the second hand, the spin of the planets and the star is integrated with a RungeKutta of order 5 within the routine. The RungeKutta integration is performed twice in a Mercury time step. To compute the spins at a time t, the positions and velocities of the planets are interpolated between t −dt and t at the time intervals needed for the RungeKutta routine. The integration scheme is illustrated in Fig. 2.
Fig. 2 Integration scheme of MercuryT. 

Open with DEXTER 
If the host body evolves, the moment of inertia of the star, given by: I_{⋆} = M_{⋆}(rg_{⋆}R_{⋆})^{2}, where rg_{⋆} is the radius of gyration (Hut 1981) is going to vary with time.
The equation of each component of the spin Ω_{⋆} is given in the following equation, written here for the z component Ω_{⋆,z}: (27)The error brought by the spin integration depends on the third power of the time step (Chambers 1999). The part of the integration that causes the more errors is the integration of the rotation induced flattening effect. The integration of the tidal torque does not require such a precise integration due to long timescales of evolution with respect to the Mercury time step (which is usually taken slightly smaller as one tenth of the inner planet’s orbital period). However, the rotation induced flattening causes changes on the spin of the planets on a much shorter timescale.
For very closein planets, the integration of this part can lead to a purely numerical diminution of the rotation period. It is therefore important to evaluate the error made during the integration of the rotationinduced flattening torque. We advise to use a time step smaller than the orbital period of the inner planet divided by 20. The time step should be chosen according to the precision needed on the rotation period of the inner planet (see Sect. 3.5 for details).
2.6. Input parameters
Planetary parameters implemented in MercuryT.
Host body parameters implemented in MercuryT.
In order to work the code requires the necessary planetary parameters which are used in all the equations of Sects. 2.1−2.3. Any Nbody integrator needs parameters such as: the masses of the planets M_{pj}, the mass of the host body M_{⋆}, the semimajor axis (SMA) of the planets, their eccentricities (ecc), their inclination (inc) and their orbital angles (argument of pericentre, longitude of the ascending node, mean anomaly). In the following tests, the orbital angles are set to 0°.
In order to calculate the evolution due to rotational flattening, one needs the fluid Love number of the star k_{2f,⋆} and of the planets k_{2,pj}.
In order to calculate the tidal evolution, one needs also the radius of the star R_{⋆} and the radius of the planets R_{pj}, the potential Love number of degree 2 of the star k_{2,⋆} and of the planets k_{2,pj}, the time lag of the star τ_{⋆} and the time lag of the planets τ_{pj}.
All of these parameters are of course changeable in our code, but we implemented some useful values and relations for commodity. These implemented values are given in Tables 1 and 2.
2.6.1. Planets model
Values of the reduced dissipation factor for the host bodies implemented in MercuryT.
Test simulations for tides: one BD and one planet.
For Earthlike planets and SuperEarths, we assume that the product of the potential Love number of degree 2 and the time lag of the planet is equal to that of the Earth: k_{2,p}Δτ_{p} = k_{2, ⊕}Δτ_{⊕}. We use the value of k_{2, ⊕}Δτ_{⊕} = 213 s given by Neron de Surgy & Laskar (1997). We assume here that the fluid Love number and the potential Love number of degree 2 are equal.
Given the mass of the planet, we offer the user two possibilities to choose the radius: either he gives a value himself or he assumes a composition and the code calculates the radius following Fortney et al. (2007). For example, a superEarth of 10 M_{⊕} would have a radius of 1.8 R_{⊕}.
For the Jupiterlike planet, we computed the time lag τ_{HJ} from the value of the dissipation parameter σ_{k} for Hot Jupiters of Hansen (2010). The notation σ_{k} was introduced by Eggleton et al. (1998) and is linked to the quantity k_{2,k}τ_{k} by: (28)where k represents either the star or a planet.
2.6.2. Host body evolution and dissipation
It is possible to use stellar evolution tracks in MercuryT to compute the evolution of planets around an evolving object. We implemented this for an evolving Jupiter (Leconte & Chabrier 2013), for evolving BDs of masses: 0.01, 0.012, 0.015, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.072, 0.075 and 0.08 M_{⊙} (Leconte et al. 2011), for a Mdwarf (dM) of mass 0.1 M_{⊙} and a Sunlike star (Chabrier & Baraffe 1997; Baraffe et al. 1998).
Table 2 shows for each type of host body which are the evolving quantities and which have implemented values. The evolving quantities are tabulated and MercuryT interpoles the values during the integration in order to have the correct radius/Love number/moment of inertia in the acceleration formula and the spin equations.
We make here the hypothesis that during their evolution the dissipation factor σ_{k} of the host body remains constant. We use the value of the dissipation of Jupiter given in Leconte et al. (2010) for the Jupiter host body. We use the dissipation factor of Bolmont et al. (2011) for BDs and for the Mdwarf, and we use Hansen (2010)’s value of stellar dissipation for the dissipation of the Sunlike star.
In Table 3, we give as an indication the value of the parameter for the different bodies, where is defined as: (29)and where is the freefall time at the surface of the star. As this definition of depends on the radius of the star, we compute its value for all the bodies of Table 2 for an age of 1 Myr and of 1 Gyr.
Test simulation for rotational flattening: one BD with 2 planets.
3. Code verification
In order to validate the tidal part of the code, we first simulated the tidal evolution of one Earthmass planet orbiting a 0.08 M_{⊙} BD with two different approaches. The first approach is to use a secular code, which solves the averaged equations of the tidal evolution of one planet (equations in semimajor axis, eccentricity, etc. often used in tidal studies such as Hut 1981; Leconte et al. 2010; Bolmont et al. 2011). The second simulation was done with the MercuryT code we developed. We first compared the outcomes for different regimes, switching on or off the planetary tide (i.e., the tide raised by the BD in the planet) and the BD tide (i.e., the tide raised by the planet in the BD) and testing the effect of the evolving radius of the BD.
The details of all simulations are listed in Table 4, where dt is the time step used for the simulation.
In order to test the rotational flattening part of the code, we compared our results with the numerical code used in Correia & Robutel (2013), hereafter denoted by “CR13 code”. This code has been developed independently of the present one and uses the ODEX integrator (e.g. Hairer et al. 1993). In Correia & Robutel (2013) the CR13 code was applied to a concrete situation, the spin evolution of trojan bodies, but it is much more general than that. It has not be made available for public use, but it has the ability to perform the same kind of simulations of the MercuryT code. Therefore, the CR13 code will be used for crosschecking some of our results.
We considered a 2planet system orbiting a 0.08 M_{⊙} BD and validated the effect of the rotational flattening of the star, of the planet and compared our results for a full simulation (effects of tides and rotational flattening). The details of the simulations are listed in Table 5, where k_{2,p1 / 2} is the Love number of degree 2 of planet 1 and planet 2, k_{2,⋆} is the Love number of degree 2 of the host body.
3.1. Nonevolving BD, effect of planetary tide
This case corresponds to case (1) of Table 4, for which we switched off the effect of the BD tide.
Fig. 3 Case (1): tidal evolution of a planet of mass 1 M_{⊕} orbiting a 0.08 M_{⊙} BD calculated with the secular code (blue dashed line) and the MercuryT code (full red line). Graph a) from top to bottom: evolution of the semimajor axis of the planet, evolution of its eccentricity and evolution of the tidal heat flux. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its rotation period (the pseudosynchronization period is represented in long red dashes), and conservation of total angular momentum. 

Open with DEXTER 
Figure 3 shows the evolution of the semimajor axis, the eccentricity and the averaged tidal heat flux ⟨ φ_{tides} ⟩ defined as follows: (30)where ⟨ Ė_{tides} ⟩ is the averaged gravitational energy lost by the system by dissipation: (31)where T_{p} is the dissipation timescale. Na1(e), Na2(e) and Ω(e) are eccentricity dependent factors defined in Bolmont et al. (2013). The tidal heat flux depends on the eccentricity and on the obliquity of the planet. If the planet has no obliquity, no eccentricity and if its rotation is synchronized, the tidal heat flux is zero.
We also show the evolution of the instantaneous tidal heat flux, computed from the instantaneous energy loss given by: (32)where is the derivative of the spin of planet j, given by Eq. (25). Contrary to ⟨ Ė_{tides} ⟩ which depends on averaged computed values like the semimajor axis and the eccentricity, Ė_{tides}(t) depends on instantaneous position, velocity and spin of planet j.
The eccentricity of the planet decreases to values below 10^{4} in 10^{7} yr. The decrease of the eccentricity of the planet is accompanied by a decrease of the semimajor axis. The evolution of these two quantities show a good agreement between the secular code and the MercuryT code. The obliquity of the planet decreases from its initial value of 11.5° to less than 10^{4} degrees in less than 500 yr. During the same time, the rotation period evolves from its initial value of 24 h to the pseudosynchronization period, which is here ~48.5 h. The evolution of obliquity and rotation period show a good agreement between the secular code and MercuryT.
After 2 × 10^{7} yr of evolution, the eccentricity obtained with MercuryT is equal to a few 10^{7}. This residual value of the eccentricity comes from the way the Mercury code calculates the orbital elements from the positions and velocities of the planets. It assumes a Keplerian potential however in this situation where the tidal forces are taken into account this is not the case.
In any case, we can assume that an eccentricity of 10^{7} can be considered as null. Furthermore, this code is designed to study multiplanet systems. In the examples we will give later, the eccentricity due to planetplanet interactions is typically bigger than 10^{7}.
This residual eccentricity is responsible for a non zero averaged tidal heat flux ≲ 10^{2} W/m^{2}. However, the instantaneous tidal heat flux reaches values as low as a few 10^{9} W/m^{2}. This low value illustrates the fact the “real” eccentricity of the planet must be much lower than what MercuryT is calculating.
Fig. 4 Case (3): tidal evolution of a planet of mass 1 M_{⊕} orbiting a 0.08 M_{⊙} BD calculated with the secular code (blue dashed line) and the MercuryTcode (full red line). Graph a) from top to bottom: evolution of semimajor axis (in red) and of the corotation distance (in black), evolution of eccentricity and evolution of tidal heat flux. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its inclination and evolution of its rotation period (in red) and the BD rotation period (in black). The pseudosynchronization period is represented in long red dashes. 

Open with DEXTER 
We also verify that each component of the total angular momentum is a conserved quantity during the evolution of the system (Eq. (20)). We define the quantity α_{i} where i is x, y or z, and α as: (33)where L_{i} is the icomponent of the total angular momentum vector and L is the norm of the vector L of Eq. (20). In this example, we only consider the effect of the planetary tide, which is equivalent to consider the BD as a point mass, so that the total angular momentum is here only the sum of the orbital angular momentum and the rotational angular momentum of the planet.
Right bottom panel of Fig. 3 shows the conservation of the total angular momentum as a function of time. For the MercuryT simulation, each component of the total angular momentum α_{i} is conserved and α reaches 10^{6} after 100 Myr of evolution. For the secular code, the total angular momentum is a little less well conserved and reaches a few 10^{5} at the end of the simulation. In this example, the orbital angular momentum is 10^{5} higher than the angular momentum of the planet, so the question remains of knowing if the spin of the planet is correctly computed.
In order to test this, we performed another simulation – case (1′) – with a Jupitermass planet to reduce the difference between the orbital angular momentum and the planet’s angular momentum. In this case, the orbital angular momentum is about 10^{3} higher than the angular momentum of the planet. We find that for this simulation the total angular momentum is conserved, with α asymptoticly reaching only 6 × 10^{6} after a 10 Myr evolution^{3}.
3.2. Nonevolving BD, effect of BD tide
This case corresponds to case (2) of Table 4, for we switched off the effect of the planetary tide and considered a very dissipative BD. The planet is initially outside the corotation radius. The initial orbital distance being larger than the corotation distance, the planet migrates outward.
In agreement with the secular code, the BD tide makes the inclination of the planet decrease from ~12° to ~4.5° in 10^{8} yr. When the planet migrates outward, the rotation period of the BD increases in agreement with the conservation of total angular momentum. We find that each of the component of the total angular momentum is conserved. In this example, the angular momentum of the BD is two orders of magnitude higher than the orbital angular momentum of the planet. As α remains below 10^{4} after 100 Myr, we can conclude here that the phenomenon is accurately reproduced.
3.3. Nonevolving BD, effect of both tides
This case corresponds to case (3) of Table 4. The planet is initially outside the corotation radius.
Fig. 5 Case (3′): tidal evolution of a Jupitermass planet orbiting a 0.08 M_{⊙} BD calculated with the secular code (blue dashed line) and the MercuryT code (full red line). Graph a) from top to bottom: evolution of semimajor axis (in red) and of the corotation distance (in black), evolution of eccentricity and conservation of angular momentum. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its inclination and evolution of its rotation period (in red) and the BD rotation period (in black). The pseudosynchronization period is represented in long red dashes. 

Open with DEXTER 
Figure 4 shows the evolution of this system. We find that the results of MercuryT are in good agreement with the secular code. All quantities vary accordingly and the quantitative agreement is very good. As the planet migrates away, the evolution timescales become longer which entails a slower evolution at late ages of the eccentricity and inclination more particularly.
The initial heat flux is very strong in comparison to the tidal heat fluxes measured for solar system bodies: 0.08 W/m^{2} for Earth (Pollack et al. 1993) and between 2.4 and 4.8 W/m^{2} for Io (Spencer et al. 2000). Such an important heat flux must have repercussions on the planet’s internal structure. The high fluxes of Fig. 4 must imply that the surface and the interior of the planet are melted and that the vertical heat transfer is very efficient, which must not be in agreement with the dissipation factor value used here. We here neglect any feedback of the dissipation on the internal structure (as in Bolmont et al. 2013).
For this system, each components of the total angular momentum α_{i} is conserved and α reaches a few 10^{7} at a time of 10 Myr for the MercuryT simulation. The total angular momentum is therefore here also conserved.
We also test the strength of our code with a more extreme case: case (3′). With an initial orbital distance of 9 × 10^{3} AU, the planet is initially inside the corotation radius and thus migrates inward. We consider here a BD of a low dissipation factor (0.01 × σ_{BD}) so that the planet’s BDtide driven inward migration is not too quick.
Figure 5 shows the evolution of this system. The planet falls onto the BD in about 3000 yr. During the fall eccentricity, obliquity and inclination decrease. In less than 2000 yr, the rotation period of the planet evolves from 240 h to the pseudosynchronization period (of about 24 h). The rotation period of the BD decreases just prior to the fall due to the angular momentum transfer from the planet’s orbit to the BD spin.
Even for this extreme case, both codes lead to the same simulated evolution. The collision time is slightly different, but is of the same order of magnitude for both simulations. Left bottom panel of Fig. 5 shows the conservation of total angular momentum for this example.
As the planet ends up by colliding with the BD, we do not expect the conservation of total angular momentum to be perfect. Indeed, Fig. 5 shows that for both simulations, α increases with time and reaches about 10^{2} when the collision occurs (4 × 10^{3} for MercuryT). The planet is initially very close to the BD and gets closer in time so tidal effects are stronger and stronger. This example allows to test the limits of our model. For closein planets, one should always verify that α is conserved.
However, the destiny of the planet is compatible with the theory. As its initial orbital distance is lower than the corotation distance, the BD tide acts to push the planet inward. The qualitative evolution is not likely to change if the code was improved, however the time of collision between the planet and the BD could change.
3.4. Evolving BD, effect of both tides
This case corresponds to case (4) of Table 4, for which we consider the evolution of the radius and radius of gyration of the BD.
Our code allows to choose the initial time of the simulations, i.e. the BD age from which we consider the tidal evolution of the planets. We consider in this work (as in Bolmont et al. 2011, where they discuss the influence of this initial time) that the initial time corresponds to the time of the dispersal of the gas protoplanetary disk. We then consider that the planets are fully formed by this time. The time indicated on the figures corresponds to the time spent after this initial time.
Fig. 6 Case (4): tidal evolution of an Earthmass planet orbiting a 0.08 M_{⊙} BD calculated with the secular code (blue dashed line) and the MercuryT code (full red line). Graph a) from top to bottom: evolution of semimajor axis (in red) and of the corotation distance (in black), evolution of eccentricity and evolution of tidal heat flux. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its rotation period (in red), the BD rotation period (in black), and the pseudosynchronization period (in long red dashes), and evolution of α. 

Open with DEXTER 
Figure 6 shows the evolution of this system. The evolution calculated with the secular code is in good agreement with the evolution calculated with MercuryT. The competition between the outward migration due to the BD tide and the inward migration due to the planetary tide is well reproduced.
The comparison between the outcomes of the two codes shows a difference of 10^{5} AU in the calculated semimajor axis when the migration direction changes. This difference remains small and tends to decrease when the precision of the secular code is increased, which shows once more that the precision of the MercuryT seems better than the secular code.
In any case, the qualitative behavior is very well reproduced even though small quantitative differences can be seen. The MercuryT code reproduces well the evolution of the spin of the BD due to the contraction of its radius (middle panel of the graph b) of Fig. 6).
The total angular momentum is well conserved as can be seen in Fig. 6. Each component of the total angular momentum as well as α remain below 3 × 10^{6}.
Those different tests show that the tidal integration part of the MercuryT code shows a good agreement with the secular code concerning the orbital evolution of the planet as well as its rotation state evolution and that of the BD. The total angular momentum is always conserved, except when the planet collides with the BD. We therefore consider that this code is valid to be used to study the evolution of tidally evolving multiplanet systems. For any simulation, one should however always make sure that the total angular momentum is conserved.
3.5. Effect of the rotational induced flattening
Where an Euler integration of the spin was enough to correctly describe the tidal evolution of the spin of the planets, we needed to implement a better integrator to accurately describe the precession of the spin axis of the planet due to its own flattening. This precession happens on much shorter timescale than the tidal evolution. For the example of case (5), the timescales are about a few 10^{1} yr.
Thus, in order to have an accurate integration with MercuryT, we need to perform a integration RungeKutta of order 5 twice in a Mercury time step (dividing the timestep in 2 inside one Mercury time step allows us to increase the precision without being too time demanding).
We tested the integration of the rotational induced flattening by comparing our code to the CR13 code. We obtain similar results for all cases with some small quantitative differences.
For the cases (6) of Table 5, only the rotation flattening of the inner planet is taken into account. The rotation period P_{p} (i.e., the norm of the spin) is not influenced by the effect of the rotational flattening. However due to the integration scheme, we observe a small drift of the rotation period of the inner planet (Fig. 7). This drift increases linearly with time and decreases when the time step is reduced. For a time step of 0.08 day, case (6), the drift is of about 8 × 10^{5} h after 100 000 yr of evolution, while for a time step of 0.01 day, case (6′′), it is less than 10^{6} h. For time steps lower than 0.01, the drift is essentially null for the 100 000 yr of evolution.
The smaller the time step, the smaller the differences. From a time step of 0.08 day to 0.01 day the improvement is visible, but we can see that there is almost no difference between the light and dark blue curves corresponding to a time step of 0.01 and 0.001 day. Of course, the execution time is longer for smaller time step, so the time step should be chosen according to the duration of the simulation and the precision needed for the spin of the inner planet of the system.
Fig. 7 Evolution of (P_{p} − P_{p,0}) /P_{p,0} of the inner planet of the system corresponding to cases (6) of Table 5. The colored lines correspond to the results of MercuryT: red with a time step of 0.08 day, green with a time step of 0.05 day, light blue with a time step of 0.01 and dark blue with a time step of 0.001 day. 

Open with DEXTER 
Fig. 8 Evolution of the obliquity of the inner planet of the system corresponding to case (6′′) of Table 5 for the last 600 yr of the evolution. The black line corresponds to the results of the CR13 code. The red line corresponds to the results of MercuryT. 

Open with DEXTER 
Stellar properties.
The semimajor axis of the planets show a perfect agreement, while the eccentricity, the obliquity, the rotation period, and to a lesser extend the inclination show a small difference of oscillation frequency.
Planetary physical parameters.
Figure 8 shows the evolution of the obliquity of the inner planet of the system corresponding to the case (6′′) of Table 5 compared with the CR13 code. The mean value of the obliquity, as well as the maximum and minimum values and well reproduced, the only thing different is the oscillation frequency.
The remaining small difference between MercuryT and the CR13 code – i.e., the rotation period drift and the difference in the frequency of the oscillations of quantities – come from the different integration schemes and numerical effects.
The system is also very sensitive to the initial conditions. For example, for a similar rotation period and obliquity, changing the initial direction of the spin of the planet Ω_{p} leads to a different mean value of the oscillations of the obliquity.
For case (7), for which all the effects are considered, the general behavior is perfectly reproduced. Both planets migrate outward due to the BD tide and enter a mean motion resonance approximatively at the same moment. Entering the resonance, eccentricities and obliquities evolve similarly.
We tested as well the conservation of energy and total angular momentum for these examples. For case (6), we find that each component of the total angular momentum is conserved and α reaches a few 10^{7} after 10^{5} yr of evolution. The total energy of the system is conserved up to a few 10^{5}. For all the other cases (6′) to (6^{′′′}), the conservation is slightly better but the orders of magnitude are the same.
In the case of a nondissipative force, our code conserves total energy and angular momentum.
Apart from very small differences, MercuryT and the CR13 code give the same results. We consider this agreement good enough for the study of exoplanets.
4. The case of Kepler62
MercuryT can be used to study hypothetical systems but it can also be used for known exoplanet systems (this code has been used in the following articles: Bolmont et al. 2013, 2014b; Quintana et al. 2014; Heller et al. 2014). We present here a study about Kepler62, which is a system hosting 5 planets (Borucki et al. 2013). Two of these planets are in the insolation habitable zone (HZ).
Planetary climate depends on many different parameters such as orbital distance, eccentricity, obliquity, rotation period and tidal heating (Milankovitch 1941; Spiegel et al. 2009, 2010; Dressing et al. 2010). All these parameters are influenced by tidal interactions so it is paramount to take into account tides in climate studies (Bolmont et al. 2014a).
We present here a dynamical study of the Kepler62 system, showing the influence of the different physical effects – tides, rotation flattening, general relativity – on the stability of the system and on the spin state evolution of the planets.
We used as initial conditions the data from Borucki et al. (2013) for semimajor axis, eccentricity, longitude of periapsis and epoch of mid transit. We used the values given for the radius of the planets, and we tested the system with different masses for the planets (which we all considered rocky).
The simulations we show here are of course possible evolutions of the system, we are aware that there are many uncertainties on many parameters, starting with the masses of the planets and their dissipation factors. However, despite those uncertainties, we aim here at showing that some general behavior can be identified in the dynamics of the system.
4.1. Dynamics and stabilization
In order to investigate the effects of tides, rotationflattening and general relativity on the dynamics of the system, we tested the system for five different cases, which are listed in Table 8.
Test simulations for stability.
Assuming the planets have a rocky composition (in this hypothesis, planet d has the maximum mass given in Borucki et al. 2013), we find that the system – hereafter called system – is unstable in case (1). After 3 Myr, planet c is ejected from the system. In this case, planet d is very massive: 14 M_{⊕}, its influence on the less massive planet c destabilizes the system. We find that system is also unstable in case (2). However, the destabilization occurs much later, after ~20 Myr of evolution. The correction for general relativity has here the effect of stabilizing the system. General relativity causes apsidal advance, this can therefore lead to situations favorable to stability or unfavorable according to initial conditions. Changing the initial orbital angles such as the longitude of periastron modify the amplitude of the eccentricity oscillations and could also lead to a more or less stable system. However, for our particular choice of initial conditions, it would seem that general relativity has here a stabilizing effect.
Using the same masses for planets b, c, e and f as in system , we tested the stability of the system in cases (1) and (2) for different masses of planet d. We found that the system is systematically stable for masses of planet d lower than ~8 M_{⊕}. However for masses bigger than ~8 M_{⊕}, most simulations lead to destabilization within 30 Myr. For simulations done with a mass of planet d higher than 8 M_{⊕}, destabilization occurs either for cases (1) or (2) or both illustrating the importance of taking into account the correction for general relativity.
We also observed that changing the masses of the planet very slightly can influence the stability of the system. Changing the masses of 5% leads to a stable system in all cases – hereafter called system ℬ. A broad study of the stability of this system is out of the scope of this paper, where we want to illustrate the possible studies allowed by the use of MercuryT. The stability should be tested for all possible masses, which leads a high number of combinations and simulations in order to have a map of the stability of the system (e.g., Laskar 1990; Correia et al. 2005; Couetdic et al. 2010; Mahajan & Wu 2014). Unstable regions would therefore correspond to unrealistic configurations. This illustrates the importance of constraining the masses of planetary systems (e.g. with HARPS in Dumusque et al. 2014).
Fig. 9 Tidal evolution of the Kepler62 system (ℬ). Top panel: evolution of the obliquities of the five planets. Bottom panel: evolution of their rotation periods in colored full lines. The dashed lines correspond to the pseudosynchronous rotation period and the full black line corresponds to the rotation period of the star. 

Open with DEXTER 
When we add tides – case (3) – and assume nominal dissipation factors for the planets, we find that system becomes stable for the duration of the simulation. In this case, tides have a stabilizing effect on the system. Indeed, in our simulations, both planetary tide and stellar tide have a damping effect on the eccentricity, therefore decreasing the probability of the system to be chaotic. System ℬ remains stable at least for the duration of the simulation of 30 Myr. Semimajor axes and eccentricities do not significantly tidally evolve during the simulation, however the obliquities and rotation period of the planets do (see Fig. 9).
We tested the evolution of system ℬ for different planetary dissipation factors. The higher the dissipation of planet j, the faster its tidal evolution. The effect is first visible on the rotation of the planets (obliquity and rotation period), which evolve much more rapidly. It also has a small effect on the eccentricity of the planets, though not visible on the graphs. In order to quantify it, we computed the mean angular momentum deficit (AMD, e.g. Laskar 1997) for a set of simulations.
We did a test varying the dissipation of each planet from a reference simulation corresponding to a dissipation factor of 0.1σ_{⊕} for all planets. We consecutively increased the dissipation of each planet from the reference value of 0.1σ_{⊕} to 10σ_{⊕} and 100σ_{⊕}.
Increasing the dissipation of a planet makes the AMD decrease. The decrease is more or less pronounced depending on the planet considered. When increasing the dissipation of planets c, e and f to 100σ_{⊕}, the AMD is ~0.03% lower than the reference case. However, when increasing the dissipation of planet d, the AMD is 0.055% lower. As planet d is very massive in the system, damping its eccentricity slightly has consequences on the whole system. When increasing the dissipation of planet b to 100σ_{⊕}, the AMD is ~0.7% lower than the reference case. Increasing the dissipation of the closest planet has the biggest effect on the dynamics of the system. Increasing the dissipation leads to a slightly less chaotic system.
When we add the effect of rotationflattening – case (4) – we find that system is stable at least for the duration of the simulation of 30 Myr. This effect stabilizes the system by changing the precession rates. The dynamics of system ℬ are not significantly changed by this effect.
When we consider all effects – case (5) – we find that system is stable at least for the duration of the simulation of 30 Myr. Compared to case (3), the addition of the rotationflattening effect only changes slightly the equilibrium values of the obliquities of the planets.
4.2. Obliquity and rotation period
Due to the planetary tide, the obliquity decreases and the rotation period evolves towards pseudosynchronization. The evolution timescales of these two quantities are shorter than the timescales of evolution of semimajor axis and eccentricity. Figure 9 shows that for Kepler62 the rotation period of the three inner planets of the system evolves towards pseudosynchronization in less than ten million years and their obliquities evolve towards small equilibrium values (< 1°).
These simulations were performed using as initial conditions the values in Borucki et al. (2013), so that our simulation shows what could be the evolution of the system in the future. However, we can draw some conclusions on the past evolution of the system from Fig. 9, or a simple evolution timescale calculation. Given as the age of the system is estimated at 7 Gyr (Borucki et al. 2013), we indeed expect that the three inner planets of the Kepler62 system are now slowly rotating (their period is higher than 100 h) and they have quasi null obliquities.
The ratio between the pseudo synchronization rate and the orbital frequency only depends on the eccentricity of the planet. But in a multiplanet system, the eccentricity of a planet is excited due to planetplanet interactions, and oscillates with a combination of frequencies that correspond to secular modes (e.g., Murray & Dermott 1999). In our simulations, the planets experience relatively big eccentricity oscillations (see next section), so that the planets’ pseudosynchronization periods also oscillate. In reality, the rotation periods of the planets are not exactly equal to the corresponding pseudosynchronization period. Figure 10 shows the evolution of the rotation period of Kepler62b compared to the pseudosynchronization period and the synchronization period. The pseudosynchronization period oscillates too fast for the rotation period to be able to follow. The instantaneous rotation period of Kepler62b oscillates out of phase with the pseudosynchronization period and with a lower amplitude.
Fig. 10 Shortterm (100 000 yr) evolution of the rotation period of Kepler62b (ℬ). Full line: rotation period, dashed line: pseudosynchronization period, dasheddotted line: synchronization period. 

Open with DEXTER 
During the 30 million yr of the simulation, the obliquities and rotation periods of the HZ planets Kepler62e and f do not evolve significantly. So we performed longer simulations for these two outer planets. Assuming an Earthlike dissipation for the two planets, we found that Kepler62e is likely today to have reached pseudosynchronization and have low obliquity. Figure 11 shows that after 3 Gyr of evolution, the obliquity has been damped and the rotation pseudosynchronized. For Kepler62f the timescales of evolution are higher and Fig. 11 shows that the rotation period is still evolving towards pseudosynchronization after 7 Gyr of evolution and that the obliquity can still be high.
The dissipation of the planets is not constrained and changing the dissipation would only shift the curves right (if the dissipation is lower) or left (if the dissipation is higher). As an Earthlike dissipation is actually probably a high dissipation value (due to the presence of oceans, e.g. Lambeck 1977), it is likely that the curves should be shifted to the right.
Fig. 11 Longterm tidal evolution of the two outer planets of the Kepler62 system. Top panel: evolution of the obliquities of the five planets. Bottom panel: evolution of their rotation periods in colored full lines. The dashed lines correspond to the pseudosynchronous rotation period and the full black line corresponds to the corotation radius. 

Open with DEXTER 
4.3. Consequence of dynamics on the potential habitability of Kepler62 e and f
Kepler62e and Kepler62f are both inside the HZ, however being in the HZ does not entail presence of surface liquid water. Surface conditions compatible with liquid water depend of course on the properties of the atmosphere considered (pressure, temperature and chemical composition) but also on orbital parameters such as semimajor axis and eccentricity, and physical parameters such as the obliquity and the rotation period of the planet (Milankovitch 1941).
Thanks to the MercuryT code, we can simulate the dynamical evolution of habitable planets within their system, taking into account the rich dynamics occurring in a multiplanet system. This allows us to provide to any kind of climate model a set of orbital and physical input parameters consistent with the real dynamics of a planetary system.
Figure 12 shows the evolution of the eccentricity of the planets for 100 000 yr. The eccentricities of Kepler62e and f oscillate respectively between 0.02 and 0.16 and between 0.05 and 0.19 with a modulated frequency. These important periodical changes in eccentricity have an effect on the climate on Kepler62e and f, such as the Milankovitch cycles had an impact on the paleoclimate of Earth (Berger et al. 1992).
Furthermore, as seen in the previous section, we can also draw some conclusions about the rotation of the planets. We have shown that it is likely that Kepler62e has a pseudosynchronous rotation (or a near pseudosynchronous rotation, as shown in Sect. 4.2) and that its obliquity is very small. Such a planet with a slow rotation (almost 3000 h or 125 day) could have large Hadley cells that would bring hot air to the poles and there would be no longitudinal circulation (Merlis & Schneider 2010; Leconte et al. 2013). However, as Kepler62e is close to the inner boundary of the HZ, the surface temperatures might not reach such low values as to allow cold traps. A study of this planet with a Global Circulation Model would be needed to test the potential of this planet to host surface liquid water.
Figure 11 shows that Kepler62f can have a high obliquity and a fast rotation period (as fast as Earth’s 24 h rotation). Of course, we don’t know the initial conditions on the spin of the planets. However, formation scenarios showed that planets are likely to have an initial fast rotation rate due collisions and have an isotropic distribution of obliquities (Kokubo & Ida 2007). For fast rotation rates, the obliquity is excited and can reach high values even if it started at a small one (see for example the red full line in Fig. 11), so there is therefore a high probability that the obliquity of Kepler62f is actually non negligible. Furthermore, its rotation period could be still quite fast: between 20 and 40 h at the supposed age of the system (see Fig. 11). It is likely that Kepler62f would have a very different type of climate than its neighbor. A nonnegligible obliquity means seasonal effects and a fast rotation means a different wind pattern with not only latitudinal winds but longitudinal winds also.
5. Conclusions
We presented here a code that computes orbital evolution for tidally evolving multiplanet systems. The theory on which this code is based is the constant time lag model, which is an equilibrium tide model. This code allows the user to compute the evolution of the orbital distance, eccentricity and inclination of planets as well as their rotation state (obliquity and rotation period). It also computes consistently the rotation period of the host star (taking account the spinup due to radius shrinking and the effects of tides).
The evolution tracks of the radius of various host bodies are implemented: brown dwarfs of mass between 0.01 and 0.08 M_{⊙}, Mdwarfs of 0.1 M_{⊙}, Sunlike stars and Jupiter. This allows the user to study the influence of a changing radius of host body on the tidal evolution of planets.
We endeavored in this work to validate our code. To this purpose, we compared the outputs of a code solving the tidal secular equations of singleplanet systems (see Bolmont et al. 2011, 2012), with the outputs of our new code. We also tested the rotational flattening effect by comparing MercuryT with the CR13 code, that was developed independently. We found that MercuryT reproduces well the secular evolution of the planets. We made sure that the total angular momentum was conserved in all of our examples.
Fig. 12 Shortterm (100 000 yr) evolution of the eccentricity of the Kepler62 system planets. 

Open with DEXTER 
We advise the potential users of this code to keep in mind that, when a planet is alone in the system, the code can produce a spurious remnant eccentricity. We also advise to verify for each simulation the conservation of total angular momentum and the robustness of the spin integration (by doing a simulation without tides and with the effect of the rotationalinduced flattening to see if there is a drift of the mean value of the obliquity. If there is a drift, the time step has to be decreased).
Some ongoing improvements of this code consist in improving the models of the planets. Indeed, the use of the constant time lag model for terrestrial planet has been criticized (e.g., Makarov & Efroimsky 2013; Efroimsky & Makarov 2013; Makarov & Berghea 2014; Correia et al. 2014) and it is probable that the planets are not evolving towards pseudosynchronization but are trapped in spinorbit resonances. Besides, in order to correctly determine the spin of planets, one needs to take into account the thermal tides (e.g. Cunha et al. 2015; Leconte et al. 2015). With a Global Circulation Model, Leconte et al. (2015) showed that this phenomenon can drive planets out of synchronization even if they have a thin atmosphere. We also intend to implement a better description of the dissipation within the star (e.g., using models as in AuclairDesrotour et al. 2014). A wind prescription will also be soon added (like in Bolmont et al. 2012). We also intend to investigate in the future the multibulge effect, i.e. the influence of the bulge raised on the star by planet j on the dynamical evolution of planet i≠j (such as was done in Touma & Wisdom 1994, for the EarthMoonSun system).
MercuryT is a very powerful tool to simulate the evolution of any kind of planetary systems. It can be used to simulate known exoplanetary systems to try to identify some trends, such as was done in this work for the Kepler62 system: investigate the stability of the system taking into account all the important physical phenomena, investigate the influence of tidal dissipation factors on the evolution of the system to maybe constrain the parameters space.
Bolmont et al. (2013) used a previous version of the code to evaluate the possible eccentricity of the transiting inner planet of the 55 Cancri system. Thanks to this code, they investigated if tidal heating could significantly contribute to 55 Cancri e’s thermal emission.
Our code also allows the user to have an idea of the spin state of the planets (as in this work or in Bolmont et al. 2014b, which focuses on the Kepler186 system and also constitutes a fine example of the use of MercuryT).
MercuryT is particularly interesting to use to simulate the orbital dynamical evolution of habitable planets because it allows to have reasonable and consistent inputs for climate models to investigate the potential of these planets to host surface liquid water and also to investigate the influence of eccentricity oscillations on such a climate.
The link to this code and the manual can be found here: http://www.emelinebolmont.com/researchinterests
Acknowledgments
The authors would like to thank John Chambers for his code Mercury and Christophe Cossou for having updated the Mercury code in fortran 90, and for his help to develop this code. E.B. acknowledges that this work is part of the F.R.S.FNRS “ExtraOrDynHa” research project. A.C. acknowledges support from CIDMA strategic project UID/MAT/04106/2013. The authors would like to thank Rosemary Mardling for a thorough referee report that helped them improve the quality of the manuscript.
References
 Alexander, M. E. 1973, Ap&SS, 23, 459 [NASA ADS] [CrossRef] [Google Scholar]
 AuclairDesrotour, P., Le PoncinLafitte, C., & Mathis, S. 2014, A&A, 561, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403 [NASA ADS] [Google Scholar]
 Batygin, K., Bodenheimer, P., & Laughlin, G. 2009, ApJ, 704, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Beaugé, C., & Nesvorný, D. 2012, ApJ, 751, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, A., Loutre, M. F., & Laskar, J. 1992, Science, 255, 560 [NASA ADS] [CrossRef] [Google Scholar]
 Bolmont, E., Raymond, S. N., & Leconte, J. 2011, A&A, 535, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bolmont, E., Raymond, S. N., Leconte, J., & Matt, S. P. 2012, A&A, 544, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bolmont, E., Selsis, F., Raymond, S. N., et al. 2013, A&A, 556, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bolmont, E., Raymond, S. N., & Selsis, F. 2014a, in SF2A2014: Proc. Annual meeting of the French Society of Astronomy and Astrophysics, eds. J. Ballet, F. Martins, F. Bournaud, R. Monier, & C. Reylé, 63 [Google Scholar]
 Bolmont, E., Raymond, S. N., von Paris, P., et al. 2014b, ApJ, 793, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Borucki, W. J., Agol, E., Fressin, F., et al. 2013, Science, 340, 587 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Chabrier, G., & Baraffe, I. 1997, A&A, 327, 1039 [NASA ADS] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Correia, A. C. M., & Robutel, P. 2013, ApJ, 779, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Rodríguez, A. 2013, ApJ, 767, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Udry, S., Mayor, M., et al. 2005, A&A, 440, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dyn. Astron., 111, 105 [NASA ADS] [CrossRef] [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., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Couetdic, J., Laskar, J., Correia, A. C. M., Mayor, M., & Udry, S. 2010, A&A, 519, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cunha, D., Correia, A. C. M., & Laskar, J. 2015, Int. J. Astrobiol., 14, 233 [CrossRef] [Google Scholar]
 Dressing, C. D., Spiegel, D. S., Scharf, C. A., Menou, K., & Raymond, S. N. 2010, ApJ, 721, 1295 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Bonomo, A. S., Haywood, R. D., et al. 2014, ApJ, 789, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Efroimsky, M., & Makarov, V. V. 2013, ApJ, 764, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Eggleton, P. P., Kiseleva, L. G., & Hut, P. 1998, ApJ, 499, 853 [NASA ADS] [CrossRef] [Google Scholar]
 Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298 [NASA ADS] [CrossRef] [Google Scholar]
 Fienga, A., Manche, H., Laskar, J., & Gastineau, M. 2008, A&A, 477,315 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661 [NASA ADS] [CrossRef] [Google Scholar]
 Hairer, E., Norsett, S. P., & Wanner, G. 1993, Solving Ordinary Differential Equations I (Springer) [Google Scholar]
 Hansen, B. M. S. 2010, ApJ, 723, 285 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R., Williams, D., Kipping, D., et al. 2014, Astrobiology, 14, 798 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
 Kidder, L. E. 1995, Phys. Rev. D, 52, 821 [NASA ADS] [CrossRef] [Google Scholar]
 Kokubo, E., & Ida, S. 2007, ApJ, 671, 2082 [NASA ADS] [CrossRef] [Google Scholar]
 Lambeck, K. 1977, Roy. Soc. London Philo. Trans. Ser. A, 287, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1997, A&A, 317, L75 [NASA ADS] [Google Scholar]
 Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., Boué, G., & Correia, A. C. M. 2012, A&A, 538, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leconte, J., & Chabrier, G. 2013, Nature Geoscience, 6, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Leconte, J., Chabrier, G., Baraffe, I., & Levrard, B. 2010, A&A, 516, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leconte, J., Lai, D., & Chabrier, G. 2011, A&A, 528, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013, Nature, 504, 268 [NASA ADS] [CrossRef] [Google Scholar]
 Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Mahajan, N., & Wu, Y. 2014, ApJ, 795, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Makarov, V. V., & Berghea, C. 2014, ApJ, 780, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Makarov, V. V., & Efroimsky, M. 2013, ApJ, 764, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Mardling, R. A. 2007, MNRAS, 382, 1768 [NASA ADS] [Google Scholar]
 Mardling, R. A. 2010, MNRAS, 407, 1048 [NASA ADS] [CrossRef] [Google Scholar]
 Mardling, R. A., & Lin, D. N. C. 2002, ApJ, 573, 829 [NASA ADS] [CrossRef] [Google Scholar]
 Merlis, T. M., & Schneider, T. 2010, J. Adv. Mod. Earth Systems, 2, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Mignard, F. 1979, Moon and Planets, 20, 301 [NASA ADS] [CrossRef] [Google Scholar]
 Milankovitch, M. 1941, Kanon der Erdebestrahlung und seine anwendung auf das eiszeitenproblem (Koniglich Serbische Akademie) [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (Cambridge University Press) [Google Scholar]
 Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2011, Nature, 473, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [NASA ADS] [Google Scholar]
 Pollack, H. N., Hurter, S. J., & Johnson, J. R. 1993, Rev. Geophys., 31, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Quintana, E. V., Barclay, T., Raymond, S. N., et al. 2014, Science, 767, 128 [Google Scholar]
 Spencer, J. R., Jessup, K. L., McGrath, M. A., Ballester, G. E., & Yelle, R. 2000, Science, 288, 1208 [NASA ADS] [CrossRef] [Google Scholar]
 Spiegel, D. S., Menou, K., & Scharf, C. A. 2009, ApJ, 691, 596 [NASA ADS] [CrossRef] [Google Scholar]
 Spiegel, D. S., Burrows, A., & Milsom, J. A. 2010, in AAS/Division for Planetary Sciences Meeting Abstracts, 42, 27.27 [Google Scholar]
 Touma, J., & Wisdom, J. 1994, AJ, 108, 1943 [NASA ADS] [CrossRef] [Google Scholar]
 Touma, J., & Wisdom, J. 1998, AJ, 115, 1653 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y., & Goldreich, P. 2002, ApJ, 564, 1024 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y., & Murray, N. 2003, ApJ, 589, 605 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Values of the reduced dissipation factor for the host bodies implemented in MercuryT.
All Figures
Fig. 1 Twodimensional schematic representing the two deformed bodies. The star is divided into three masses: the central one of mass M_{⋆} − 2μ_{⋆} at S, and two bulges of mass μ_{⋆} at S′ and S′′. Planet j is divided into three masses: the central one of mass M_{pj} − 2μ_{pj} at P_{j}, and two bulges of mass μ_{pj} at and . Ω_{⋆} is the star rotation vector (its norm is Ω_{⋆}, the star rotation frequency), Ω_{pj} is planet j rotation vector (its norm is Ω_{pj}, the planet rotation frequency), and is a vector collinear with the orbital angular momentum of planet j (its norm is equal to the derivative of the true anomaly). e_{rj} is the radial vector. 

Open with DEXTER  
In the text 
Fig. 2 Integration scheme of MercuryT. 

Open with DEXTER  
In the text 
Fig. 3 Case (1): tidal evolution of a planet of mass 1 M_{⊕} orbiting a 0.08 M_{⊙} BD calculated with the secular code (blue dashed line) and the MercuryT code (full red line). Graph a) from top to bottom: evolution of the semimajor axis of the planet, evolution of its eccentricity and evolution of the tidal heat flux. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its rotation period (the pseudosynchronization period is represented in long red dashes), and conservation of total angular momentum. 

Open with DEXTER  
In the text 
Fig. 4 Case (3): tidal evolution of a planet of mass 1 M_{⊕} orbiting a 0.08 M_{⊙} BD calculated with the secular code (blue dashed line) and the MercuryTcode (full red line). Graph a) from top to bottom: evolution of semimajor axis (in red) and of the corotation distance (in black), evolution of eccentricity and evolution of tidal heat flux. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its inclination and evolution of its rotation period (in red) and the BD rotation period (in black). The pseudosynchronization period is represented in long red dashes. 

Open with DEXTER  
In the text 
Fig. 5 Case (3′): tidal evolution of a Jupitermass planet orbiting a 0.08 M_{⊙} BD calculated with the secular code (blue dashed line) and the MercuryT code (full red line). Graph a) from top to bottom: evolution of semimajor axis (in red) and of the corotation distance (in black), evolution of eccentricity and conservation of angular momentum. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its inclination and evolution of its rotation period (in red) and the BD rotation period (in black). The pseudosynchronization period is represented in long red dashes. 

Open with DEXTER  
In the text 
Fig. 6 Case (4): tidal evolution of an Earthmass planet orbiting a 0.08 M_{⊙} BD calculated with the secular code (blue dashed line) and the MercuryT code (full red line). Graph a) from top to bottom: evolution of semimajor axis (in red) and of the corotation distance (in black), evolution of eccentricity and evolution of tidal heat flux. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its rotation period (in red), the BD rotation period (in black), and the pseudosynchronization period (in long red dashes), and evolution of α. 

Open with DEXTER  
In the text 
Fig. 7 Evolution of (P_{p} − P_{p,0}) /P_{p,0} of the inner planet of the system corresponding to cases (6) of Table 5. The colored lines correspond to the results of MercuryT: red with a time step of 0.08 day, green with a time step of 0.05 day, light blue with a time step of 0.01 and dark blue with a time step of 0.001 day. 

Open with DEXTER  
In the text 
Fig. 8 Evolution of the obliquity of the inner planet of the system corresponding to case (6′′) of Table 5 for the last 600 yr of the evolution. The black line corresponds to the results of the CR13 code. The red line corresponds to the results of MercuryT. 

Open with DEXTER  
In the text 
Fig. 9 Tidal evolution of the Kepler62 system (ℬ). Top panel: evolution of the obliquities of the five planets. Bottom panel: evolution of their rotation periods in colored full lines. The dashed lines correspond to the pseudosynchronous rotation period and the full black line corresponds to the rotation period of the star. 

Open with DEXTER  
In the text 
Fig. 10 Shortterm (100 000 yr) evolution of the rotation period of Kepler62b (ℬ). Full line: rotation period, dashed line: pseudosynchronization period, dasheddotted line: synchronization period. 

Open with DEXTER  
In the text 
Fig. 11 Longterm tidal evolution of the two outer planets of the Kepler62 system. Top panel: evolution of the obliquities of the five planets. Bottom panel: evolution of their rotation periods in colored full lines. The dashed lines correspond to the pseudosynchronous rotation period and the full black line corresponds to the corotation radius. 

Open with DEXTER  
In the text 
Fig. 12 Shortterm (100 000 yr) evolution of the eccentricity of the Kepler62 system planets. 

Open with DEXTER  
In the text 