© ESO, 2015

1. Introduction

More than 1400 exoplanets have now been detected and about 20% of them are part of multi-planet systems¹. Many of these systems are compact and host close-in 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 multi-planet systems, by their effect on the eccentricities but also by their effect on precession rates.

We present here a new code, Mercury-T², which is based on the N-body code Mercury (Chambers 1999). It allows us to calculate the evolution of semi-major 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 non-evolving object (provided we know its mass, radius, dissipation factor and rotation period), but also evolving brown dwarfs (BDs), an evolving M-dwarf of 0.1 M_⊙, an evolving Sun-like star, and an evolving Jupiter.

The dynamics of multi-planet 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 semi-averaged (Mardling & Lin 2002) or non-averaged 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 Kepler-62 system (Borucki et al. 2013).

2. Model description

The major addition to Mercury provided in Mercury-T is the addition of the tidal forces and torques. But we also added the effect of general relativity and rotation-induced 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 N-body 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_jS′′ is the vector , defined in Fig. 1.

Let us define F_tr (for tides radial), P_to,⋆ and P_to,pj (for tides ortho-radial): (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 re-written 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 post-Newtonian 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₂, 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 Mercury-T 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, rotation-induced 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 Runge-Kutta of order 5 within the routine. The Runge-Kutta 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 Runge-Kutta routine. The integration scheme is illustrated in Fig. 2.

If the host body evolves, the moment of inertia of the star, given by: I_⋆ = M_⋆(rg_⋆R_⋆)², 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 close-in 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 rotation-induced 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

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 N-body integrator needs parameters such as: the masses of the planets M_pj, the mass of the host body M_⋆, the semi-major 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

For Earth-like planets and Super-Earths, 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 super-Earth of 10 M_⊕ would have a radius of 1.8 R_⊕.

For the Jupiter-like 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 Mercury-T 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 M-dwarf (dM) of mass 0.1 M_⊙ and a Sun-like 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 Mercury-T 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 M-dwarf, and we use Hansen (2010)’s value of stellar dissipation for the dissipation of the Sun-like 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 free-fall 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.

3. Code verification

In order to validate the tidal part of the code, we first simulated the tidal evolution of one Earth-mass 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 semi-major 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 Mercury-T 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 Mercury-T code. Therefore, the CR13 code will be used for cross-checking some of our results.

We considered a 2-planet 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. Non-evolving 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.

Figure 3 shows the evolution of the semi-major 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 semi-major 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⁷ yr. The decrease of the eccentricity of the planet is accompanied by a decrease of the semi-major axis. The evolution of these two quantities show a good agreement between the secular code and the Mercury-T 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 pseudo-synchronization period, which is here ~48.5 h. The evolution of obliquity and rotation period show a good agreement between the secular code and Mercury-T.

After 2 × 10⁷ yr of evolution, the eccentricity obtained with Mercury-T 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 multi-planet systems. In the examples we will give later, the eccentricity due to planet-planet interactions is typically bigger than 10^-7.

This residual eccentricity is responsible for a non zero averaged tidal heat flux ≲ 10^-2 W/m². However, the instantaneous tidal heat flux reaches values as low as a few 10^-9 W/m². This low value illustrates the fact the “real” eccentricity of the planet must be much lower than what Mercury-T is calculating.

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 i-component 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 Mercury-T 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⁵ 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 Jupiter-mass 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³ 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.2. Non-evolving 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⁸ 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. Non-evolving BD, effect of both tides

This case corresponds to case (3) of Table 4. The planet is initially outside the corotation radius.

Figure 4 shows the evolution of this system. We find that the results of Mercury-T 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² for Earth (Pollack et al. 1993) and between 2.4 and 4.8 W/m² 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 Mercury-T 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 BD-tide 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 pseudo-synchronization 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 Mercury-T). 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 close-in 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.

Figure 6 shows the evolution of this system. The evolution calculated with the secular code is in good agreement with the evolution calculated with Mercury-T. 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 semi-major 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 Mercury-T 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 Mercury-T 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 Mercury-T 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 multi-planet 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¹ yr.

Thus, in order to have an accurate integration with Mercury-T, we need to perform a integration Runge-Kutta 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.

The semi-major 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.

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 Mercury-T 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⁵ 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 non-dissipative force, our code conserves total energy and angular momentum.

Apart from very small differences, Mercury-T and the CR13 code give the same results. We consider this agreement good enough for the study of exoplanets.

4. The case of Kepler-62

Mercury-T 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 Kepler-62, 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 Kepler-62 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 semi-major 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, rotation-flattening and general relativity on the dynamics of the system, we tested the system for five different cases, which are listed in Table 8.

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 Mercury-T. 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).

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. Semi-major 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 rotation-flattening – 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 rotation-flattening 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 pseudo-synchronization. The evolution timescales of these two quantities are shorter than the timescales of evolution of semi-major axis and eccentricity. Figure 9 shows that for Kepler-62 the rotation period of the three inner planets of the system evolves towards pseudo-synchronization 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 Kepler-62 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 multi-planet system, the eccentricity of a planet is excited due to planet-planet 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’ pseudo-synchronization periods also oscillate. In reality, the rotation periods of the planets are not exactly equal to the corresponding pseudo-synchronization period. Figure 10 shows the evolution of the rotation period of Kepler-62b compared to the pseudo-synchronization period and the synchronization period. The pseudo-synchronization period oscillates too fast for the rotation period to be able to follow. The instantaneous rotation period of Kepler-62b oscillates out of phase with the pseudo-synchronization period and with a lower amplitude.

During the 30 million yr of the simulation, the obliquities and rotation periods of the HZ planets Kepler-62e and f do not evolve significantly. So we performed longer simulations for these two outer planets. Assuming an Earth-like dissipation for the two planets, we found that Kepler-62e is likely today to have reached pseudo-synchronization and have low obliquity. Figure 11 shows that after 3 Gyr of evolution, the obliquity has been damped and the rotation pseudo-synchronized. For Kepler-62f the timescales of evolution are higher and Fig. 11 shows that the rotation period is still evolving towards pseudo-synchronization 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 Earth-like 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.

4.3. Consequence of dynamics on the potential habitability of Kepler-62 e and f

Kepler-62e and Kepler-62f 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 semi-major axis and eccentricity, and physical parameters such as the obliquity and the rotation period of the planet (Milankovitch 1941).

Thanks to the Mercury-T code, we can simulate the dynamical evolution of habitable planets within their system, taking into account the rich dynamics occurring in a multi-planet 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 Kepler-62e 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 Kepler-62e 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 Kepler-62e has a pseudo-synchronous rotation (or a near pseudo-synchronous 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 Kepler-62e 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 Kepler-62f 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 Kepler-62f 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 Kepler-62f would have a very different type of climate than its neighbor. A non-negligible 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 multi-planet 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 spin-up 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_⊙, M-dwarfs of 0.1 M_⊙, Sun-like 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 single-planet systems (see Bolmont et al. 2011, 2012), with the outputs of our new code. We also tested the rotational flattening effect by comparing Mercury-T with the CR13 code, that was developed independently. We found that Mercury-T reproduces well the secular evolution of the planets. We made sure that the total angular momentum was conserved in all of our examples.

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 rotational-induced 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 pseudo-synchronization but are trapped in spin-orbit 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 Auclair-Desrotour 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 multi-bulge 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 Earth-Moon-Sun system).

Mercury-T 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 Kepler-62 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 Kepler-186 system and also constitutes a fine example of the use of Mercury-T).

Mercury-T 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.

---

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

All Tables

All Figures