Issue 
A&A
Volume 566, June 2014



Article Number  A137  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201423676  
Published online  25 June 2014 
Resonance breaking due to dissipation in planar planetary systems
^{1} ASD, IMCCECNRS UMR 8028, Observatoire de Paris, UPMC, 77 Av. DenfertRochereau, 75014 Paris, France
email: delisle@imcce.fr
^{2} Departamento de Física, I3N, Universidade de Aveiro, Campus Universitário de Santiago, 3810193 Aveiro, Portugal
Received: 19 February 2014
Accepted: 25 April 2014
We study the evolution of two planets around a star, in meanmotion resonance and undergoing tidal effects. We derive an integrable analytical model of meanmotion resonances of any order which reproduce the main features of the resonant dynamics. Using this simplified model, we obtain a criterion showing that, depending on the balance of the tidal dissipation in both planets, their final period ratio may stay at the resonant value, increase above, or decrease below the resonant value. Applying this criterion to 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 the innermost planet is gaseous (slow dissipation) while the second one is rocky (faster dissipation). We perform Nbody simulations with tidal dissipation to confirm the results of our analytical model. We also apply our criterion on GJ 581b, c (5:2 MMR) and reproduce the current period ratio (2.4) if the inner planet is gaseous and the outer is rocky (as in the case of GJ 163). Finally, we apply our model to the Kepler mission’s statistics. We show that the excess of planets pairs close to firstorder MMRs but in external circulation, i.e., with period ratios P_{out}/P_{in}> (p + 1) /p for the resonance (p + 1):p, can be reproduced by tidal dissipation in the inner planet. There is no need for any other dissipative mechanism, provided that these systems left the resonance with nonnegligible eccentricities.
Key words: celestial mechanics / planets and satellites: general
© ESO, 2014
1. Introduction
It has been shown that planets in firstorder meanmotion 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 multiplanetary 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 firstorder resonances, this integrable approximation is easily obtained (see Sessin & FerrazMello 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 nonintegrable Hamiltonian model for any MMR. Then, we reproduce the wellknown method used to obtain an integrable model for firstorder 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_{0}m_{i}/ (m_{0} + 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, semimajor 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 threebody problem reads (1)where is the Keplerian part (starplanets interactions) and is the perturbative part (planetplanet 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 Δ_{12} = r_{1} − r_{2}. 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 (nonresonant) combinations of these angles circulate rapidly. The longterm evolution of the orbits is accurately described by the averaged Hamiltonian over the nonresonant 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 semimajor axes, e.g., Michtchenko & FerrazMello 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 (Δϖ = σ_{2} − σ_{1}), 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 nonintegrable since it presents two degrees of freedom.
2.2. Firstorder resonances
For a firstorder 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 & FerrazMello 1984; Henrard et al. 1986; Wisdom 1986; Batygin & Morbidelli 2013a). Introducing R and φ such that and the new coordinates u_{1}, u_{2} 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_{2} 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_{1}, ) and the problem is integrable.
However, if one introduces secondorder terms in the Hamiltonian of a firstorder 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 firstorder 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_{1} ≠ 0, u_{2} = 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 righthand 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_{1} (u_{2} = 0). We introduce the actionangle 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 firstorder resonances, the angle θ_{2} no longer appears in the Hamiltonian, but here we still have a dependency on θ_{2} 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_{2} is negligible compared to D_{1}. 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_{2}, but only on u_{1}. This is equivalent to imposing u_{2} = 0 (D_{2} = 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_{2} = 0 is a strong hypothesis. This implies that I_{2}/I_{1} ≈ Λ_{2}/ Λ_{1}(e_{2}/e_{1})^{2} stays locked at its value at the libration center (I_{2,ell}/I_{1,ell} = tan^{2}φ). 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_{0} − G, and G_{0}, Λ_{i,0} the values at the nominal resonance and zero eccentricities (see Delisle et al. 2012). For a resonant system, the semimajor 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_{2} = 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 secondorder 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_{1} = 0, whereas for q = 1 (firstorder MMR) this fixed point never exists.
For u_{1} ≠ 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 θ_{1} = kπ/q and . We thus have to find the positive roots of (53)If δ is big enough, we can write D_{1} ~ δ. At the first order in R the solution would be (54)Finally, the position of the fixed point in terms of u_{1} can be approximated with (55)Figure 1 shows the bifurcation diagram of fixed points for a secondorder 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 phasespace exhibits only (external) circulation around this fixed point. This clearly corresponds to a nonresonant (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 secondorder resonance (q = 2) in the simplified model (Hamiltonian (49)) as a function of δ. We give the real (top) and imaginary (bottom) parts of u_{1} 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 secondorder resonance. We plot the energy levels in the complex plane (real and imaginary parts of u_{1}), 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 C_{e} (respectively C_{i}) to external (respectively internal) circulation. 
Higher order MMRs exhibit very similar behavior: successive bifurcations from the fixed point at zero. However, the case of firstorder 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 firstorder 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 δ^{2} ≫ Rδ^{q/ 2}; see Fig. 2, bottom). We will assume here that the mode u_{1} (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 _{1} 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 δ^{2} ≫ 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 longterm, 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 semimajor 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_{1} 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_{0}. We thus have (72)Since (p + q) /p> 1, γ_{1}>γ_{2}. Therefore, γ decreases with τ, from γ_{1} (τ = 0, dissipation in planet 1 dominates) to γ_{2} (τ = + ∞, dissipation in planet 2 dominates).
More precisely, it can be shown that γ_{1} ∈ [2,2(p + q) /p] and γ_{2} ∈ [2p/ (p + q),2]. The exact values depend on the masses of both planets. The lower bounds correspond to the case m_{1} ≪ m_{2} and the upper bounds to m_{2} ≪ m_{1}. We thus always have (see Eq. (60)) (73)whereas γ_{2} can either be greater or smaller than γ_{c} depending on the considered resonance and the masses of the planets. If γ_{2}>γ_{c}, then γ>γ_{c} regardless of the balance of dissipation in both planets (τ). If γ_{2}<γ_{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 semimajor axis ratio, α = a_{1}/a_{2}, just after the crossing can be estimated from Eq. (63), (75)We introduce (76)Then, if τ>τ_{α}, the semimajor axis ratio decreases (diverging orbits), whereas if τ<τ_{α} the semimajor 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_{1}/T_{2}) >τ_{c}, the amplitude of libration decreases while the eccentricities are being damped. For firstorder 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_{2}/P_{1}< (p + q) /p),

if τ<τ_{α}, the orbits are diverging (P_{2}/P_{1}> (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_{1}/e_{2} 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_{2}/P_{1}) since the semimajor 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.
Orbital parameters of the GJ 163 system (from Bonfils et al. 2013, fit with tidal constraint), and the GJ 581 system (from Forveille et al. 2011).
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 Mdwarf 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_{2}/P_{1}< 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_{2}/P_{1}< (p + q) /p, the ratio τ = T_{1}/T_{2} 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_{2}/ Δt_{1} 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_{1}/e_{2} ≈ 1.1 at the libration center. The criterion τ>τ_{α} implies Δt_{2}/κΔt_{1} ≳ 300. The criterion τ<τ_{c} implies Δt_{2}/κΔt_{1} ≲ 1450. Thus, the scenario that we described should be possible if (84)
4.1.3. Nbody 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^{3} − 10^{4}, 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^{3} − 10^{4} 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_{1} ~ 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_{1} = 10^{7}, 10^{6}, 10^{5}, and 10^{4} 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_{1} = 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_{1} = 0°, the system is initially at the libration center while for M_{1} = 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 antialigned, and M_{2} = 0°. We took a_{1} = 0.062 such that when eccentricities go to zero the system ends up approximately at its current position (a_{1} ≈ 0.0607).
Figure 4 shows the final period ratio of the planets as a function of the lag time ratio with M_{1} = 100° (top) and 140° (bottom). In both cases, we superimposed the results obtained with the four different dissipative timescales (Δt_{1} = 10^{7}, 10^{6}, 10^{5}, and 10^{4} 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_{2}/κΔt_{1}. 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_{2}/κΔt_{1} ≈ 250 for both initial amplitudes of libration (M_{1} = 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_{2}/κΔt_{1} ≈ 750 for M_{1} = 100°, and 900 for M_{1} = 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_{1} → 180°).
Fig. 3 Diagram of initial conditions used in numerical simulations of GJ 163 (3:1 MMR). M_{1} 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_{1}). 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_{1} = 10^{7} 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_{2}/κΔt_{1} ≈ 250.
Fig. 4 Final period ratio of GJ 163’s inner planets as a function of the dissipation balance between both planets (Δt_{2}/κΔt_{1}). We fix M_{1} = 100° (left), 140° (right), and compare the results obtained with Δt_{1} = 10^{7}, 10^{6}, 10^{5}, and 10^{4} 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 Δt_{1} (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 Δt_{1} = 10^{7} s. We observe only small variations when changing the timescale of the dissipation (Δt_{1}) and no systematic trend. The third planet also does not seem to have a strong effect. The vertical black lines mark the range of Δt_{2}/κΔt_{1} 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 Δt_{2}/κΔt_{1} = 260 in the case M_{1} = 100° (left) and 240 for M_{1} = 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_{2}/κΔt_{1} ~ 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_{1}/e_{2} 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_{1} = 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_{2}/ Δt_{1} ~ 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_{1}/R_{2})^{5}. The Love numbers of both planets should be of the same order of magnitude^{1}. 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 xaxis of Fig. 6 is scaled using this value of κ. Finally, we conclude that Δt_{2}/ Δt_{1} ~ 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 M_{1} initially set to 100° and Δt_{2}/κΔt_{1} = 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 Δt_{1} = 10^{5} s in all these simulations. For a more realistic dissipation (Δt_{1} = 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 Mdwarf and both planets b and c, which we are interested in, are uncontested. We reproduced the orbital parameters of the fourplanet 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_{2}/κΔt_{1} ∈ [3,52] (with e_{1}/e_{2} ≈ 0.25). As we did for GJ 163, we performed numerical simulations of the system for different lag time ratios (Δt_{2}/κΔt_{1}) and different initial amplitude of libration (M_{1}). 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 antialigned, and M_{2} = −36°. We took a_{1} = 0.042 such that when eccentricities are damped the system ends up approximately at its current position (a_{1} ≈ 0.041). For M_{1} = 0°, the system is initially at the center of libration, while for M_{1} = 90° the system begins on the separatrix. The system is integrated during 5 × 10^{4} yr with Δt_{1} = 10^{7} s. The age of the system is estimated to be about 8 Gyr, thus our simulations are approximately equivalent to Δt_{1} ~ 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_{2}/ Δt_{1} ~ 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_{1}/e_{2} ≈ 1.1, while here we have e_{1}/e_{2} ≈ 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_{1}/e_{2} increases with m_{2}/m_{1} (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 (Δt_{2}/κΔt_{1}, bottom xaxis) and initial amplitude of libration (M_{1}). The coefficient κ = k_{2,1}/k_{2,2}(R_{1}/R_{2})^{5} is unknown, but can be estimated using a massradius power law. The top xaxis (Δt_{2}/ Δt_{1}) is scaled using the value κ ≈ 2.2 obtained in Eq. (86). When M_{1} is set to 0°, the system begins at the libration center, whereas when M_{1} is set to 180° it starts at the separatrix. We fixed Δt_{1} = 10^{7} 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 Δt_{1} (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 Δt_{2}/κΔt_{1} 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 xaxis (Δt_{2}/ Δt_{1}) is scaled using the value κ ≈ 15 obtained in Eq. (87). We note that M_{1} varies only between 0° (center of libration) and 90° (separatrix), because the resonant combination is 5M_{2} − 2M_{1} (factor 2 in front of M_{1}). 
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 xaxis of Fig. 7 is scaled using this value of κ. The lag time ratio of this system should be Δt_{2}/ Δt_{1} ~ 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
Multiplanetary systems detected by the Kepler mission present an excess of planet pairs close to firstorder 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 firstorder 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_{1} undergoes an exponential decrease e_{1} = e_{1,0}e^{− t/T1}, with e_{1,0} the eccentricity of the inner planet when the system crosses the separatrix, and T_{1} 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 (e_{1,0}). is the period ratio at the exact resonance. 
Lee et al. (2013) estimated that tidal dissipation would need at least t ≳ 50T_{1} 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_{1}.
Fig. 9 Time required to produce a departure of P_{2}/P_{1} of 0.03 from the resonant value as a function of the eccentricity of the inner planet when the system leaves the resonance (e_{1,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/T_{1} ≳ 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_{1}) 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_{1}). On the contrary, when e_{1,0} ≳ 0.15, the system reaches the desired period ratio on a much shorter timescale: t/T_{1} ~ 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 nearresonant 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 nonnegligible 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 twoplanet meanmotion 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_{1}/T_{2} or Δt_{2}/ Δt_{1}) and the position of the libration center (especially the ratio e_{1}/e_{2}).
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_{2}/ Δt_{1} ~ 1000. Using Nbody simulations with dissipation we reach the same conclusion with slightly refined bounds for Δt_{2}/ Δt_{1}. 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_{2}/ Δt_{1} ~ 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 highorder resonances, which may show large chaotic areas. Moreover, we make all our estimates using a constant eccentricity ratio (e_{1}/e_{2}) 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 nonnegligible (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 firstorder 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 nonnegligible eccentricities (e_{1} ≳ 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.
In the solar system, the Love number (k_{2}) of the Earth is 0.29 (Kozai 1968), Jupiter 0.379, Saturn 0.341 (Gavrilov & Zharkov 1977).
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 PNPCNRS, CS of Paris Observatory, PICS05998 FrancePortugal program, and FCTPortugal (PEstC/CTM/LA0025/2011).
References
 Antoniadou, K. I., & Voyatzis, G. 2014, Astrophys. Space Sci., 349, 657 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Morbidelli, A. 2013a, A&A, 556, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Batygin, K., & Morbidelli, A. 2013b, AJ, 145, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Beaugé, C., Michtchenko, T. A., & FerrazMello, S. 2006, MNRAS, 365, 1160 [Google Scholar]
 Bonfils, X., Forveille, T., Delfosse, X., et al. 2005, A&A, 443, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonfils, X., Lo Curto, G., Correia, A. C. M., et al. 2013, A&A, 556, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Callegari, N., Michtchenko, T. A., & FerrazMello, S. 2004, Celest. Mech. Dyn. Astron., 89, 201 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Callegari, N., FerrazMello, S., & Michtchenko, T. A. 2006, Celest. Mech. Dyn. Astron., 94, 381 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dyn. Astron., 111, 105 [Google Scholar]
 Delisle, J.B., Laskar, J., Correia, A. C. M., & Boué, G. 2012, A&A, 546, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2012, ApJ, submitted [arXiv:1202.6328] [Google Scholar]
 Forveille, T., Bonfils, X., Delfosse, X., et al. 2011, A&A, submitted [arXiv:1109.2505] [Google Scholar]
 Gavrilov, S. V., & Zharkov, V. N. 1977, Icarus, 32, 443 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Schlichting, H. E. 2014, AJ, 147, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Hadjidemetriou, J. D. 2002, Celest. Mech. Dyn. Astron., 83, 141 [Google Scholar]
 Hairer, E., Nørsett, S. P., & Wanner, G. 2010, Solving Ordinary Differential Equations I: Nonstiff Problems (Springer) [Google Scholar]
 Henrard, J., & Lemaitre, A. 1983, Celest. Mech., 30, 197 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Henrard, J., Milani, A., Murray, C. D., & Lemaitre, A. 1986, Celest. Mech., 38, 335 [Google Scholar]
 Kozai, Y. 1968, PASJ, 20, 24 [NASA ADS] [Google Scholar]
 Laskar, J. 2000, Phys. Rev. Lett., 84, 3240 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lee, M. H., Fabrycky, D., & Lin, D. N. C. 2013, ApJ, 774, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Lithwick, Y., & Wu, Y. 2012, ApJ, 756, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Bonfils, X., Forveille, T., et al. 2009, A&A, 507, 487 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michtchenko, T. A., & FerrazMello, S. 2001, Icarus, 149, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Michtchenko, T. A., Beaugé, C., & FerrazMello, S. 2006, Celest. Mech. Dyn. Astron., 94, 411 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Mignard, F. 1979, Moon and Planets, 20, 301 [Google Scholar]
 Novak, G. S., Lai, D., & Lin, D. N. C. 2003, in Scientific Frontiers in Research on Extrasolar Planets, eds. D. Deming, & S. Seager, ASP Conf. Ser., 294, 177 [Google Scholar]
 Papaloizou, J. C. B. 2011, Celest. Mech. Dyn. Astron., 111, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Terquem, C. 2010, MNRAS, 405, 573 [NASA ADS] [Google Scholar]
 Rein, H. 2012, MNRAS, 427, L21 [NASA ADS] [Google Scholar]
 Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sessin, W., & FerrazMello, S. 1984, Celest. Mech., 32, 307 [Google Scholar]
 Singer, S. F. 1968, Geophys. J. R. Astron. Soc., 15, 205 [CrossRef] [Google Scholar]
 Udry, S., Bonfils, X., Delfosse, X., et al. 2007, A&A, 469, L43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Wisdom, J. 1986, Celest. Mech., 38, 175 [Google Scholar]
Appendix A: Evolution of the amplitude of libration
The aim of this appendix is to prove Prop. 3.1. Assuming δ^{2} ≫ Rδ^{q/ 2}, the Hamiltonian (50) can be approximated with the pendulumlike Hamiltonian (A.1)Noting D_{1} = δ + ϵ, 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} = δ^{2} + 2Rδ^{q/ 2}. The hyperbolic fixed point is at qθ_{1,hyp} = π, ϵ_{hyp} = 0, ℋ_{hyp} = δ^{2} − 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 θ_{1} 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 longterm evolution of A. We need to compute the mean values of ϵ, ϵ^{2}, and ϵ^{3}. Because of the symmetry of the problem, ϵ and ϵ^{3} 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
Orbital parameters of the GJ 163 system (from Bonfils et al. 2013, fit with tidal constraint), and the GJ 581 system (from Forveille et al. 2011).
All Figures
Fig. 1 Position of the fixed points of a secondorder resonance (q = 2) in the simplified model (Hamiltonian (49)) as a function of δ. We give the real (top) and imaginary (bottom) parts of u_{1} 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 secondorder resonance. We plot the energy levels in the complex plane (real and imaginary parts of u_{1}), 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 C_{e} (respectively C_{i}) to external (respectively internal) circulation. 

In the text 
Fig. 3 Diagram of initial conditions used in numerical simulations of GJ 163 (3:1 MMR). M_{1} 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 (Δt_{2}/κΔt_{1}). We fix M_{1} = 100° (left), 140° (right), and compare the results obtained with Δt_{1} = 10^{7}, 10^{6}, 10^{5}, and 10^{4} 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 Δt_{1} (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 Δt_{1} = 10^{7} s. We observe only small variations when changing the timescale of the dissipation (Δt_{1}) and no systematic trend. The third planet also does not seem to have a strong effect. The vertical black lines mark the range of Δt_{2}/κΔt_{1} 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 Δt_{2}/κΔt_{1} = 260 in the case M_{1} = 100° (left) and 240 for M_{1} = 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 M_{1} initially set to 100° and Δt_{2}/κΔt_{1} = 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 Δt_{1} = 10^{5} s in all these simulations. For a more realistic dissipation (Δt_{1} = 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 (Δt_{2}/κΔt_{1}, bottom xaxis) and initial amplitude of libration (M_{1}). The coefficient κ = k_{2,1}/k_{2,2}(R_{1}/R_{2})^{5} is unknown, but can be estimated using a massradius power law. The top xaxis (Δt_{2}/ Δt_{1}) is scaled using the value κ ≈ 2.2 obtained in Eq. (86). When M_{1} is set to 0°, the system begins at the libration center, whereas when M_{1} is set to 180° it starts at the separatrix. We fixed Δt_{1} = 10^{7} 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 Δt_{1} (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 Δt_{2}/κΔt_{1} 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 xaxis (Δt_{2}/ Δt_{1}) is scaled using the value κ ≈ 15 obtained in Eq. (87). We note that M_{1} varies only between 0° (center of libration) and 90° (separatrix), because the resonant combination is 5M_{2} − 2M_{1} (factor 2 in front of M_{1}). 

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 (e_{1,0}). is the period ratio at the exact resonance. 

In the text 
Fig. 9 Time required to produce a departure of P_{2}/P_{1} of 0.03 from the resonant value as a function of the eccentricity of the inner planet when the system leaves the resonance (e_{1,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/T_{1} ≳ 50 − 1000). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.