© ESO, 2014

1. Introduction

It has been shown that planets in first-order mean-motion resonances (MMRs) that undergo tidal dissipation, naturally leave the resonant configuration by moving away from each other (Papaloizou & Terquem 2010; Papaloizou 2011; Lithwick & Wu 2012; Delisle et al. 2012; Batygin & Morbidelli 2013b). The tidal dissipation first induces a decrease of both eccentricities (as expected) and the system initially stays in resonance. However, when eccentricities reach low values, the ratio between the orbital periods of the outer planet and the inner one begins to increase (diverging orbits) as the eccentricities continue to decrease. If the timescale of the dissipation is sufficiently short (compared to the age of the system), the period ratio can significantly depart from the resonant value.

It is important to note that during this process, the system never crosses the resonance separatrix. The separatrix simply disappears at low eccentricities, and the system ends up with a period ratio P_out/P_in greater than the resonant value (e.g., Delisle et al. 2012). However, if the amplitude of libration in the resonance becomes sufficiently high, the system may cross the separatrix before it disappears and may end up either in the internal or the external circulation areas (e.g., Novak et al. 2003; Goldreich & Schlichting 2014). External circulation refers to the configuration where planets are close to a MMR (p + q):p, but with a period ratio greater than the resonant value (P_out/P_in> (p + q) /p). On the contrary, internal circulation refers to the configuration P_out/P_in< (p + q) /p.

In this study, we obtain a simple criterion for the dissipation undergone by planetary systems in a MMR of any order to end up in the resonant area, the internal, or the external circulations areas. In Sect. 2 we present a method for obtaining a simple integrable model for resonances of any order. In Sect. 3 we consider the evolution of the dynamics under dissipation. We focus on the evolution of the amplitude of libration in the resonance depending on the balance of dissipation in both planets. In Sect. 4 we apply our reasoning to the GJ 163 (3:1 MMR) and GJ 581 (5:2 MMR) planetary systems. Finally, in Sect. 5 we discuss the impact of resonance breaking induced by tides on the Kepler mission’s multi-planetary systems statistics.

2. An integrable simplified model of resonances

In order to obtain simple criteria on the dissipation we need an analytical model of MMRs that is as simple as possible and that still captures the main characteristics of the resonant dynamics. In particular, the model should be integrable (one degree of freedom). For first-order resonances, this integrable approximation is easily obtained (see Sessin & Ferraz-Mello 1984; Henrard et al. 1986; Wisdom 1986; Batygin & Morbidelli 2013a). However, for higher order resonances, the situation is more complex and we need to make some simplifying hypotheses.

First, we introduce our notations and the classical method for obtaining a non-integrable Hamiltonian model for any MMR. Then, we reproduce the well-known method used to obtain an integrable model for first-order resonances as a template for higher order resonances.

2.1. Hamiltonian of a (p + q):p MMR

We refer to the star as body 0, to the inner planet as body 1, and to the outer planet as body 2. Noting the masses of the three bodies m_i, we introduce for both planets and β_i = m₀m_i/ (m₀ + m_i), where is the gravitational constant. Let r_i be the position vectors of the planets with respect to the star and the canonically conjugated momenta (in astrocentric coordinates, see Laskar & Robutel 1995). As usual in the literature, semi-major axes are noted a_i, eccentricities e_i, mean longitudes λ_i, and longitudes of periastron ϖ_i. It should be noted that we only consider the planar case in this study. The Hamiltonian of the three-body problem reads (1)where is the Keplerian part (star-planets interactions) and is the perturbative part (planet-planet interactions). The Keplerian part is given by (2)where is the circular angular momentum of planet i, (3)The perturbative part can be decomposed in direct and indirect interactions, with Δ₁₂ = ||r₁ − r₂||. This perturbation can be expressed as a function of elliptical orbital elements by expanding it in Fourier series of the angles λ_i,ϖ_i (e.g., Laskar & Robutel 1995).

For a given MMR (p + q):p between both planets, the corresponding combination of the mean longitudes undergoes slow variations (7)whereas other (non-resonant) combinations of these angles circulate rapidly. The long-term evolution of the orbits is accurately described by the averaged Hamiltonian over the non-resonant combinations of the mean longitudes.

By performing this averaging and the classical angular momentum reduction, one obtains a two degrees of freedom problem (e.g., Delisle et al. 2012) and two constants of motion. The first constant of motion is the total angular momentum, (8)The second, which comes from the averaging, is a combination of the circular angular momenta (or semi-major axes, e.g., Michtchenko & Ferraz-Mello 2001), (9)As shown in Delisle et al. (2012), the constant Γ can be used as a scaling factor and does not influence the dynamics of the system except by changing the scales of the problem (in space, energy, and time). The elimination of Γ is achieved by performing the following change of coordinates (see Delisle et al. 2012) while angle coordinates are unchanged. Using these new coordinates, the dynamics of the system depends on only one parameter, G, the renormalized angular momentum.

The remaining two degrees of freedom can be represented by both resonant angles, (15)and both deficits of angular momentum (Laskar 2000), which are canonically conjugated to the resonant angles, (16)We can also introduce rectangular coordinates, (17)It should be noted that for small eccentricities, | x_i | ∝ e_i. The averaged Hamiltonian takes the form (18)where is the secular part of the Hamiltonian depending on the difference of longitudes of periastron (Δϖ = σ₂ − σ₁), but not on mean longitudes of the planets, and ℛ is the resonant part. These two parts can be expanded as power series of eccentricities (more precisely of x_i). The method used to obtain these expansions is presented in Laskar & Robutel (1995); Delisle et al. (2012).

The Keplerian part can be expressed as a function of the momenta I_i by substituting the expressions of Λ_i in Eq. (2), with the total angular momentum deficit (Laskar 2000), (21)The secular part contains terms of degree two and more in eccentricities while the resonant part contains terms of degree q and more. Thus, the simplest model of the resonance should take into account at least those terms of order q in eccentricities in the perturbative part, (22)where is the secular part truncated at degree q and R_k are constant coefficient (see Delisle et al. 2012).

This problem is much simpler than the initial four degrees of freedom problem. However, in general it is still non-integrable since it presents two degrees of freedom.

2.2. First-order resonances

For a first-order MMR (such as the 2:1 or 3:2 resonances) the simplest Hamiltonian reads (23)where there are no secular terms since they only appear at degree two. It is well known that the Hamiltonian (23) is integrable (see Sessin & Ferraz-Mello 1984; Henrard et al. 1986; Wisdom 1986; Batygin & Morbidelli 2013a). Introducing R and φ such that and the new coordinates u₁, u₂ such that (26)with R_φ the rotation matrix (27)the Hamiltonian (23) reads (28)with (29)In these coordinates, the Hamiltonian (28) does not depend on the angle associated with u₂ but only on the action which is thus a new constant of motion of the system. We are then left with only one degree of freedom (u₁, ) and the problem is integrable.

However, if one introduces second-order terms in the Hamiltonian of a first-order MMR, this simplification no longer occurs. It does not occur either, in the case of higher order MMRs, even at the minimal degree of development of the Hamiltonian. Moreover, the presence of chaotic motion proves that no reduction to an integrable system is possible (see Wisdom 1986).

For those higher order MMRs, we have to make some additional hypothesis and/or to neglect some terms in the Hamiltonian in order to obtain an integrable system that still captures most of the characteristics of the resonant motion. Nevertheless, we can get some inspiration from the method described for first-order resonances. The main idea of this method is to perform a rotation of coordinates and introduce some kinds of proper modes of the motion (u_i). This rotation is chosen such that the fixed point corresponding to the center of libration of the resonance lies in the direction of the first mode (i.e., u₁ ≠ 0, u₂ = 0 at the libration center). Thus, all the resonant dynamics is concentrated in the first mode. The second mode only adds some circulation around the libration center.

The aim of our method for higher order resonance is to obtain a similar rotation that concentrates the main characteristics of the resonant motion in one mode. Of course, we will have to make some simplifications and the model will not be able to reproduce the chaotic motion that is observed in the numerical simulations.

2.3. Higher order resonances

We now consider a resonance (p + q):p of order q> 1 and the Hamiltonian (22). The position of the center of libration can be obtained by solving Hamilton’s equations with zero as the right-hand member (e.g., Delisle et al. 2012), (30)Let I_i,ell, σ_i,ell be the coordinates of this libration center. We then introduce the angle φ such that (31)and the diagonal matrix M_σ with (M_σ)_i,i = e^iσ_i,ell. It can easily be checked that the canonical change of coordinates (32)puts the libration center in the direction of u₁ (u₂ = 0). We introduce the action-angle coordinates D_i,θ_i such that . The new Hamiltonian has exactly the same form as the Hamiltonian (22), (33)where we still have . In the case of first-order resonances, the angle θ₂ no longer appears in the Hamiltonian, but here we still have a dependency on θ₂ in the secular and the resonant parts. In order to go further we need a simplifying approximation. We note that by construction, at the libration center, we have D_2,ell = 0 while D_1,ell ≠ 0. Close to this elliptic fixed point, D₂ is negligible compared to D₁. The Hamiltonian (33) can be expanded in a power series of , which is a small parameter as long as the system is close to the center of libration. We approximate this Hamiltonian with the zeroth order part of the expansion that does not depend on u₂, but only on u₁. This is equivalent to imposing u₂ = 0 (D₂ = 0). In this approximation, the Hamiltonian is much simpler and integrable, (34)Moreover, we have and the secular part becomes simply a polynomial of degree ⌊ q/ 2 ⌋ in . We note that imposing D₂ = 0 is a strong hypothesis. This implies that I₂/I₁ ≈ Λ₂/ Λ₁(e₂/e₁)² stays locked at its value at the libration center (I_2,ell/I_1,ell = tan²φ). Therefore, the eccentricities stay very close to their values at the libration center. This hypothesis is reasonable when the system is close to the libration center. When the system is close to the separatrix, the dynamics is more complex, but we keep this simplified model as a first approximation.

2.4. Development of the Keplerian part

Since we obtained an integrable Hamiltonian there is no need for a new simplification. However, if we study the dynamics at low eccentricities, the Keplerian part can be approximated with a polynomial in by a Taylor expansion. This allows one to obtain polynomials equations of motion.

We can first rewrite Eqs. (19) and (20) in the form with ΔG = G₀ − G, and G₀, Λ_i,0 the values at the nominal resonance and zero eccentricities (see Delisle et al. 2012). For a resonant system, the semi-major axis ratio (or the period ratio) should be close to the nominal resonant value. Therefore, we should have Λ_i ≈ Λ_i,0 and thus . ΔG can thus be interpreted as an estimate (or a generalization) of the total angular momentum deficit (Laskar 2000). It should be observed, however, that ΔG is a constant of motion while is not conserved. Since, in the resonant case is small, we develop the Keplerian part in power series of . The constant terms of this development can be discarded because they have no impact on the dynamics. Moreover, it can be shown that, because of the definition of Λ_i,0, there is no contribution at the first order. The development should then be done at least at the second order, (37)where is a constant that will be ignored in the following and (38)where n_1,0 is the Keplerian mean motion at the nominal resonance, (39)We note that the degree of development of the Keplerian part should be consistent with the development of the perturbative part. For a resonance of order q, we developed the perturbative part at degree q, thus the Keplerian part should also be developed at order q in eccentricity (order ⌊ q/ 2 ⌋ in ). For the sake of simplicity we only consider here resonances of order q< 6 for which the Keplerian part can be approximated with (40)where we eliminated constant terms and used the hypothesis (D₂ = 0).

2.5. Final form of the Hamiltonian

The only remaining step is to express explicitly the secular part. Assuming q< 6, has the form (41)where we dropped constant terms, and the Hamiltonian reads (42)We define (43)such that (44)One can also rescale the Hamiltonian (and the timescale of the dynamics) with the transformation The final form we obtain for the Hamiltonian of a MMR of order q is (48)In the following, we will drop the stars for the sake of brevity and will study the Hamiltonian (49)or equivalently (50)For the sake of simplicity we will assume that R> 0. If R< 0, one only needs to perform the rotation to meet this hypothesis.

It is important to note that the change of coordinates x → u may depend on the parameter ΔG since the direction of the libration center can change with ΔG. However, for second-order resonances and at leading order in eccentricities, the direction always stays the same and the change of coordinates is constant. Moreover, for any MMR a constant change of coordinates can always be used as a first approximation at low eccentricities (e.g. Beaugé et al. 2006; Michtchenko et al. 2006, which gives the position of this libration center for most common MMRs and a variety of mass ratios). We thus suppose in the following that the change of coordinates is constant.

2.6. Dynamics in the integrable model

Using Hamilton equations and the Hamiltonian (49) we obtain (51)For q> 1 there is an obvious fixed point at u₁ = 0, whereas for q = 1 (first-order MMR) this fixed point never exists.

For u₁ ≠ 0, the position of the fixed points are given by (52)The first part of this equation is obviously real, thus is also real. This means that θ₁ = kπ/q and . We thus have to find the positive roots of (53)If δ is big enough, we can write D₁ ~ δ. At the first order in R the solution would be (54)Finally, the position of the fixed point in terms of u₁ can be approximated with (55)Figure 1 shows the bifurcation diagram of fixed points for a second-order resonance with R = 0.1 and Fig. 2 shows the three different cases for the phase space depending on the number of fixed points (1, 3, or 5). For δ< − R (Fig. 2, top), the system admits only one elliptical (stable) fixed point at zero (zero eccentricities) and the phase-space exhibits only (external) circulation around this fixed point. This clearly corresponds to a non-resonant (or secular) dynamics. For −R<δ<R (Fig. 2, middle), the fixed point at zero becomes hyperbolic (unstable) and two symmetrical elliptical fixed points bifurcate from it on the real line. Each one is associated with a resonant area. The circulation area still exists, surrounding both resonant areas. The separatrix between both types of motion passes through the hyperbolic point at zero. Finally, for δ>R (Fig. 2, bottom), the central fixed point becomes stable again and two symmetrical hyperbolic fixed points appear on the imaginary line. A new circulation area appears (internal circulation), between both resonant areas, around the stable fixed point at zero.

Fig. 1 Position of the fixed points of a second-order resonance (q = 2) in the simplified model (Hamiltonian (49)) as a function of δ. We give the real (top) and imaginary (bottom) parts of u1 as given by Eq. (55) with R = 0.1. Continuous lines correspond to stable branches while dashed lines to unstable ones.

Fig. 2 The three typical phase spaces of the simplified model (Hamiltonian (49)) for a second-order resonance. We plot the energy levels in the complex plane (real and imaginary parts of u1), for R = 0.1 and δ = −0.2 (top), 0 (middle), 0.2 (bottom). Elliptical (stable) fixed points are marked with dots, while hyperbolic (unstable) ones are marked with crosses. R corresponds to a resonant area and Ce (respectively Ci) to external (respectively internal) circulation.

Higher order MMRs exhibit very similar behavior: successive bifurcations from the fixed point at zero. However, the case of first-order MMRs is different. The fixed point at zero eccentricities does not exist for q = 1, and the bifurcation between purely secular motion and a resonant phase space significantly differs (see Henrard & Lemaitre 1983; Delisle et al. 2012). In particular, at the bifurcation, the center of external circulation becomes the center of libration of the resonance. This means that under dissipation (that induces an evolution of the phase space), the system can evolve from external circulation to resonant motion (and vice versa) without crossing the separatrix of the resonance. This phenomenon was invoked to explain the excess of systems in Kepler data that rely close to first-order MMRs (2:1 and 3:2), but in external circulation (see Delisle et al. 2012). For higher order resonances, the same process does not occur because the size of the resonant areas tends to zero when reaching the bifurcation (see Figs. 1 and 2). Thus, a system always needs to cross the separatrix to pass from resonant motion to circulation.

3. Dissipation in a resonant system

3.1. Amplitude of libration

The aim of this section is to study the evolution of a system that is initially trapped in resonance and undergoes some dissipation, such as the tidal effect, that damps the eccentricities of the planets. At the beginning of the process, the system is assumed to be in the resonant area of the phase space (with δ² ≫ Rδ^{q/ 2}; see Fig. 2, bottom). We will assume here that the mode u₁ (defined in Eq. (32)) is damped on a timescale T_d, (56)and that this dissipative force induces a decrease of the parameter δ, which is in first approximation proportional to the action , (57)We present in Sect. 3.2 estimates of T_d and γ in the case of tidal effect in both planets.

We now introduce a measure of the amplitude of libration in the resonant area, (58)where θ_1,max is the maximum value reached by the resonant angle during a libration. We have A = 0 at the center of libration, A = 1 at the separatrix, and 0 <A< 1 in between. The evolution of the amplitude A is governed by the following proposition.

Proposition 1. Assuming that the conservative evolution of u ₁ is given by the Hamiltonian (49), that the dissipation affects the system as described by Eqs. (56) and (57), and that the system is in the regime δ² ≫ Rδ^{q/ 2}, the amplitude of libration in the resonance ( A , Eq. (58)) follows (demonstration in Appendix A) (59)

Therefore, if γ< 8 / (4 + q) the amplitude of libration decreases with respect to the size of the resonance, while if γ> 8 / (4 + q) the amplitude increases. Let γ_c be this critical value, (60)In the first case (γ<γ_c), the system evolves closer and closer to the fixed point. In the long-term, the system will follow very closely the stable branch while eccentricities tend to zero. On the contrary, if γ>γ_c, then A grows exponentially and, depending on the initial amplitude, the system can cross the separatrix before eccentricities are completely damped.

After having crossed the separatrix, the system is no longer locked in resonance and the semi-major axis ratio (or the period ratio) evolution depends on the relative strength of the dissipation in each planet.

3.2. Outcome of tidal dissipation in resonance

From our previous computations, it is clear that the most important parameter to determine is γ and especially how it compares with γ_c. The value of T_d is also important since it gives the timescale of the dissipative evolution and must be compared to the age of the considered system.

The parameter γ is mainly influenced by the balance of the dissipation in both planets. The more the first planet dissipates (compared to the second one), the greater γ is.

At leading order in eccentricities the tidal effect induces (e.g., Correia et al. 2011) where T_i is the timescale of the tidal dissipation in planet i. For instance, if we take a constant lag time model (Singer 1968; Mignard 1979), this timescale is given by (e.g., Bonfils et al. 2013) (64)where Δt_i, k_2,i, and R_i are the time lag, the second Love number, and the radius of planet i. The evolution of u₁ can be deduced from Eq. (62), (65)The dissipation timescale T_d is thus (66)The impact of the dissipation on the parameter ΔG (at leading order) is given by (see Delisle et al. 2012) (67)from which we deduce (68)with where G can be approximated with G₀. We thus have (72)Since (p + q) /p> 1, γ₁>γ₂. Therefore, γ decreases with τ, from γ₁ (τ = 0, dissipation in planet 1 dominates) to γ₂ (τ = + ∞, dissipation in planet 2 dominates).

More precisely, it can be shown that γ₁ ∈ [2,2(p + q) /p] and γ₂ ∈ [2p/ (p + q),2]. The exact values depend on the masses of both planets. The lower bounds correspond to the case m₁ ≪ m₂ and the upper bounds to m₂ ≪ m₁. We thus always have (see Eq. (60)) (73)whereas γ₂ can either be greater or smaller than γ_c depending on the considered resonance and the masses of the planets. If γ₂>γ_c, then γ>γ_c regardless of the balance of dissipation in both planets (τ). If γ₂<γ_c, then there exists a critical value of τ, (74)corresponding to the critical value γ_c. If τ>τ_c, then γ<γ_c, and the amplitude of oscillation decreases, whereas if τ<τ_c the amplitude increases.

In the case τ>τ_c, the crossing of the separatrix is possible, and the evolution of the semi-major axis ratio, α = a₁/a₂, just after the crossing can be estimated from Eq. (63), (75)We introduce (76)Then, if τ>τ_α, the semi-major axis ratio decreases (diverging orbits), whereas if τ<τ_α the semi-major axis ratio increases (converging orbits).

These two criteria can be rewritten in a more explicit form. We note that φ gives the direction of the libration center of the resonance. More precisely, we have (77)where all quantities must be evaluated at the libration center. Thus, τ_c and τ_α can be approximated with where (80)To sum up:

• If τ ( = T₁/T₂) >τ_c, the amplitude of libration decreases while the eccentricities are being damped. For first-order resonances, the system leaves the resonance with diverging orbits when planets reach low eccentricities (e.g., Delisle et al. 2012). For higher order resonances the system stays in the resonant configuration.

• If τ<τ_c, the amplitude of libration increases and the system may cross the separatrix depending on the initial amplitude. The smaller τ is, the more quickly the amplitude grows. If the initial amplitude is small and the growth is slow, the eccentricities are being damped before the system reaches the separatrix and it stays in resonance. Thus, τ_c is an upper estimate of the maximum value of τ needed to cross the separatrix. The true maximum value depends on the initial eccentricities of the planets and on the initial amplitude of libration. If the system crosses the separatrix, the subsequent evolution is given by the comparison of τ and τ_α:

  • if τ>τ_α, the orbits are converging (P₂/P₁< (p + q) /p),

  • if τ<τ_α, the orbits are diverging (P₂/P₁> (p + q) /p).

It should be noted that both τ_c and τ_α have very simple expressions (Eqs. (78) and (79)) that do not depend on the coefficient R of the resonant part of the Hamiltonian.

The only difficulty in estimating τ_c and τ_α is estimating the eccentricity ratio e₁/e₂ at the libration center. This ratio (or equivalently the angle φ) can be estimated from the development of the Hamiltonian at leading degree, but a better estimate can be obtained by searching the fixed points of the Hamiltonian developed at a higher degree in eccentricities (e.g., Delisle et al. 2012), or using numerical averaging methods (e.g., Beaugé et al. 2006; Michtchenko et al. 2006), or computations of periodic orbits (e.g. Hadjidemetriou 2002; Antoniadou & Voyatzis 2014).

4. Application to observed planetary systems

In this section we present applications of our model to observed planetary systems. It should be noted that the dissipation timescale has a strong dependency on the distance to the star (see Eq. (64)). Therefore, in most planetary systems, the tidal dissipation is dominated by the contribution of the inner planet, while the tidal effect in the outer planet can be neglected. According to our model, in this case the amplitude of libration always increases with time (see Eq. (73)). Therefore, if the initial amplitude is large enough, these systems can leave the resonance by crossing the separatrix before eccentricities are damped. Moreover, after these systems leave the resonance, the dissipation in the inner planet induces an increase of the period ratio (P₂/P₁) since the semi-major axis of the inner planet decreases (see Eq. (63)). We thus conclude that for most resonant planetary systems, the final outcome of the tidal dissipation process should be external circulation if the initial amplitude of libration was large enough, or resonant motion if the initial amplitude was small.

However, a system can end up in internal circulation if the tidal effect is much more efficient in the outer planet than in the inner one (e.g., in the case of a gaseous inner planet and a rocky outer one), and/or if the eccentricity of the outer planet is much larger than that of the inner planet. Systems observed close to MMRs but in internal circulation are thus of particular interest because we can obtain strong constraints on the nature of the planets (rocky or gaseous). We looked in the exoplanet.eu database (Schneider et al. 2011) for systems observed in internal circulation and with estimated masses compatible with a gaseous inner planet and a rocky outer one. Only a few systems correspond to these criteria. We selected GJ 163b, c (3:1 MMR) and GJ 581b, c (5:2 MMR) to illustrate our model.

4.1. Application to GJ 163 b, c (3:1 MMR)

4.1.1. The detected system

The star GJ 163 is an M-dwarf that hosts three planets (see Bonfils et al. 2013 and Table 1 for the orbital parameters). The two inner planets are close to a 3:1 MMR with a period ratio of 2.97 (internal circulation). The innermost planet’s minimum mass (msini) is estimated to be 10.7 M_⊕, while the other one is about 7.3 M_⊕. The radii of these planets have not been estimated so their density and nature (rocky or gaseous) is unknown. However, it seems reasonable to suppose that the inner planet is gaseous while the second one could be either rocky or gaseous (see Bonfils et al. 2013).

4.1.2. Scenario

The main question we want to answer is: Is there a natural explanation for having a system very close to but outside of the 3:1 MMR, with a period ratio P₂/P₁< 3 (2.97)? It is of course possible that the system is close to the 3:1 MMR just by chance and that it was never locked in this resonance. However, it seems more probable that the system was locked in the 3:1 MMR in the past (for instance due to convergent migration in the protoplanetary disk), and afterwards slightly diverged from this resonant ratio to lower values.

We investigate here the possibility that tidal dissipation in planets may have induced this resonant departure as described in Sect. 3. Moreover, we deduce constraints on the GJ 163 system in order for this scenario to be valid. From Sect. 3 it follows that for the system to leave the resonance with P₂/P₁< (p + q) /p, the ratio τ = T₁/T₂ of tidal dissipation timescales in both planets must verify τ_α<τ<τ_c (see Eqs. (79) and (78) for definitions of τ_α and τ_c). For numerical applications, we use a constant lag time model (Singer 1968; Mignard 1979). Criteria on τ transcribe in criteria on the lag time ratio Δt₂/ Δt₁ by using the expression of T_i given in Eq. (64), (81)with (82)and (83)For this application, the eccentricity ratio is computed using the simplest analytical model (degree 2). We obtain e₁/e₂ ≈ 1.1 at the libration center. The criterion τ>τ_α implies Δt₂/κΔt₁ ≳ 300. The criterion τ<τ_c implies Δt₂/κΔt₁ ≲ 1450. Thus, the scenario that we described should be possible if (84)

4.1.3. N-body simulations

We performed numerical simulations of GJ 163, starting in resonance with different lag time ratios for the tidal effect in both planets. We used a constant Δt model for the dissipation (Singer 1968; Mignard 1979) and the ODEX integrator (e.g., Hairer et al. 2010) as described in Bonfils et al. (2013). Gaseous planets typically have dissipation quality factors Q ~ 10³ − 10⁴, while rocky planets have quality factors in the range Q ~ 10 − 100. This corresponds to lag times in the range Δt ~ 10 − 100 s for gaseous planets, and Δt ~ 10³ − 10⁴ s for rocky ones (Q ≈ 1 / (nΔt), where n is the mean motion of the planet). Since the innermost planet is probably gaseous, its lag time should be in the range Δt₁ ~ 10 − 100 s. However, in order to speed up the simulations, we used a higher tidal lag time. As explained by Bonfils et al. (2013), the timescale of the evolution is roughly inversely proportional to Δt. In order to check that this approximation does not dramatically affect our results we performed simulations with Δt₁ = 10⁷, 10⁶, 10⁵, and 10⁴ s (respectively 5, 4, 3, and 2 orders of magnitude higher than the expected value). The system is integrated during respectively 0.1, 1, 10, and 100 Myr which would roughly correspond to 10 Gyr with Δt₁ = 100 s. As a comparison, the age of the system is estimated to be in the range 1 − 10 Gyr (Bonfils et al. 2013).

We chose the initial elliptical elements of the planets such that the system is initially in resonance but not exactly at the libration center. As explained in Sect. 3, when the initial amplitude of libration is small, the system may not cross the separatrix before eccentricities are very small even if the amplitude increases. This initial amplitude is generated by varying the mean anomaly of the inner planet while all other elliptical elements correspond to an ACR (libration center of the resonance). For M₁ = 0°, the system is initially at the libration center while for M₁ = 180°, it starts at the separatrix (see Fig. 3 for a diagram of initial conditions). The initial eccentricities of the planets are set to about 0.16 and 0.11 in order to be close to the ACR, the perihelia are anti-aligned, and M₂ = 0°. We took a₁ = 0.062 such that when eccentricities go to zero the system ends up approximately at its current position (a₁ ≈ 0.0607).

Figure 4 shows the final period ratio of the planets as a function of the lag time ratio with M₁ = 100° (top) and 140° (bottom). In both cases, we superimposed the results obtained with the four different dissipative timescales (Δt₁ = 10⁷, 10⁶, 10⁵, and 10⁴ s). We see in Fig. 4 that the four curves show very similar features and exhibit only small variations between them. We do not observe a particular trend when varying the timescale of the dissipation. Thus, we can assume that taking a dissipation timescale several order of magnitudes higher than realistic values does not affect the results of this study much, while it speeds up the computations.

As described by our model, we observe in Fig. 4 the three possible final states (external/internal circulation or resonant motion) for the system depending on the value of Δt₂/κΔt₁. We plot in Fig. 5 an example of simulation for each of the three final configurations. For each case, we plot the evolution of the period ratio, the eccentricities, and the resonant angles. The limit between external and internal circulation (Fig. 4) is around Δt₂/κΔt₁ ≈ 250 for both initial amplitudes of libration (M₁ = 100 and 140°). Our analytical estimate (300, see Eq. (84)) is fairly close to this numerical computation. The limit between internal circulation and resonant motion depends (as expected) on the initial amplitude and occurs at Δt₂/κΔt₁ ≈ 750 for M₁ = 100°, and 900 for M₁ = 140°. The analytical model value (1450, see Eq. (84)) is thus significantly higher than these numerical results. However, as explained in Sect. 3 and observed in the simulation, this limit strongly depends on the initial amplitude and/or initial eccentricities. In our simple model we do not take into account these parameters and obtain an upper estimate of this limit which corresponds to a system initially in resonance but close to the separatrix (M₁ → 180°).

Fig. 3 Diagram of initial conditions used in numerical simulations of GJ 163 (3:1 MMR). M1 varies between 0° (center of libration) and 180° (separatrix of the resonance) along the purple line.

In order to illustrate this dependency we performed numerical simulations of the system with different initial amplitude of libration (M₁). The final period ratio of the system as a function of the lag time ratio and the initial amplitude is shown in Fig. 6. For this figure we used Δt₁ = 10⁷ s since we showed that this value does not strongly affect our results and it allows us to speed up the computations.

For a small initial amplitude of libration, we observe that the system always ends up in resonance (green zone in Fig. 6). For higher initial amplitudes, internal (in blue) and external (in yellow/orange) circulation are possible. We observe that the limit between external and internal circulation always occurs for Δt₂/κΔt₁ ≈ 250.

Fig. 4 Final period ratio of GJ 163’s inner planets as a function of the dissipation balance between both planets (Δt2/κΔt1). We fix M1 = 100° (left), 140° (right), and compare the results obtained with Δt1 = 107, 106, 105, and 104 s (top) with an integration time of respectively 0.1, 1, 10, and 100 Myr, which roughly corresponds to 10 Gyr for a realistic value of Δt1 (100 s). We also compare these results with integrations taking into account the third planet that has been detected in the system (bottom; see Bonfils et al. 2013). For this last comparison we use Δt1 = 107 s. We observe only small variations when changing the timescale of the dissipation (Δt1) and no systematic trend. The third planet also does not seem to have a strong effect. The vertical black lines mark the range of Δt2/κΔt1 that may conduct to internal circulation according to our analytical model. The green lines highlight the same range obtained with the simulations. The lower bound (solid green line) is at Δt2/κΔt1 = 260 in the case M1 = 100° (left) and 240 for M1 = 140° (right). It does not vary much between the experiments and it is close to the analytical value (300). The upper bound (dashed green line) is at 750 (left) and 900 (right) while the analytical value is 1450. The difference between the numerical results, and between the numerical and analytical values can easily be explained (see Sect. 4.1.3).

The limit between internal circulation and resonance depends on the initial amplitude of libration. At very high initial amplitudes, the limit tends to Δt₂/κΔt₁ ~ 1100. This is fairly close to the analytically predicted value (1450, see Eq. (84)). The remaining difference between analytical and numerical results probably comes from our simplifying hypothesis. Analytical results are obtained by assuming that e₁/e₂ stays locked at its value at the libration center. While this is reasonable when the system is close to the libration center, when it is close to the separatrix the eccentricities can significantly differ from their values at the libration center.

It should be noted that around M₁ = 120°, some systems, which were expected to stay in resonance, end up in internal circulation (blue points in the green area). This is due to the presence of the separatrix of the secular resonance inside the MMR. Resonances of this type were already studied in the case of the 2:1 and 3:2 MMR by Callegari et al. (2004, 2006). As described by these studies, the secular frequency (dominating the motion of Δϖ) falls to zero near the separatrix and the angle Δϖ librates in the opposite direction from one side to the other of the separatrix. Of course, our simplified model cannot predict this kind of phenomenon. However, it still captures the main features of the resonant motion under dissipation and allows a better understanding of the mechanisms that lead to the three different final states in our simulations.

It is important to notice that what we call initial amplitude of libration is the initial condition of our specific simulations. If we had chosen lower initial eccentricities, the amplitude of libration would have less time to grow, and a larger initial amplitude would be necessary to cross the separatrix of the resonance before reaching zero eccentricities. With higher initial eccentricities, the minimal initial amplitude of libration needed to cross the separatrix would be smaller. Moreover, at moderate eccentricities (~0.3), the phase space of the resonance exhibits bifurcations (e.g., Michtchenko et al. 2006) and when the system crosses these bifurcations the amplitude of libration undergoes jumps. This means that an even smaller initial amplitude of libration would be necessary if the system was initially at larger eccentricities.

If the system was initially in a 3:1 resonance, our model implies that Δt₂/ Δt₁ ~ 500κ. In order to come to a conclusion about the nature of both planets, we need to estimate κ = k_2,1/k_2,2(R₁/R₂)⁵. The Love numbers of both planets should be of the same order of magnitude¹. The radius ratio, raised to the fifth power, has a more significant impact. Since the planets’ radii are unknown we estimate them from the masses using the empirical power law obtained by Weiss et al. (2013), (85)Using this estimate for GJ 163b, c, we obtain (86)The upper x-axis of Fig. 6 is scaled using this value of κ. Finally, we conclude that Δt₂/ Δt₁ ~ 1000. Such a ratio is reasonable only if the first planet is gaseous while the second one is telluric (see Bonfils et al. 2013). It should be noted that the minimum mass of the outer planet (msini) is about 7 M_⊕. If it is a rocky planet, its mass is probably close to this value, and its inclination should be close to 90°.

4.2. Application to GJ 581 b, c (5:2 MMR)

The model and mechanisms we describe in this study are very general and are not limited to the GJ 163 planetary system. Our conclusions are valid for a MMR of any order q. However, the criterion we obtain on the lag time ratio depends on the considered MMR, on the masses and radii of the planets, etc. Thus, the final outcome of the dissipation (internal/external circulation or resonant motion) must be determined for each system individually.

Fig. 5 Evolution of the period ratio (left), the eccentricities (middle), and the resonant angles (right) of GJ 163b and c, for M1 initially set to 100° and Δt2/κΔt1 = 100 (top, departure from resonance to external circulation), 400 (middle, departure from resonance to internal circulation), and 1000 (bottom, staying in resonance). In the middle and right columns, red dots correspond to planet b and green dots to planet c. Time is given in Myr, but Δt1 = 105 s in all these simulations. For a more realistic dissipation (Δt1 = 100 s), the time should be read as Gyr.

In order to illustrate the generality of our mechanism, we performed a similar analysis on the system GJ 581. This system has raised much discussion about the number of detected planets. Most studies agree on the presence of four planets around this M-dwarf and both planets b and c, which we are interested in, are uncontested. We reproduced the orbital parameters of the four-planet system given in Forveille et al. (2011) in Table 1. Planets b and c (see Bonfils et al. 2005; Udry et al. 2007; Mayor et al. 2009; Forveille et al. 2011) have a period ratio of 2.4 which is close to a 5:2 MMR (internal circulation). The inner planet (b) has a minimum mass of about 15.9 M_⊕ while planet c has a minimum mass of 5.3 M_⊕. Applying our analytical criterion to this system we deduce that internal circulation is possible for Δt₂/κΔt₁ ∈ [3,52] (with e₁/e₂ ≈ 0.25). As we did for GJ 163, we performed numerical simulations of the system for different lag time ratios (Δt₂/κΔt₁) and different initial amplitude of libration (M₁). We set the initial conditions to be close to the center of libration of the resonance. The eccentricities of the planets are initially set to about 0.05 and 0.2, the perihelia are anti-aligned, and M₂ = −36°. We took a₁ = 0.042 such that when eccentricities are damped the system ends up approximately at its current position (a₁ ≈ 0.041). For M₁ = 0°, the system is initially at the center of libration, while for M₁ = 90° the system begins on the separatrix. The system is integrated during 5 × 10⁴ yr with Δt₁ = 10⁷ s. The age of the system is estimated to be about 8 Gyr, thus our simulations are approximately equivalent to Δt₁ ~ 10−100 s on the age of the system (which corresponds to a gaseous planet).

The final outcome of these simulations are represented in Fig. 7. We observe a good agreement between analytical estimates and numerical results. In both cases, the current configuration of GJ 581b and c is obtained if Δt₂/ Δt₁ ~ 20κ. It is interesting to note that this lag time ratio is significantly (one order of magnitude) smaller than what we obtained for GJ 163. The main reason for this is the difference in the position of the libration center in terms of eccentricity. For GJ 163 we have e₁/e₂ ≈ 1.1, while here we have e₁/e₂ ≈ 0.25 at the libration center. Since the lag time ratio is proportional to the eccentricity ratio squared, this explains the difference between the results. It should be observed that for most resonances, at low eccentricities, e₁/e₂ increases with m₂/m₁ (e.g., Michtchenko et al. 2006).

Fig. 6 Final period ratio of GJ 163’s inner planets (3:1 MMR) as a function of the dissipation balance between both planets (Δt2/κΔt1, bottom x-axis) and initial amplitude of libration (M1). The coefficient κ = k2,1/k2,2(R1/R2)5 is unknown, but can be estimated using a mass-radius power law. The top x-axis (Δt2/ Δt1) is scaled using the value κ ≈ 2.2 obtained in Eq. (86). When M1 is set to 0°, the system begins at the libration center, whereas when M1 is set to 180° it starts at the separatrix. We fixed Δt1 = 107 s in order to speed up the simulations and integrated the system during 0.1 Myr which roughly corresponds to 10 Gyr for a realistic value of Δt1 (100 s). The horizontal gray lines correspond to both sets of initial conditions used in Fig. 4. The vertical white lines mark the range of Δt2/κΔt1 that may conduct to internal circulation according to our analytical model.

Fig. 7 Same as Fig. 6 but for GJ 581b, and c (5:2 MMR). The top x-axis (Δt2/ Δt1) is scaled using the value κ ≈ 15 obtained in Eq. (87). We note that M1 varies only between 0° (center of libration) and 90° (separatrix), because the resonant combination is 5M2 − 2M1 (factor 2 in front of M1).

As in the case of GJ 163, the radii and Love numbers of both planets are unknown. We suppose equal Love numbers for both planets and estimate the radii using the same empirical power law as for GJ 163 (Weiss et al. 2013), (87)The top x-axis of Fig. 7 is scaled using this value of κ. The lag time ratio of this system should be Δt₂/ Δt₁ ~ 300. As in the case of GJ 163, we conclude that if GJ 581b and c formed by tidal dissipation in the 5:2 MMR, the inner planet should be gaseous and the outer one rocky. The minimum masses (respectively 15.86 M_⊕ and 5.34 M_⊕ for the inner and the outer planets) are compatible with this conclusion.

5. Kepler’s statistics

Multi-planetary systems detected by the Kepler mission present an excess of planet pairs close to first-order MMRs (2:1 and 3:2), but in external circulation (see Lissauer et al. 2011; Fabrycky et al. 2012). Since tidal dissipation in planets involved in a first-order MMR can induce a departure from the resonance to external circulation when eccentricities reach very low values, this scenario has been proposed to explain Kepler’s statistics (Lithwick & Wu 2012; Delisle et al. 2012; Batygin & Morbidelli 2013b). However, the timescale of this dissipation might be too long compared to the age of the systems to explain their present configurations (see Lee et al. 2013; and also Rein 2012). It should be noted that the scenario considered in these studies assumes that the planets leave the resonance with very low eccentricities and without crossing the separatrix of the resonance. The departure of the period ratio from the resonant value is very slow since the tidal dissipation decreases with eccentricities (see Eqs. (61)–(63)).

Nevertheless, a resonant system that initially has a large enough amplitude of libration can cross the separatrix of the resonance while the eccentricities of the planets are still high. The departure of the period ratio from the resonant value is much faster in this case since the eccentricities are higher. As we observed in the introduction of Sect. 4, in most planetary systems the tidal dissipation mainly occurs in the inner planet and the dissipation in the outer one can be neglected. In this case, the amplitude of libration increases with time, and when a system crosses the resonance separatrix, the period ratio increases. Therefore, this mechanism also produces the excess of planets in external circulation observed in Kepler data, but on a shorter timescale.

Neglecting the dissipation in the outer planet, the evolution of the period ratio () after the system crossed the separatrix is given by (88)For the sake of simplicity we will neglect the secular interactions between the planets and suppose that e₁ undergoes an exponential decrease e₁ = e_1,0e^{− t/T₁}, with e_1,0 the eccentricity of the inner planet when the system crosses the separatrix, and T₁ the dissipation timescale (see Eq. (64)). In this approximation, the temporal evolution of the period ratio follows (89)with .

Using Eq. (89), one can obtain the maximum value reached by the period ratio, (90)Figure 8 shows this maximum value as a function of the eccentricity e_1,0. It should be noted that when eccentricities reach low values the mechanism invoked previously (e.g., Delisle et al. 2012) dominates and the period ratio continues to increase but much more slowly.

Fig. 8 Maximum value reached by the period ratio (when t → ∞, and according to Eq. (90)) as a function of the inner planet’s eccentricity when the system leaves the resonance (e1,0). is the period ratio at the exact resonance.

Lee et al. (2013) estimated that tidal dissipation would need at least t ≳ 50T₁ for the system to reach , both for the 2:1 and the 3:2 resonances. Depending on the considered resonance and on the masses of both planets, this can even reach t ≳ 1000 T₁.

Fig. 9 Time required to produce a departure of P2/P1 of 0.03 from the resonant value as a function of the eccentricity of the inner planet when the system leaves the resonance (e1,0). This value is computed using Eq. (89). It has to be compared with the estimate by Lee et al. (2013) considering a scenario of departure from the resonance at very low eccentricities (t/T1 ≳ 50 − 1000).

Using Eq. (89) we can estimate the time needed for the system to reach the same configuration, but supposing it left the resonance with significant eccentricities. We plot in Fig. 9 this estimate (t/T₁) as a function of the eccentricity e_1,0. We see in Fig. 9 a vertical asymptote at e_1,0 ≈ 0.1, because for e_1,0 ≲ 0.1, the limit is smaller than 0.03. In this case, the system can eventually reach 0.03, on a timescale shorter but comparable to the one estimated by Lee et al. (2013) (50–1000 T₁). On the contrary, when e_1,0 ≳ 0.15, the system reaches the desired period ratio on a much shorter timescale: t/T₁ ~ 0.01 − 0.1. This corresponds to a gain of 3–5 orders of magnitude with this alternate scenario.

Lee et al. (2013) discarded many near-resonant systems because tidal dissipation seemed too slow to explain their current configuration. However, this study only considered a scenario of resonant departure at low eccentricities. Some of these discarded systems might actually have formed by crossing the resonance separatrix with non-negligible eccentricities (e_1,0 ≳ 0.15) due to the increase of the amplitude of libration induced by the tidal dissipation. In that case, the evolution of the period ratio after the resonance breaking is about 3–5 orders of magnitude more rapid, and the current configuration can be obtained on a more reasonable timescale.

6. Conclusion

We presented an integrable model of two-planet mean-motion resonances of any order. This model is highly simplified and cannot reproduce all the features of the resonant dynamics. However, it allows us to deduce a very simple criterion on the tidal dissipation undergone by both planets to end up inside the resonance, or on one side or the other of the resonance. The main factors that enter into account are the balance of tidal dissipation between both planets (T₁/T₂ or Δt₂/ Δt₁) and the position of the libration center (especially the ratio e₁/e₂).

Using this criterion on the two inner planets orbiting GJ 163 we deduce that the current period ratio (2.97) could be the outcome of dissipation in the 3:1 MMR provided that Δt₂/ Δt₁ ~ 1000. Using N-body simulations with dissipation we reach the same conclusion with slightly refined bounds for Δt₂/ Δt₁. Both methods clearly imply that the inner planet should be gaseous and the outer planet should be rocky. The minimum masses of both planets (respectively 10.7 M_⊕ and 7.3 M_⊕) are compatible with this hypothesis, but since the inclinations and the radii are currently unknown, some uncertainty remains. We also applied this model to GJ 581b, c and could reproduce the current configuration with tidal dissipation in the 5:2 MMR if Δt₂/ Δt₁ ~ 300. As in the case of GJ 163, we conclude that the inner planet should be gaseous and the outer planet should be rocky, which is compatible with the minimum masses of both planets (respectively 15.86 M_⊕ and 5.34 M_⊕ for the inner and the outer planets).

As we noted in the case of GJ 163, some secondary resonances can affect the outcome of the considered system. Our integrable model of resonances is not able to predict such complex behavior, nor chaotic motion. This might be a limitation for high-order resonances, which may show large chaotic areas. Moreover, we make all our estimates using a constant eccentricity ratio (e₁/e₂) which is computed at the center of libration of the resonance. As eccentricities are being damped, the position of the libration center evolves and the eccentricity ratio is not constant. Depending on the resonance and on the considered range of eccentricities, the changes on the eccentricity ratio at the libration center can be non-negligible (e.g., Michtchenko et al. 2006). Besides, when the amplitude of libration is small, eccentricities of both planets should be close to the values at the libration center, but when the system reaches the separatrix and leaves the resonance, the eccentricities can be significantly different. Therefore, this estimate of the eccentricity ratio is the main limit in our model and in the computation of criteria on the lag time ratio. However, as we observed in the cases of GJ 163 and GJ 581, with this approximation we still obtain a good estimate of the order of magnitude of the lag time ratio and a better understanding of the mechanisms that are at stake in determining the outcome of the dissipative process.

The most interesting cases to study are those in internal circulation because we can obtain strong constraints on the nature of the planets for our scenario to be possible (as we did for GJ 163 and GJ 581). However, our mechanism also applies to many systems that are observed in external circulation. In most cases, the tidal dissipation in the outer planet is negligible compared to the dissipation in the inner planet, and the most probable outcome for the system is external circulation. For first-order MMRs, external circulation can also be obtained when eccentricities are very low without crossing the separatrix of the resonance (the separatrix simply disappears at low eccentricities, e.g., Delisle et al. 2012). However, the subsequent evolution of the period ratio is very slow because of the low eccentricities. Lee et al. (2013) showed that for many systems the evolution of the period ratio is too slow to reach the current value on a reasonable timescale. We show that a gain of 3–5 orders of magnitude on the timescale is obtained by considering a scenario of resonance breaking due to tides at non-negligible eccentricities (e₁ ≳ 0.15). This allows us to explain the presence of an excess of planets in external circulation in Kepler data without introducing any other mechanism than tidal dissipation.

Acknowledgments

We thank the anonymous referee for pointing out the importance of planets’ radii in the estimate of the tidal dissipation and for other constructive comments that improved the quality of this article. This work has been supported by PNP-CNRS, CS of Paris Observatory, PICS05998 France-Portugal program, and FCT-Portugal (PEst-C/CTM/LA0025/2011).

References

Appendix A: Evolution of the amplitude of libration

The aim of this appendix is to prove Prop. 3.1. Assuming δ² ≫ Rδ^{q/ 2}, the Hamiltonian (50) can be approximated with the pendulum-like Hamiltonian (A.1)Noting D₁ = δ + ϵ, we obtain (A.2)Hamilton’s equations give (in the conservative case) The elliptical fixed point corresponds to qθ_1,ell = 0, ϵ_ell = 0, ℋ_ell = δ² + 2Rδ^{q/ 2}. The hyperbolic fixed point is at qθ_1,hyp = π, ϵ_hyp = 0, ℋ_hyp = δ² − 2Rδ^{q/ 2}. Noting Δℋ = ℋ_ell − ℋ_hyp = 4Rδ^{q/ 2}, one can verify that (A.5)For a given conservative resonant trajectory, the maximum value of θ₁ is reached when . This corresponds to (A.6)This clearly provides a measure of the amplitude of libration which we call A, (A.7)We now consider the evolution of this amplitude under dissipation. The dissipation affects the system as described by Eqs. (56) and (57). From these expressions we derive The instantaneous derivative of the amplitude A is given by (see Eq. (A.7)) We average this instantaneous derivative over one libration period in order to evaluate the long-term evolution of A. We need to compute the mean values of ϵ, ϵ², and ϵ³. Because of the symmetry of the problem, ϵ and ϵ³ average out to zero (odd powers of ϵ). We thus simply have (A.15)with (A.16)and Both integrals can be expressed using elliptic integrals and we obtain where K and E are the complete elliptic integrals of the first and second kinds. Finally, at leading order, the evolution of the amplitude A is governed by (A.21)

All Tables

All Figures

	Fig. 1 Position of the fixed points of a second-order resonance (q = 2) in the simplified model (Hamiltonian (49)) as a function of δ. We give the real (top) and imaginary (bottom) parts of u1 as given by Eq. (55) with R = 0.1. Continuous lines correspond to stable branches while dashed lines to unstable ones.
In the text

Fig. 2 The three typical phase spaces of the simplified model (Hamiltonian (49)) for a second-order resonance. We plot the energy levels in the complex plane (real and imaginary parts of u1), for R = 0.1 and δ = −0.2 (top), 0 (middle), 0.2 (bottom). Elliptical (stable) fixed points are marked with dots, while hyperbolic (unstable) ones are marked with crosses. R corresponds to a resonant area and Ce (respectively Ci) to external (respectively internal) circulation.

In the text

	Fig. 3 Diagram of initial conditions used in numerical simulations of GJ 163 (3:1 MMR). M1 varies between 0° (center of libration) and 180° (separatrix of the resonance) along the purple line.
In the text

Fig. 4 Final period ratio of GJ 163’s inner planets as a function of the dissipation balance between both planets (Δt2/κΔt1). We fix M1 = 100° (left), 140° (right), and compare the results obtained with Δt1 = 107, 106, 105, and 104 s (top) with an integration time of respectively 0.1, 1, 10, and 100 Myr, which roughly corresponds to 10 Gyr for a realistic value of Δt1 (100 s). We also compare these results with integrations taking into account the third planet that has been detected in the system (bottom; see Bonfils et al. 2013). For this last comparison we use Δt1 = 107 s. We observe only small variations when changing the timescale of the dissipation (Δt1) and no systematic trend. The third planet also does not seem to have a strong effect. The vertical black lines mark the range of Δt2/κΔt1 that may conduct to internal circulation according to our analytical model. The green lines highlight the same range obtained with the simulations. The lower bound (solid green line) is at Δt2/κΔt1 = 260 in the case M1 = 100° (left) and 240 for M1 = 140° (right). It does not vary much between the experiments and it is close to the analytical value (300). The upper bound (dashed green line) is at 750 (left) and 900 (right) while the analytical value is 1450. The difference between the numerical results, and between the numerical and analytical values can easily be explained (see Sect. 4.1.3).

In the text

Fig. 5 Evolution of the period ratio (left), the eccentricities (middle), and the resonant angles (right) of GJ 163b and c, for M1 initially set to 100° and Δt2/κΔt1 = 100 (top, departure from resonance to external circulation), 400 (middle, departure from resonance to internal circulation), and 1000 (bottom, staying in resonance). In the middle and right columns, red dots correspond to planet b and green dots to planet c. Time is given in Myr, but Δt1 = 105 s in all these simulations. For a more realistic dissipation (Δt1 = 100 s), the time should be read as Gyr.

In the text

Fig. 6 Final period ratio of GJ 163’s inner planets (3:1 MMR) as a function of the dissipation balance between both planets (Δt2/κΔt1, bottom x-axis) and initial amplitude of libration (M1). The coefficient κ = k2,1/k2,2(R1/R2)5 is unknown, but can be estimated using a mass-radius power law. The top x-axis (Δt2/ Δt1) is scaled using the value κ ≈ 2.2 obtained in Eq. (86). When M1 is set to 0°, the system begins at the libration center, whereas when M1 is set to 180° it starts at the separatrix. We fixed Δt1 = 107 s in order to speed up the simulations and integrated the system during 0.1 Myr which roughly corresponds to 10 Gyr for a realistic value of Δt1 (100 s). The horizontal gray lines correspond to both sets of initial conditions used in Fig. 4. The vertical white lines mark the range of Δt2/κΔt1 that may conduct to internal circulation according to our analytical model.

In the text

	Fig. 7 Same as Fig. 6 but for GJ 581b, and c (5:2 MMR). The top x-axis (Δt2/ Δt1) is scaled using the value κ ≈ 15 obtained in Eq. (87). We note that M1 varies only between 0° (center of libration) and 90° (separatrix), because the resonant combination is 5M2 − 2M1 (factor 2 in front of M1).
In the text

	Fig. 8 Maximum value reached by the period ratio (when t → ∞, and according to Eq. (90)) as a function of the inner planet’s eccentricity when the system leaves the resonance (e1,0). is the period ratio at the exact resonance.
In the text

	Fig. 9 Time required to produce a departure of P2/P1 of 0.03 from the resonant value as a function of the eccentricity of the inner planet when the system leaves the resonance (e1,0). This value is computed using Eq. (89). It has to be compared with the estimate by Lee et al. (2013) considering a scenario of departure from the resonance at very low eccentricities (t/T1 ≳ 50 − 1000).
In the text