© ESO, 2017

1. Introduction

Mean motion resonances (MMR) between two planets are a natural outcome of the convergent migration of planets in a gas-disk (e.g., Weidenschilling & Davis 1985). The planets initially form farther away from each other, and planet-disk interactions induce a migration of the planets. The period ratio between the planets decreases until they get captured in a MMR. The planets then continue to migrate whilst maintaining their period ratio at a rational value (2/1, 3/2, etc.). The eccentricities increase due to the resonant interactions, until they reach an equilibrium between the migration torque and the eccentricity damping exerted by the disk. The argument of the resonance, which is a combination of the mean longitudes of the two planets, enters into libration (oscillations around an equilibrium value).

For systems of three and more planets, once a pair of planets has been captured in a MMR, the other planets might also join this couple to form a chain of resonances. Each time a planet gets captured in the chain, it enters into a MMR (and thus maintains a constant and rational period ratio) with each of the other planets of the chain. The eccentricities of the planets and the resonant arguments of each pair find a new equilibrium. Such multi-planetary resonant chains are expected from simulations of planet migration (e.g., Cresswell & Nelson 2006). Recently, Mills et al. (2016) showed that the four planets in the Kepler-223 system are in a 3:4:6:8 resonant chain (period ratios of 4/3, 3/2, and 4/3 between consecutive pairs of planets). Using transit timing variations (TTVs), the authors observed that the Laplace angles of the system are librating with small amplitudes. The Laplace angles are combinations of the mean longitudes of three planets in the chain, and the observation of their libration is evidence that the system is indeed captured in the resonant chain. Using numerical simulations, Mills et al. (2016) showed that the observed orbital configuration is very well reproduced by a smooth convergent migration of the planets.

In this article, I present an analytical model of resonant chains. Analytical models have already been proposed, in particular to study the dynamics of the Laplace resonance (1:2:4 chain) between Io, Europa, and Ganymede (e.g., Henrard 1983). However, while several numerical studies have been dedicated to the capture of planets in various resonant chains (e.g., Cresswell & Nelson 2006; Papaloizou & Terquem 2010; Libert & Tsiganis 2011; Papaloizou 2016), general analytical models have not yet been proposed. Recently, Papaloizou (2015) proposed a semi-analytical model of three-planet resonances taking into account only the interactions between consecutive planets in the chain, with a particular focus on the Kepler-60 system (12:15:20 resonant chain, see also Steffen et al. 2013; Goździewski et al. 2016). This model is very similar to the studies of the Laplace resonance between the Galilean moons, but is not well suited in the general case. For instance, four-planet (or more) resonances are not considered. Moreover, for some three-planet resonances, the interactions between non-consecutive planets cannot be neglected. For instance, in a 3:4:6 resonant chain, each planet is locked in a first-order resonance with each of the other planets. In particular, the innermost and outermost planets are involved in a 2/1 MMR that strongly influences the dynamics of the system. I describe here a general model of resonant chains, with any number of planets, valid for any resonance order. I particularly focus on finding the equilibrium configurations (eccentricities, resonant arguments, etc.) around which a resonant system should librate. While a real system may be observed with significant amplitude of libration around the equilibrium, or could even have some angles circulating, the position of the equilibria still provides useful insights into the dynamics of the system. In Sect. 2, I describe this analytical model, and the method I use to find the equilibrium configurations. In Sect. 3, I apply the model to Kepler-223. I show that six equilibrium configurations exist for this resonant chain, and that the system is observed to be librating around one of them. I also show that knowing the current configuration of the system allows for interesting constraints to be put on its migration scenario, and in particular on the order in which the planets have been captured in the chain.

2. Model

I consider a planetary system with n planets (which I denote with indices 1,...,n from the innermost to the outermost) orbiting around a star (index 0). I assume that the system is coplanar and is locked in a chain of resonances. In such a resonant chain, each pair of planets is locked in a MMR. For two planets i<j, I denote by k_j,i/k_i,j the resonant ratio, such that ${\mathit{k}}_{\mathit{j,i}}{\mathit{n}}_{\mathit{j}}\mathrm{-}{\mathit{k}}_{\mathit{i,j}}{\mathit{n}}_{\mathit{i}}\mathrm{\approx }\mathrm{0}\mathit{,}$(1)where n_i (n_j) is the mean motion of planet i (j). I also introduce the degree of the resonance between planet i and planet j${\mathit{q}}_{\mathit{i,j}}\mathrm{=}{\mathit{k}}_{\mathit{j,i}}\mathrm{-}{\mathit{k}}_{\mathit{i,j}}\mathit{.}$(2)At low eccentricities, resonances of a lower degree have a stronger influence on the dynamics of the system.

In order to study the dynamics of these resonant chains, I generalize to n planets the method developped in the case of two-planet resonances (Delisle et al. 2012, 2014). The Hamiltonian of the system takes the form (Laskar 1991) $\mathrm{\mathscr{H}}\mathrm{=}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{𝒢}\frac{{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathit{i}}}{\mathrm{2}{\mathit{a}}_{\mathit{i}}}\mathrm{+}\sum _{\mathrm{1}\mathrm{\le }\mathit{i}\mathit{<}\mathit{j}\mathrm{\le }\mathit{n}}\left(\mathrm{-}\mathrm{𝒢}\frac{{\mathit{m}}_{\mathit{i}}{\mathit{m}}_{\mathit{j}}}{\mathrm{\left|}\mathrm{\right|}{{r}}_{\mathit{i}}\mathrm{-}{{r}}_{\mathit{j}}\mathrm{\left|}\mathrm{\right|}}\mathrm{+}\frac{{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathit{i}}\mathrm{·}{\begin{array}{c}{˜}\\ {r}\end{array}}_{\mathit{j}}}{{\mathit{m}}_{\mathrm{0}}}\right)\mathit{,}$(3)where is the gravitational constant, m_i is the mass of body i, a_i is the semi-major axis, r_i the position vector, and ${\begin{array}{c}{{˜}}_{}\\ {r}\end{array}}_{\mathit{i}}$ the canonically conjugated momentum of planet i (in astrocentric coordinates, see Laskar 1991). The first sum on the right-hand side of Eq. (3) is the Keplerian part of the Hamiltonian (planet-star interactions), while the second sum is the perturbative part (planet-planet interactions).

In the coplanar case (which I assume here) the system has 2n degrees of freedom (d.o.f.), with 2 d.o.f. (4 coordinates) associated to each planet. As for two-planet resonances (e.g., Delisle et al. 2012), the number of d.o.f. can be reduced by using the conservation of the total angular momentum (1 d.o.f.), and by averaging over the fast angles (1 d.o.f.). Therefore, the problem can be reduced to 2(n−1) d.o.f. Even with these reductions, the phase space is still very complex, especially for systems of many planets such as Kepler-223 (chain of 4 planets, 6 d.o.f.), and the problem is, in most cases, non-integrable. In this study, I focus on finding the fixed point of the averaged problem, which provides useful insight into the dynamics of the system, and especially into the values around which the angles of a resonant system should librate. The method described in the following is a generalization of the method presented in Delisle et al. (2012) which focuses on finding the fixed points for two-planet MMR.

I denote by λ_i and ϖ_i the mean longitude and longitude of periastron of planet i (in astrocentric coordinates), respectively. The actions canonically conjugated to the angles λ_i and −ϖ_i are the circular angular momentum Λ_i and the angular momentum deficit (AMD, see Laskar 2000) D_i, respectively. These actions are defined as follows $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathit{i}}& \mathrm{=}& {\mathit{\beta }}_{\mathit{i}}\sqrt{{\mathit{\mu }}_{\mathit{i}}{\mathit{a}}_{\mathit{i}}}\mathit{,}\\ {\mathit{D}}_{\mathit{i}}& \mathrm{=}& {\mathrm{\Lambda }}_{\mathit{i}}\mathrm{-}{\mathit{G}}_{\mathit{i}}\mathrm{=}{\mathrm{\Lambda }}_{\mathit{i}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}\right)\mathit{,}\end{array}$where ${\mathit{G}}_{\mathit{i}}\mathrm{=}{\mathrm{\Lambda }}_{\mathit{i}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}$ is the angular momentum of planet i, β_i = m_im₀/ (m₀ + m_i), . At low eccentricities the deficit of angular momentum D_i is proportional to ${\mathit{e}}_{\mathit{i}}^{\mathrm{2}}$. The Hamiltonian (Eq. (3)) can be expressed using these action-angle coordinates $\mathrm{\mathscr{H}}\mathrm{=}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{\mathit{\mu }}_{\mathit{i}}^{\mathrm{2}}{\mathit{\beta }}_{\mathit{i}}^{\mathrm{3}}}{\mathrm{2}{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{2}}}\mathrm{+}\sum _{\mathrm{1}\mathrm{\le }\mathit{i}\mathit{<}\mathit{j}\mathrm{\le }\mathit{n}}{\mathrm{\mathscr{H}}}_{\mathit{i,j}}\mathrm{\left(}{\mathrm{\Lambda }}_{\mathit{i}}\mathit{,}{\mathrm{\Lambda }}_{\mathit{j}}\mathit{,}{\mathit{D}}_{\mathit{i}}\mathit{,}{\mathit{D}}_{\mathit{j}}\mathit{,}{\mathit{\lambda }}_{\mathit{i}}\mathit{,}{\mathit{\lambda }}_{\mathit{j}}\mathit{,}{\mathit{\varpi }}_{\mathit{i}}\mathit{,}{\mathit{\varpi }}_{\mathit{j}}\mathrm{\right)}\mathit{,}$(6)where the first sum is the Keplerian part, which depends only on Λ_i (or equivalently a_i), and ℋ_i,j is the perturbation between the planets i and j, which depends on the eight coordinates associated to i and j. I follow the method described in Laskar & Robutel (1995) to compute H_i,j as a power series of the eccentricities (or equivalently of $\sqrt{{\mathit{D}}_{\mathit{i}}}$ and $\sqrt{{\mathit{D}}_{\mathit{j}}}$), and a Fourrier series of the angles, where the coefficients are functions of Λ_i and Λ_j (i.e., of the semi-major axes). For a system that is close to the resonance or resonant, the semi-major axes remain close to the nominal resonant values (Kepler’s third law) $\frac{{\mathit{a}}_{\mathit{i}}}{{\mathit{a}}_{\mathit{j}}}\mathrm{\approx }\frac{{\mathit{a}}_{\mathit{i,}\mathrm{0}}}{{\mathit{a}}_{\mathit{j,}\mathrm{0}}}\mathrm{=}{\left(\frac{{\mathit{k}}_{\mathit{i,j}}}{{\mathit{k}}_{\mathit{j,i}}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}{\left(\frac{{\mathit{\mu }}_{\mathit{i}}}{{\mathit{\mu }}_{\mathit{j}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}$(7)I introduce $\mathrm{\Delta }{\mathrm{\Lambda }}_{\mathit{i}}\mathrm{=}{\mathrm{\Lambda }}_{\mathit{i}}\mathrm{-}{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}\mathit{,}$(8)where ${\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}\mathrm{=}{\mathit{\beta }}_{\mathit{i}}\sqrt{{\mathit{\mu }}_{\mathit{i}}{\mathit{a}}_{\mathit{i,}\mathrm{0}}}\mathit{,}$(9)and expand the Keplerian part at degree 2, and the perturbative part at degree 0 in ΔΛ_i$\begin{array}{ccc}\mathrm{\mathscr{H}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{n}}_{\mathit{i,}\mathrm{0}}\mathrm{\Delta }{\mathrm{\Lambda }}_{\mathit{i}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\frac{{\mathit{n}}_{\mathit{i,}\mathrm{0}}}{{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}}\mathrm{\Delta }{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{2}}\\ & & \mathrm{+}\sum _{\mathrm{1}\mathrm{\le }\mathit{i}\mathit{<}\mathit{j}\mathrm{\le }\mathit{n}}{\mathrm{\mathscr{H}}}_{\mathit{i,j}}\mathrm{\left(}{\mathit{D}}_{\mathit{i}}\mathit{,}{\mathit{D}}_{\mathit{j}}\mathit{,}{\mathit{\lambda }}_{\mathit{i}}\mathit{,}{\mathit{\lambda }}_{\mathit{j}}\mathit{,}{\mathit{\varpi }}_{\mathit{i}}\mathit{,}{\mathit{\varpi }}_{\mathit{j}}\mathrm{\right)}\mathit{,}\end{array}$(10)where n_i,0 is the nominal mean motion of planet i, such that $\frac{{\mathit{n}}_{\mathit{i,}\mathrm{0}}}{{\mathit{n}}_{\mathit{j,}\mathrm{0}}}\mathrm{=}\frac{{\mathit{k}}_{\mathit{j,i}}}{{\mathit{k}}_{\mathit{i,j}}}\mathrm{·}$(11)The perturbative part does not depend on Λ_i anymore, but is simply evaluated at Λ_i,0.

In order to perform the reductions associated to the conservation of angular momentum and to the averaging, I first change the system of coordinates. For the sake of readability, I present the general case (with any number of planets, in any resonance of any degree) in Appendix A, and take here the example of a system of four planets in a 3:4:6:8 resonant chain (as is the case for Kepler-223). I introduce new canonically conjugated angles and actions as follows (see Eqs. (A.1) and (A.2)) $\begin{array}{ccc}& & {\mathit{\phi }}_{\mathrm{1}}\mathrm{=}{\mathit{\lambda }}_{\mathrm{1}}\mathrm{+}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{-}\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\mathit{,} {\mathit{L}}_{\mathrm{1}}\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{1}}\mathit{,}\\ & & {\mathit{\phi }}_{\mathrm{2}}\mathrm{=}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{\lambda }}_{\mathrm{4}}\mathrm{-}\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}\mathit{,} {\mathit{L}}_{\mathrm{2}}\mathrm{=}\mathrm{2}{\mathrm{\Lambda }}_{\mathrm{1}}\mathrm{+}{\mathrm{\Lambda }}_{\mathrm{2}}\mathit{,}\\ & & {\mathit{\phi }}_{\mathrm{3}}\mathrm{=}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{-}{\mathit{\lambda }}_{\mathrm{4}}\mathit{,} \mathrm{\Gamma }\mathrm{=}\mathrm{8}{\mathrm{\Lambda }}_{\mathrm{1}}\mathrm{+}\mathrm{6}{\mathrm{\Lambda }}_{\mathrm{2}}\mathrm{+}\mathrm{4}{\mathrm{\Lambda }}_{\mathrm{3}}\mathrm{+}\mathrm{3}{\mathrm{\Lambda }}_{\mathrm{4}}\mathit{,}\\ & & {\mathit{\phi }}_{\mathrm{4}}\mathrm{=}\mathrm{4}{\mathit{\lambda }}_{\mathrm{4}}\mathrm{-}\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}\mathit{,} \mathit{G}\mathrm{=}{\mathit{G}}_{\mathrm{1}}\mathrm{+}{\mathit{G}}_{\mathrm{2}}\mathrm{+}{\mathit{G}}_{\mathrm{3}}\mathrm{+}{\mathit{G}}_{\mathrm{4}}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathrm{1}}\mathrm{=}\mathrm{4}{\mathit{\lambda }}_{\mathrm{4}}\mathrm{-}\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{1}}\mathit{,} {\mathit{D}}_{\mathrm{1}}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathrm{2}}\mathrm{=}\mathrm{4}{\mathit{\lambda }}_{\mathrm{4}}\mathrm{-}\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{2}}\mathit{,} {\mathit{D}}_{\mathrm{2}}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathrm{3}}\mathrm{=}\mathrm{4}{\mathit{\lambda }}_{\mathrm{4}}\mathrm{-}\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{3}}\mathit{,} {\mathit{D}}_{\mathrm{3}}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathrm{4}}\mathrm{=}\mathrm{4}{\mathit{\lambda }}_{\mathrm{4}}\mathrm{-}\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{4}}\mathit{,} {\mathit{D}}_{\mathrm{4}}\mathit{.}\end{array}$(12)G is the total angular momentum, and is a conserved quantity; Its canonically conjugated angle (φ₄) does not appear in the Hamiltonian. The angle φ₃ is the only fast angle, and the averaging of the Hamiltonian is done over this angle. Therefore, its conjugated action (Γ) is constant in the average problem. The averaging is simply done by discarding all terms that depends on φ₃ in the Fourrier expansion of the Hamiltonian¹. The angle φ₁ is the argument of the Laplace resonance between the three innermost planets. The angle φ₂ is the argument of the Laplace resonance between the three outermost planets.

The two outer planets (3 and 4) may seem to play an important role in Eq. (12), but this is only due to the arbitrary choice of canonical coordinates (many other choices are possible). Any two-planet resonant angle can be expressed as a combination of the angles of Eq. (12). The arguments of the 4/3 resonance between the two outermost planets are σ₃ and σ₄. The arguments of the 3/2 resonance between planets 2 and 3 are σ₂−2φ₂ = 3λ₃−2λ₂−ϖ₂ and σ₃−2φ₂ = 3λ₃−2λ₂−ϖ₃. The arguments of the 4/3 resonance between planets 1 and 2 are σ₁−2φ₂−3φ₁ = 4λ₂−3λ₁−ϖ₁ and σ₂−2φ₂−3φ₁ = 4λ₂−3λ₁−ϖ₂. Arguments of resonances between non-consecutive pairs can also be expressed in the same way. For instance, for the 2/1 resonance between planets 1 and 3, the arguments are σ₁−2φ₂−φ₁ = 2λ₃−λ₁−ϖ₁ and σ₃−2φ₂−φ₁ = 2λ₃−λ₁−ϖ₃. For a system of n planets captured in a resonant chain, φ_i (i ≤ n−2) and σ_i (i ≤ n) (and all their linear combinations) librate around equilibrium values. All the actions also oscillate around equilibria. These equilibrium values correspond to stable fixed points of the average problem.

For this example, I expand the perturbative part at first order in eccentricities ($\sqrt{{\mathit{D}}_{\mathit{i}}}$), and obtain an expression of the form $\begin{array}{ccc}\mathrm{\mathscr{H}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{4}}\frac{{\mathit{n}}_{\mathit{i,}\mathrm{0}}}{{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}}\mathrm{\Delta }{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{2}}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{1}\mathit{,}\mathrm{2}}\sqrt{{\mathit{D}}_{\mathrm{1}}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathrm{1}}\mathrm{-}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{-}\mathrm{3}{\mathit{\phi }}_{\mathrm{1}}\mathrm{)}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{1}}\sqrt{{\mathit{D}}_{\mathrm{2}}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{-}\mathrm{3}{\mathit{\phi }}_{\mathrm{1}}\mathrm{)}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{1}\mathit{,}\mathrm{3}}\sqrt{{\mathit{D}}_{\mathrm{1}}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathrm{1}}\mathrm{-}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{-}{\mathit{\phi }}_{\mathrm{1}}\mathrm{)}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{1}}\sqrt{{\mathit{D}}_{\mathrm{3}}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathrm{3}}\mathrm{-}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{-}{\mathit{\phi }}_{\mathrm{1}}\mathrm{)}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{3}}\sqrt{{\mathit{D}}_{\mathrm{2}}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{2}}\sqrt{{\mathit{D}}_{\mathrm{3}}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathrm{3}}\mathrm{-}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{)}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{4}}\sqrt{{\mathit{D}}_{\mathrm{2}}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathrm{2}}\mathrm{-}{\mathit{\phi }}_{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{C}}_{\mathrm{4}\mathit{,}\mathrm{2}}\sqrt{{\mathit{D}}_{\mathrm{4}}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathrm{4}}\mathrm{-}{\mathit{\phi }}_{\mathrm{2}}\mathrm{)}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{4}}\sqrt{{\mathit{D}}_{\mathrm{3}}}\cos \mathrm{\left(}{\mathit{\sigma }}_{\mathrm{3}}\mathrm{\right)}\mathrm{+}{\mathit{C}}_{\mathrm{4}\mathit{,}\mathrm{3}}\sqrt{{\mathit{D}}_{\mathrm{4}}}\cos \mathrm{\left(}{\mathit{\sigma }}_{\mathrm{4}}\mathrm{\right)}\mathit{,}\end{array}$(13)where the first term of Eq. (10) vanishes (because Γ is constant), and C_i,j are constant coefficients that depend on the masses and nominal semi-major axes (a_i,0). I provide explicit formulas for the case of the 3:4:6:8 resonant chain in Appendix B. Since the Hamiltonian (Eq. (13)) is developed at first order in eccentricities, only first-order resonances appear. In particular, the 3/8 resonance between planets 1 and 4 is neglected since it would only appear at order 5 in eccentricities. In order to use a consistent (canonical) set of coordinates, ΔΛ_i should be replaced in Eq. (13) by $\begin{array}{ccc}\mathrm{\Delta }{\mathrm{\Lambda }}_{\mathrm{1}}& \mathrm{=}& \mathrm{\Delta }{\mathit{L}}_{\mathrm{1}}\mathit{,}\\ \mathrm{\Delta }{\mathrm{\Lambda }}_{\mathrm{2}}& \mathrm{=}& \mathrm{\Delta }{\mathit{L}}_{\mathrm{2}}\mathrm{-}\mathrm{2}\mathrm{\Delta }{\mathit{L}}_{\mathrm{1}}\mathit{,}\\ \mathrm{\Delta }{\mathrm{\Lambda }}_{\mathrm{3}}& \mathrm{=}& \mathrm{\Delta }{\mathit{L}}_{\mathrm{1}}\mathrm{-}\mathrm{3}\mathrm{(}\mathrm{\Delta }{\mathit{L}}_{\mathrm{2}}\mathrm{+}\mathit{ϵ}\mathrm{)}\mathit{,}\\ \mathrm{\Delta }{\mathrm{\Lambda }}_{\mathrm{4}}& \mathrm{=}& \mathrm{2}\mathrm{\Delta }{\mathit{L}}_{\mathrm{2}}\mathrm{+}\mathrm{4}\mathit{ϵ,}\end{array}$(14)where $\mathit{ϵ}\mathrm{=}\mathit{D}\mathrm{-}\mathit{\delta }\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{\Delta }{\mathrm{\Lambda }}_{\mathit{i}}\mathit{,}$(15)which measures the distance of the system to the exact resonance, where$\mathit{D}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{D}}_{\mathit{i}}\mathit{,}$(16)is the total deficit of angular momentum, and $\mathit{\delta }\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}\mathrm{-}\mathit{G}$(17)is the nominal total deficit of angular momentum (at exact resonance). Since Λ_i,0 and G are constants, δ is also a conserved quantity and can be used as a parameter instead of G. The other parameter Γ does not appear explicitly in the Hamiltonian, but is hidden in the values of Λ_i,0, a_i,0, and n_i,0. Indeed, the nominal semi-major axis ratios are fixed at the resonant values (Eq. (7)), but Γ sets the global scale of the system (Γ = 8Λ_1,0 + 6Λ_2,0 + 4Λ_3,0 + 3Λ_4,0, for a 3:4:6:8 chain). The value of Γ does not influence the dynamics of the system appart from changing the scales of distance, time, and energy (Delisle et al. 2012). Therefore, one only need to vary the value of δ/ Γ to study the evolution of the phase space (in particular the positions of fixed points).

Fig. 1 Location of the (six) fixed points for the Kepler-223 resonant chain (4 planets), as a function of the parameter δ/ Γ (see Sect. 2). I plot the eccentricities of the planets (left), the difference of longitudes of periastron between two successive planets, as well as the two-planet resonant angle σ4 (center), and the Laplace angles (right). Each fixed point is represented by the same color in all the plots. The gray bands (on the left and right columns) show the observed parameters of the planets (taken from Mills et al. 2016). For the eccentricities (left), the gray bands correspond to the 1σ uncertainties. For the Laplace angles, they correspond to the observed libration amplitudes (which are well constrained by TTVs, see Mills et al. 2016). The observed positions of the planets (and especially the Laplace angles), are consistent with a libration around the red fixed point (see also Fig. 2).

For a given value of δ/ Γ, the fixed points are found by solving the following system of equations $\begin{array}{ccc}\mathit{\sigma ̇}\mathit{i}& \mathrm{=}& \frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }{\mathit{D}}_{\mathit{i}}}\mathrm{=}\mathrm{0} \mathrm{(}\mathit{i}\mathrm{\le }\mathit{n}\mathrm{)}\mathit{,}\\ \mathit{Ḋ}\mathit{i}& \mathrm{=}& \mathrm{-}\frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }{\mathit{\sigma }}_{\mathit{i}}}\mathrm{=}\mathrm{0} \mathrm{(}\mathit{i}\mathrm{\le }\mathit{n}\mathrm{)}\mathit{,}\\ \mathit{\phi ̇}\mathit{i}& \mathrm{=}& \frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }{\mathit{L}}_{\mathit{i}}}\mathrm{=}\mathrm{0} \mathrm{(}\mathit{i}\mathrm{\le }\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{,}\\ \mathit{L̇}\mathit{i}& \mathrm{=}& \mathrm{-}\frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }{\mathit{\phi }}_{\mathit{i}}}\mathrm{=}\mathrm{0} \mathrm{(}\mathit{i}\mathrm{\le }\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{.}\end{array}$(18)This is a system of 4(n−1) equations, with 4(n−1) unknowns (2(n−1) d.o.f.), which in general possesses a finite number of solutions. These solutions can correspond to elliptical (stable) fixed points or hyperbolic (unstable) ones. To assess the stability of fixed points, I compute the eigenvalues of the linearized equations of motions around the fixed point. Stable fixed points have purely imaginary eigenvalues, while the eigenvalues around unstable fixed points have a non-zero real part.

3. Application to Kepler-223

In this section, I apply my model to the four planets orbiting in a 3:4:6:8 resonant chain around Kepler-223. In Sect. 3.1, I focus on the comparison between the positions of the stable fixed points and the observed configuration of the system (values of resonant angles, eccentricities). In Sect. 3.2, I derive constraints on the order in which the planets have been captured in the resonant chain, from the observation of the equilibrium around which the system is currently librating.

3.1. Equilibrium configurations

I follow the method described in Sect. 2, and expand the Hamiltonian at first order in eccentricities (Eq. (13)), to solve for the positions of the fixed points (Eq. (18)) as a function of the parameter δ/ Γ (see Appendix C for more details). I only consider here the stable (elliptical) fixed points that correspond to the libration in resonance. In the case of Kepler-223, I find six families of stable fixed points (parameterized by δ/ Γ), corresponding to six possible areas of libration for the system. I show in Fig. 1 the positions of these fixed points, and compare them with the observed values (taken from Mills et al. 2016). As shown by Mills et al. (2016), the TTVs constrain the values of φ₁ and φ₂ very well, while the eccentricities are only roughly determined. The two-planet resonant angles σ_i (which depend on the longitudes of periastron) are not well constrained (Mills et al. 2016), so the theoretical values cannot be compared to the observations. As δ/ Γ increases, the eccentricities increase, but the angles (φ_i, σ_i) remain constant. Since the best observational constraints are on the values of φ₁ and φ₂, I compare in Fig. 2 the observed values of these angles, with the six possible equilibrium values (which are independent of δ/ Γ). The observations correspond very well to a libration of the system around one of the six configurations (the red one). This result confirms that the system is captured in the resonant chain, and that my model is correct in first approximation. It is interesting to wonder why the system was captured around this particular configuration and not one of the five others. This might simply be by chance (probability of 1/6 for each configuration), but might also greatly depend on the migration scenario. Observing around which equilibrium the resonant chain is librating could provide interesting constraints on the migration undergone by the planets.

3.2. Constraints on migration scenarios

As explained in Sect. 1, when two planets are captured in a MMR, their eccentricities increase until they reach an equilibrium (e.g., Delisle et al. 2012), their period ratio remains locked at the resonant value, and the arguments of the resonance librate. For a resonance between two planets, and at low eccentricities, the equilibrium configuration is unique. When additional planets join the chain, the system evolves toward a new equilibrium (eccentricities, two-planet resonant angles, Laplace angles). However, when more than two planets are involved in the chain, there can be more than one equilibrium (see Sect. 3.1), and the probability of capture around each of the new equilibria is not necessarily equal. Moreover, this probability might depend on the order in which the planets are captured in the chain. Indeed, the planets that are already captured in resonance have their eccentricities excited, and the angles associated to the already formed resonances are librating around an equilibrium, while planets that are not yet captured should have lower eccentricities, and their mean longitudes should be randomly distributed. Therefore, the initial conditions (angles, eccentricities) at the moment of the capture of a planet in the chain greatly depend on which planets are already in the chain.

Fig. 2 Comparison of the observed libration amplitudes of the Laplace angles in the Kepler-223 system (gray rectangle, values taken from Mills et al. 2016), with the positions of the (six) fixed points as determined from the analytical model. I use the same colors as in Fig. 1. The observations are consistent with a libration around the red fixed point.

In this section, I investigate how the order in which the planets are captured in the resonant chain influences the probability of capture around each of the six equilibria found in the case of Kepler-223 (see Sect. 3.1). This problem is very complex (the phase space has 6 d.o.f.). In particular, each time an additional planet joins the chain, the system might cross one or several separatrices before being captured around one of the new equilibria. I do not attempt to treat this problem analytically, but rather numerically estimate the probabilities of capture by running 6000 N-body simulations including prescriptions for the migration, with varying initial conditions. I use the same integrator as in Delisle et al. (2015), and set constant timescales for the migration torque and the eccentricity damping of each planet. The migration timescales (T_mig.,i) are set to 10⁶, 3 × 10⁵, 2 × 10⁵, and 1.5 × 10⁵ yr (from the innermost to the outermost planet). The eccentricity damping timescales (T_ecc.,i) are set such that for each planet T_mig.,i/T_ecc.,i = 50. The planets start with circular and coplanar orbits. The innermost planet is at 1 AU, and λ₁ = 0, and the other planets’ semi-major axes and mean longitudes are randomly drawn. The mean longitudes are drawn from a uniform distribution between 0 and 2π. The semi-major axis ratio a₂/a₁ is uniformly drawn in the range [1.23,1.3], a₃/a₂ in the range [1.32,1.38], and a₄/a₃ in the range [1.23,1.28]. This allows the planet pairs to be captured in MMR in any possible order. The orbits are integrated for 5 × 10⁴ yr.

Among the 6000 simulations, 5758 are captured in the 3:4:6:8 resonant chain, while 242 are captured in other chains (or not captured). All of the 5758 captured simulations are observed to librate around one of the six equilibria predicted by my analytical model. This confirms that this first-order model is correctly describing the dynamics of the resonant chain. For each simulation, I check in which order the planets were captured, and around which equilibrium the system ended.

The results are shown in Table 1. I use “A” to denote the capture of the two innermost planets (1, 2) in the 4/3 MMR, “B” the capture of the intermediate ones (2, 3) in the 3/2 MMR, and “C” the capture of the outermost planets (3, 4) in the 4/3 MMR. Thus, ABC means that the planets were captured in the chain from the innermost to the outermost one, and so on. For each possible order of capture (ordering of A, B, and C), I compute the percentage of capture around each equilibrium. If the captures were equally probable, the percentages would all be of about 16.7% (1/6). This is clearly not the case for Kepler-223. For instance, the black equilibrium is difficult to reach (low probability) except in the cases ACB and CAB (see Table 1). This means that the two innermost planets and the two outermost ones must first be captured in two independent two-planet resonances, and then the two pairs join to form the four-planet resonant chain.

Since the system is currently observed around the red equilibrium (see Figs. 1 and 2), it is interesting to look at the corresponding capture probabilities. The most favorable case is BCA, with a probability of 52% to capture the system around the observed equilibrium (Table 1). This case corresponds to a capture of planets 2 and 3 in the 3/2 MMR, then planet 4 joins the chain (2:3:4), and finally planet 1 is captured to form the 3:4:6:8 chain. Figure 3 shows an example of a simulation following this scenario, and reproducing the observed configuration of the Kepler-223 system. Conversely, in the case ABC (capture from the innermost to the outermost planet), the probability to reproduce the observed configuration is only 4.1%. The mean capture probability around the observed configuration (assuming equal probability for each capture order), is the highest value (25.2%). It is thus not surprising to observe the system in this configuration, rather than in one of the five others.

Fig. 3 Example of a simulation that successfully reproduces the observed configuration of the Kepler-223 system. The scenario of capture is of type BCA (see Sect. 3.2). I plot the period ratio between each consecutive pair of planets (top), the planets’ eccentricities (middle), and the two Laplace angles (bottom). The two horizontal black lines in the bottom plot represent the equilibrium values of the Laplace angles expected from the analytical model (see Sect. 3.1).

4. Discussion

In this article I describe an analytical model of resonant chains. The model is valid for any number of planets involved in the chain, and for any resonances (of any order). In particular, I use it to determine the equilibrium configurations around which a resonant chain librates. I show that contrarily to two-planet MMR, multiple equilibria may exist, even at low eccentricities, when three or more planets are involved.

I specifically study the case of the four planets around Kepler-223 which have been confirmed to be captured in a 3:4:6:8 resonant chain (using TTVs, see Mills et al. 2016). Using the analytical model expanded at first order in eccentricities, I show that six equilibrium configurations exist for this system, and the planets might have been captured around any of these six equilibria. However, the capture probabilities are not the same for each equilibrium, and depend on the order in which the planets have been captured in the chain. Using N-body integrations including migration prescriptions, I show that the scenario the most capable of reproducing the observed configuration of the system is to first capture the intermediate planets (2 and 3) in the 3/2 MMR, then to capture the outermost planet (4) to form a 2:3:4 chain between the three outermost planets, and finally to capture the innermost planet (1). This scenario of capture reproduces the observed configuration in 52% of the simulations, while capturing the planets from the innermost to the outermost reproduces the observed configuration for only 4.1% of the simulations. It should be noted that several hypotheses are made to compute these statistics. The planets are initially outside the resonances, with period ratios slightly higher than the resonant values. The migration and eccentricity damping timescales are fixed for each planet, such that the period ratio between each pair decreases (convergent migration). The planets are initially on circular and coplanar orbits. I only vary the initial semi-major axes and initial mean longitudes of the planets. The statistics would probably slightly change with a different set-up for the simulations, but this would not change the two main results of this study:

• 1.

  Resonant chains can be captured around several equilibriumconfigurations (six for Kepler-223).

• 2.

  observing a system around one of the possible equilibria provides useful constraints on the scenario of formation and migration of the planets.

Other properties of the resonant chain might provide useful additional constraints on such a scenario. For instance, in the case of Kepler-223, the two inner planets, as well as the two outer ones, are in a compact 4/3 resonance. If these planets formed with much wider separations, they must have avoided permanent capture into the 2/1 and 3/2 resonances to reach the currently observed 4/3 resonance. However, determining a complete scenario for the formation of the Kepler-223 system is a highly degenerated problem and is beyond the scope of this article.

Acknowledgments

I thank the anonymous referee for his/her useful comments. I acknowledge financial support from the Swiss National Science Foundation (SNSF). This work has, in part, been carried out within the framework of the National Centre for Competence in Research PlanetS supported by SNSF.

References

Appendix A: Change of coordinates in the general case

In this appendix I describe the change of coordinates corresponding to Eq. (12) in the general case. I introduce the angles $\begin{array}{ccc}& & {\mathit{\phi }}_{\mathit{i}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{i}}}(\frac{{\mathit{k}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}{{\mathit{q}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}{\mathit{\lambda }}_{\mathit{i}}\mathrm{+}\frac{{\mathit{k}}_{\mathit{i}\mathrm{+}\mathrm{2}\mathit{,i}\mathrm{+}\mathrm{1}}}{{\mathit{q}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{,i}\mathrm{+}\mathrm{2}}}{\mathit{\lambda }}_{\mathit{i}\mathrm{+}\mathrm{2}}\\ & & \mathrm{-}\left(\frac{{\mathit{k}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{,i}}}{{\mathit{q}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}\mathrm{+}\frac{{\mathit{k}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{,i}\mathrm{+}\mathrm{2}}}{{\mathit{q}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{,i}\mathrm{+}\mathrm{2}}}\right){\mathit{\lambda }}_{\mathit{i}\mathrm{+}\mathrm{1}}) \mathrm{(}\mathit{i}\mathrm{\le }\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{,}\\ & & {\mathit{\phi }}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{=}{\mathit{\lambda }}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{-}{\mathit{\lambda }}_{\mathit{n}}\mathit{,}\\ & & {\mathit{\phi }}_{\mathit{n}}\mathrm{=}\frac{{\mathit{k}}_{\mathit{n,n}\mathrm{-}\mathrm{1}}{\mathit{\lambda }}_{\mathit{n}}\mathrm{-}{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}{\mathit{\lambda }}_{\mathit{n}\mathrm{-}\mathrm{1}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathit{i}}\mathrm{=}{\mathit{\phi }}_{\mathit{n}}\mathrm{-}{\mathit{\varpi }}_{\mathit{i}}\mathit{,}\end{array}$(A.1)where the coefficients c_i are renormalizing factors given by Eq. (A.6). The actions canonically conjugated to σ_i are the AMD of the planets (D_i), while the actions canonically conjugated to φ_i are $\begin{array}{ccc}& & {\mathit{L}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{i}}{\mathit{c}}_{\mathit{i}}{\mathit{l}}_{\mathit{i,j}}{\mathrm{\Lambda }}_{\mathit{j}} \mathrm{(}\mathit{i}\mathrm{\le }\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{,}\\ & & \mathrm{\Gamma }\mathrm{=}\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\left({\mathrm{\Lambda }}_{\mathit{n}}\mathrm{+}\frac{{\mathit{k}}_{\mathit{n,n}\mathrm{-}\mathrm{1}}}{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\left({\mathrm{\Lambda }}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{+}\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}\mathrm{-}\mathrm{2}}}{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{2}\mathit{,n}\mathrm{-}\mathrm{1}}}\left(\mathit{...}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{2}\mathit{,}\mathrm{1}}}{{\mathit{k}}_{\mathrm{1}\mathit{,}\mathrm{2}}}{\mathrm{\Lambda }}_{\mathrm{1}}\right)\right)\right)\mathit{,}\\ & & \mathit{G}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{\Lambda }}_{\mathit{i}}\mathrm{-}{\mathit{D}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{G}}_{\mathit{i}}\mathit{,}\end{array}$(A.2)where ${\mathit{l}}_{\mathit{i,j}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{k}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}\sum _{\mathit{r}\mathrm{=}\mathit{j}}^{\mathit{i}}{\mathit{q}}_{\mathit{r,r}\mathrm{+}\mathrm{1}}\underset{\mathit{s}\mathrm{=}\mathit{r}\mathrm{+}\mathrm{1}}{\overset{\mathit{i}}{\mathrm{􏽙}}}\frac{{\mathit{k}}_{\mathit{s}\mathrm{+}\mathrm{1}\mathit{,s}}}{{\mathit{k}}_{\mathit{s}\mathrm{-}\mathrm{1}\mathit{,s}}} \mathrm{(}\mathit{j}\mathrm{\le }\mathit{i}\mathrm{\le }\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{.}$(A.3)It can be shown that the Hamiltonian does not depend on the angle φ_n, and the conjugated action G (total angular momentum of the system) is thus conserved (as expected). Since φ_n−1 is the only fast angle in the new system of coordinates, the averaging is done over this angle, and Γ (conjugate of φ_n−1) is thus constant in the averaged problem. The remaining 2(n−1) d.o.f. (once G and Γ are fixed), are thus defined by the angles σ_i (i ≤ n) and φ_i (i ≤ n−2), and their conjugated actions D_i (i ≤ n) and L_i (i ≤ n−2). The angles φ_i (i ≤ n−2) are the Laplace resonant angles between each group of three consecutive planets. The angles σ_n−1 and σ_n are the classical two-planet resonant angles between the two outermost planets, and any two-planet (not necessarily consecutive) resonant angle can be expressed as a combination of the angles σ_i (i ≤ n) and φ_i (i ≤ n−2).

In order to express the Hamiltonian in these new coordinates, one first needs to invert the change of coordinates (i.e., express Λ_i, λ_i, etc. as functions of L_i, φ_i, etc.). For the angles, the inverse transformation reads $\begin{array}{ccc}& & {\mathit{\varpi }}_{\mathit{i}}\mathrm{=}{\mathit{\phi }}_{\mathit{n}}\mathrm{-}{\mathit{\sigma }}_{\mathit{i}}\mathit{,}\\ & & {\mathit{\lambda }}_{\mathit{i}}\mathrm{=}{\mathit{\phi }}_{\mathit{n}}\mathrm{+}\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\underset{\mathit{j}\mathrm{=}\mathit{i}}{\overset{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}\frac{{\mathit{k}}_{\mathit{j}\mathrm{+}\mathrm{1}\mathit{,j}}}{{\mathit{k}}_{\mathit{j,j}\mathrm{+}\mathrm{1}}}{\mathit{\phi }}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{+}\sum _{\mathit{j}\mathrm{=}\mathit{i}}^{\mathit{n}\mathrm{-}\mathrm{2}}{\mathit{c}}_{\mathit{j}}{\mathit{l}}_{\mathit{j,i}}{\mathit{\phi }}_{\mathit{j}}\mathit{.}\end{array}$(A.4)Then, the Fourier expansion of the perturbative part of the averaged Hamiltonian exhibits angles of the form $\mathit{d}\mathrm{(}{\mathit{k}}_{\mathit{j,i}}{\mathit{\lambda }}_{\mathit{j}}\mathrm{-}{\mathit{k}}_{\mathit{i,j}}{\mathit{\lambda }}_{\mathit{i}}\mathrm{)}\mathrm{+}\mathit{p}{\mathit{\varpi }}_{\mathit{i}}\mathrm{+}\mathrm{(}\mathit{d}{\mathit{q}}_{\mathit{i,j}}\mathrm{-}\mathit{p}\mathrm{)}{\mathit{\varpi }}_{\mathit{j}}\mathit{.}$(A.5)Since all these angles should be expressed as combinations of σ_i (i ≤ n) and φ_i (i ≤ n−2), one has to make sure that only integers appear in the combinations. Indeed, if one of the coefficients is a fraction, then the corresponding angle will not be 2π-periodic. For instance, if one of the angles appearing in the Fourrier series is 3/2φ₁ + ..., then φ₁ = 0 is not equivalent to φ₁ = 2π. One should also make the coefficients as small as possible, to avoid introducing an unnecessary periodicity in the Hamiltonian (and thus an artificial duplication of fixed points). I thus choose the renormalizing factors c_i to obtain the smallest possible integers in the Fourier expansion of the perturbative part. ${\mathit{c}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{lc}{\mathrm{m}}_{\mathrm{r}\mathrm{\le }\mathrm{i}\mathit{,}\mathrm{r}\mathit{<}\mathrm{s}}\mathrm{\left(}{\mathrm{D}}_{\mathrm{i}\mathit{,}\mathrm{r}\mathit{,}\mathrm{s}}\mathrm{\right)}}{{\gcd }_{\mathit{r}\mathrm{\le }\mathit{i,r}\mathit{<}\mathit{s}}\mathrm{\left(}{\mathit{N}}_{\mathit{i,r,s}}\mathrm{\right)}} \mathrm{(}\mathit{i}\mathrm{\le }\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{,}$(A.6)where $\frac{{\mathit{N}}_{\mathit{i,r,s}}}{{\mathit{D}}_{\mathit{i,r,s}}}\mathrm{=}{\mathit{k}}_{\mathit{s,r}}{\mathit{l}}_{\mathit{i,s}}\mathrm{-}{\mathit{k}}_{\mathit{r,s}}{\mathit{l}}_{\mathit{i,r}} \mathrm{(}{\mathit{l}}_{\mathit{i,s}}\mathrm{=}\mathrm{0}\mathrm{for}\mathit{s}\mathit{>}\mathit{i}\mathrm{)}$(A.7)is the fraction reduced to its simplest form. The Hamiltonian is expanded in power series of the eccentricities, and truncated at a given degree d_max. Only the combinations for which the degree q_r,s = k_s,r−k_r,s ≤ d_max appear in the truncated Hamiltonian. Therefore, only these combinations should be considered in the computation of the coefficient c_i in Eq. (A.6).

For the actions, the inverse change of coordinates reads $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathit{i}}& \mathrm{=}& \frac{{\mathit{k}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}{{\mathit{c}}_{\mathit{i}}{\mathit{q}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}{\mathit{L}}_{\mathit{i}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{i}\mathrm{-}\mathrm{2}}}\frac{{\mathit{k}}_{\mathit{i,i}\mathrm{-}\mathrm{1}}}{{\mathit{q}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,i}}}{\mathit{L}}_{\mathit{i}\mathrm{-}\mathrm{2}}\\ & & \mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{i}\mathrm{-}\mathrm{1}}}\left(\frac{{\mathit{k}}_{\mathit{i,i}\mathrm{-}\mathrm{1}}}{{\mathit{q}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,i}}}\mathrm{+}\frac{{\mathit{k}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}{{\mathit{q}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}\right){\mathit{L}}_{\mathit{i}\mathrm{-}\mathrm{1}} \mathrm{(}\mathit{i}\mathrm{\le }\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{,}\\ {\mathrm{\Lambda }}_{\mathit{n}\mathrm{-}\mathrm{1}}& \mathrm{=}& \mathrm{\Gamma }\mathrm{-}\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\mathrm{(}\mathit{G}\mathrm{+}\mathit{D}\mathrm{)}\mathrm{+}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{n}\mathrm{-}\mathrm{3}}}\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}\mathrm{-}\mathrm{2}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{2}\mathit{,n}\mathrm{-}\mathrm{1}}}{\mathit{L}}_{\mathit{n}\mathrm{-}\mathrm{3}}\\ & & \mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{n}\mathrm{-}\mathrm{2}}}\left(\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}\mathrm{-}\mathrm{2}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{2}\mathit{,n}\mathrm{-}\mathrm{1}}}\mathrm{+}\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\right){\mathit{L}}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathit{,}\\ {\mathrm{\Lambda }}_{\mathit{n}}& \mathrm{=}& \frac{{\mathit{k}}_{\mathit{n,n}\mathrm{-}\mathrm{1}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\mathrm{(}\mathit{G}\mathrm{+}\mathit{D}\mathrm{)}\mathrm{-}\mathrm{\Gamma }\mathrm{+}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{n}\mathrm{-}\mathrm{2}}}\frac{{\mathit{k}}_{\mathit{n,n}\mathrm{-}\mathrm{1}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}{\mathit{L}}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathit{,}\end{array}$(A.8)where $\mathit{D}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{D}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{\Lambda }}_{\mathit{i}}\mathrm{-}\mathit{G}$(A.9)is the total deficit of angular momentum. I introduce $\mathit{\delta }\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}\mathrm{-}\mathit{G,}$(A.10)which is the nominal total deficit of angular momentum (at exact resonance). Since Λ_i,0 and G are constants, δ is also a conserved quantity. I additionally define $\mathit{ϵ}\mathrm{=}\mathit{D}\mathrm{-}\mathit{\delta }\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{\Delta }{\mathrm{\Lambda }}_{\mathit{i}}\mathit{,}$(A.11)which measures the distance of the system to the exact resonance, such that $\begin{array}{ccc}\mathrm{\Delta }{\mathrm{\Lambda }}_{\mathit{i}}& \mathrm{=}& \frac{{\mathit{k}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}{{\mathit{c}}_{\mathit{i}}{\mathit{q}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}\mathrm{\Delta }{\mathit{L}}_{\mathit{i}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{i}\mathrm{-}\mathrm{2}}}\frac{{\mathit{k}}_{\mathit{i,i}\mathrm{-}\mathrm{1}}}{{\mathit{q}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,i}}}\mathrm{\Delta }{\mathit{L}}_{\mathit{i}\mathrm{-}\mathrm{2}}\\ & & \mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{i}\mathrm{-}\mathrm{1}}}\left(\frac{{\mathit{k}}_{\mathit{i,i}\mathrm{-}\mathrm{1}}}{{\mathit{q}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,i}}}\mathrm{+}\frac{{\mathit{k}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}{{\mathit{q}}_{\mathit{i,i}\mathrm{+}\mathrm{1}}}\right)\mathrm{\Delta }{\mathit{L}}_{\mathit{i}\mathrm{-}\mathrm{1}} \mathrm{(}\mathit{i}\mathrm{\le }\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{,}\\ \mathrm{\Delta }{\mathrm{\Lambda }}_{\mathit{n}\mathrm{-}\mathrm{1}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\mathit{ϵ}\mathrm{+}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{n}\mathrm{-}\mathrm{3}}}\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}\mathrm{-}\mathrm{2}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{2}\mathit{,n}\mathrm{-}\mathrm{1}}}\mathrm{\Delta }{\mathit{L}}_{\mathit{n}\mathrm{-}\mathrm{3}}\\ & & \mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{n}\mathrm{-}\mathrm{2}}}\left(\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}\mathrm{-}\mathrm{2}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{2}\mathit{,n}\mathrm{-}\mathrm{1}}}\mathrm{+}\frac{{\mathit{k}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\right)\mathrm{\Delta }{\mathit{L}}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathit{,}\\ \mathrm{\Delta }{\mathrm{\Lambda }}_{\mathit{n}}& \mathrm{=}& \frac{{\mathit{k}}_{\mathit{n,n}\mathrm{-}\mathrm{1}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\mathit{ϵ}\mathrm{+}\frac{\mathrm{1}}{{\mathit{c}}_{\mathit{n}\mathrm{-}\mathrm{2}}}\frac{{\mathit{k}}_{\mathit{n,n}\mathrm{-}\mathrm{1}}}{{\mathit{q}}_{\mathit{n}\mathrm{-}\mathrm{1}\mathit{,n}}}\mathrm{\Delta }{\mathit{L}}_{\mathit{n}\mathrm{-}\mathrm{2}}\mathit{,}\end{array}$(A.12)with ΔL_i = L_i− ∑ _{j ≤ i}c_il_i,jΛ_j,0.

Appendix B: Coefficients of the Hamiltonian for a 3:4:6:8 resonant chain

In this appendix, I provide the expressions of the coefficients C_i,j that appear in the Hamiltonian of a 3:4:6:8 resonant chain at first order in eccentricities (see Eq. (13)). I follow the method described in Laskar & Robutel (1995). Each coefficient is a sum of two components, coming from the direct and indirect part of the perturbation $\begin{array}{ccc}{\mathit{C}}_{\mathit{i,j}}& \mathrm{=}& \mathrm{-}\mathrm{𝒢}\frac{{\mathit{m}}_{\mathit{i}}{\mathit{m}}_{\mathit{j}}}{{\mathit{a}}_{\mathit{j,}\mathrm{0}}}\sqrt{\frac{\mathrm{2}}{{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}}}\left({\mathit{C}}_{\mathit{i,j}}^{\mathrm{dir}}\mathrm{-}\frac{\mathrm{𝒢}{\mathit{m}}_{\mathrm{0}}}{\sqrt{{\mathit{\mu }}_{\mathit{i}}{\mathit{\mu }}_{\mathit{j}}{\mathit{\alpha }}_{\mathit{i,j}}}}{\mathit{C}}_{\mathit{i,j}}^{\mathrm{ind}}\right)\mathit{,}\\ {\mathit{C}}_{\mathit{j,i}}& \mathrm{=}& \mathrm{-}\mathrm{𝒢}\frac{{\mathit{m}}_{\mathit{i}}{\mathit{m}}_{\mathit{j}}}{{\mathit{a}}_{\mathit{j,}\mathrm{0}}}\sqrt{\frac{\mathrm{2}}{{\mathrm{\Lambda }}_{\mathit{j,}\mathrm{0}}}}\left({\mathit{C}}_{\mathit{j,i}}^{\mathrm{dir}}\mathrm{-}\frac{\mathrm{𝒢}{\mathit{m}}_{\mathrm{0}}}{\sqrt{{\mathit{\mu }}_{\mathit{i}}{\mathit{\mu }}_{\mathit{j}}{\mathit{\alpha }}_{\mathit{i,j}}}}{\mathit{C}}_{\mathit{j,i}}^{\mathrm{ind}}\right)\mathit{,}\end{array}$(B.1)where i<j, α_i,j = a_i,0/a_j,0, and $\begin{array}{ccc}{\mathit{C}}_{\mathit{i,j}}^{\mathrm{dir}}& \mathrm{=}& \frac{\mathrm{2}}{\mathrm{3}}{\mathit{\alpha }}_{\mathit{i,j}}^{-1}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{7}}{\mathrm{6}}{\mathit{\alpha }}_{\mathit{i,j}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{5}}{\mathrm{2}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{2}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{5}}{\mathrm{3}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{3}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\\ & \mathrm{\approx }& \mathrm{-}\mathrm{1.190494}\mathit{,}\\ {\mathit{C}}_{\mathit{j,i}}^{\mathrm{dir}}& \mathrm{=}& \mathrm{-}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{5}}{\mathrm{2}}{\mathit{\alpha }}_{\mathit{i,j}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{2}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\\ & \mathrm{\approx }& \mathrm{1.688311}\mathit{,}\\ {\mathit{C}}_{\mathit{i,j}}^{\mathrm{ind}}& \mathrm{=}& \mathrm{0}\mathit{,}\\ {\mathit{C}}_{\mathit{j,i}}^{\mathrm{ind}}& \mathrm{=}& \mathrm{1}\mathit{,}\end{array}$(B.2)for a 2/1 resonance between i and j (planets pairs 1, 3 and 2, 4), $\begin{array}{ccc}{\mathit{C}}_{\mathit{i,j}}^{\mathrm{dir}}& \mathrm{=}& \frac{\mathrm{4}}{\mathrm{5}}{\mathit{\alpha }}_{\mathit{i,j}}^{-2}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{6}}{\mathrm{5}}{\mathit{\alpha }}_{\mathit{i,j}}^{-1}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{31}}{\mathrm{30}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{11}}{\mathrm{10}}{\mathit{\alpha }}_{\mathit{i,j}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\\ & & \mathrm{+}\frac{\mathrm{23}}{\mathrm{15}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{2}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{16}}{\mathrm{5}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{3}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{32}}{\mathrm{15}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{4}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\\ & \mathrm{\approx }& \mathrm{-}\mathrm{2.025223}\mathit{,}\\ {\mathit{C}}_{\mathit{j,i}}^{\mathrm{dir}}& \mathrm{=}& \mathrm{-}{\mathit{\alpha }}_{\mathit{i,j}}^{-1}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\alpha }}_{\mathit{i,j}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{+}\mathrm{3}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{2}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{-}\mathrm{2}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{3}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\\ & \mathrm{\approx }& \mathrm{2.484005}\mathit{,}\\ {\mathit{C}}_{\mathit{i,j}}^{\mathrm{ind}}& \mathrm{=}& \mathrm{0}\mathit{,}\\ {\mathit{C}}_{\mathit{j,i}}^{\mathrm{ind}}& \mathrm{=}& \mathrm{0}\mathit{,}\end{array}$(B.3)for a 3/2 resonance (planets pair 2, 3), and $\begin{array}{ccc}{\mathit{C}}_{\mathit{i,j}}^{\mathrm{dir}}& \mathrm{=}& \frac{\mathrm{32}}{\mathrm{35}}{\mathit{\alpha }}_{\mathit{i,j}}^{-3}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{48}}{\mathrm{35}}{\mathit{\alpha }}_{\mathit{i,j}}^{-2}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{36}}{\mathrm{35}}{\mathit{\alpha }}_{\mathit{i,j}}^{-1}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{36}}{\mathrm{35}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\\ & & \mathrm{+}\frac{\mathrm{17}}{\mathrm{14}}{\mathit{\alpha }}_{\mathit{i,j}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{93}}{\mathrm{70}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{2}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{64}}{\mathrm{35}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{3}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\\ & & \mathrm{-}\frac{\mathrm{132}}{\mathrm{35}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{4}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{88}}{\mathrm{35}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{5}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\\ & \mathrm{\approx }& \mathrm{-}\mathrm{2.840432}\mathit{,}\\ {\mathit{C}}_{\mathit{j,i}}^{\mathrm{dir}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{16}}{\mathrm{15}}{\mathit{\alpha }}_{\mathit{i,j}}^{-2}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{8}}{\mathrm{5}}{\mathit{\alpha }}_{\mathit{i,j}}^{-1}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{19}}{\mathrm{15}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{13}}{\mathrm{10}}{\mathit{\alpha }}_{\mathit{i,j}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\\ & & \mathrm{-}\frac{\mathrm{53}}{\mathrm{30}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{2}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{18}}{\mathrm{5}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{3}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{12}}{\mathrm{5}}{\mathit{\alpha }}_{\mathit{i,j}}^{\mathrm{4}}{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\\ & \mathrm{\approx }& \mathrm{3.283257}\mathit{,}\\ {\mathit{C}}_{\mathit{i,j}}^{\mathrm{ind}}& \mathrm{=}& \mathrm{0}\mathit{,}\\ {\mathit{C}}_{\mathit{j,i}}^{\mathrm{ind}}& \mathrm{=}& \mathrm{0}\mathit{,}\end{array}$(B.4)for a 4/3 resonance (planets pairs 1, 2 and 3, 4). The coefficients ${\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ and ${\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ are the Laplace coefficients (e.g., Laskar & Robutel 1995), evaluated at α_i,j.

Appendix C: Fixed points at first order

In this appendix I describe in more detail how to determine the position of the fixed points of the averaged problem, at first order in eccentricities. In this case, the Hamiltonian takes the form $\mathrm{\mathscr{H}}\mathrm{=}{\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{\left(}\mathit{ϵ,}\mathrm{\Delta }{\mathit{L}}_{\mathit{i}}\mathrm{\right)}\mathrm{+}\sum _{\mathit{i,j}\ne \mathit{i}}{\mathit{C}}_{\mathit{i,j}}\sqrt{{\mathit{D}}_{\mathit{i}}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathit{i}}\mathrm{+}{{p}}_{\mathit{i,j}}\mathrm{·}{\phi }\mathrm{)}\mathit{,}$(C.1)where p_j,i = p_i,j are vectors of n−2 known integer coefficients (which depend on the considered resonances), and φ is the vector of φ_i (i ≤ n−2). The fixed points of the averaged problem are solutions of the following set of equations (see Eq. (18)) $\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }{\mathit{D}}_{\mathit{i}}}\mathrm{=}\frac{\mathit{\partial }{\mathrm{\mathscr{H}}}_{\mathrm{0}}}{\mathit{\partial ϵ}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}\sqrt{{\mathit{D}}_{\mathit{i}}}}\sum _{\mathit{j}\ne \mathit{i}}{\mathit{C}}_{\mathit{i,j}}\cos \mathrm{(}{\mathit{\sigma }}_{\mathit{i}}\mathrm{+}{{p}}_{\mathit{i,j}}\mathrm{·}{\phi }\mathrm{)}\\ \mathrm{0}& \mathrm{=}& \mathrm{-}\frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }{\mathit{\sigma }}_{\mathit{i}}}\mathrm{=}\sqrt{{\mathit{D}}_{\mathit{i}}}\sum _{\mathit{j}\ne \mathit{i}}{\mathit{C}}_{\mathit{i,j}}\sin \mathrm{(}{\mathit{\sigma }}_{\mathit{i}}\mathrm{+}{{p}}_{\mathit{i,j}}\mathrm{·}{\phi }\mathrm{)}\\ \mathrm{0}& \mathrm{=}& \frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }\mathrm{\Delta }{\mathit{L}}_{\mathit{i}}}\mathrm{=}\frac{\mathit{\partial }{\mathrm{\mathscr{H}}}_{\mathrm{0}}}{\mathit{\partial }\mathrm{\Delta }{\mathit{L}}_{\mathit{i}}}\\ \mathrm{0}& \mathrm{=}& \mathrm{-}\frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }{\phi }}\mathrm{=}\sum _{\mathit{i,j}\ne \mathit{i}}{{p}}_{\mathit{i,j}}{\mathit{C}}_{\mathit{i,j}}\sqrt{{\mathit{D}}_{\mathit{i}}}\sin \mathrm{(}{\mathit{\sigma }}_{\mathit{i}}\mathrm{+}{{p}}_{\mathit{i,j}}\mathrm{·}{\phi }\mathrm{)}\mathit{,}\end{array}$from which I deduce $\begin{array}{ccc}& & \sqrt{{\mathit{D}}_{\mathit{i}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{\mathit{\sigma }}_{\mathit{i}}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{\partial }{\mathrm{\mathscr{H}}}_{\mathrm{0}}}{\mathit{\partial ϵ}}\right)}^{-1}\sum _{\mathit{j}\ne \mathit{i}}{\mathit{C}}_{\mathit{i,j}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{p}}_{\mathit{i,j}}\mathrm{·}{\phi }}\mathit{,}\\ & & \mathrm{0}\mathrm{=}\sum _{\mathit{i,j}\ne \mathit{i,r}\ne \mathit{i}}{{p}}_{\mathit{i,j}}{\mathit{C}}_{\mathit{i,j}}{\mathit{C}}_{\mathit{i,r}}\sin \mathrm{(}\mathrm{(}{{p}}_{\mathit{i,j}}\mathrm{-}{{p}}_{\mathit{i,r}}\mathrm{)}\mathrm{·}{{\phi }}^{\mathrm{)}}\mathit{.}\end{array}$The parameter δ/ Γ only appears in these equations through the value of $\frac{\mathit{\partial }{\mathrm{\mathscr{H}}}_{\mathrm{0}}}{\mathit{\partial ϵ}}$. Therefore, at first order in eccentricities, all the angles, as well as the eccentricity ratios are independent of δ/ Γ. The value of this parameter only changes a factor common to all eccentricities (see Eq. (C.6)). Equation (C.7) provides a set of n−2 equations on the n−2 angles φ_i. For instance, for a 3:4:6:8 resonant chain such as Kepler-223, I obtain (see Eqs. (13) and (C.7)) $\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \mathrm{2}{\mathit{C}}_{\mathrm{1}\mathit{,}\mathrm{2}}{\mathit{C}}_{\mathrm{1}\mathit{,}\mathrm{3}}\sin \mathrm{\left(}\mathrm{2}{\mathit{\phi }}_{\mathrm{1}}\mathrm{\right)}\\ & & \mathrm{+}\mathrm{3}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{1}}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{3}}\sin \mathrm{\left(}\mathrm{3}{\mathit{\phi }}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}\mathrm{3}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{1}}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{4}}\sin \mathrm{(}\mathrm{3}{\mathit{\phi }}_{\mathrm{1}}\mathrm{+}{\mathit{\phi }}_{\mathrm{2}}\mathrm{)}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{1}}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{2}}\sin \mathrm{\left(}{\mathit{\phi }}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{1}}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{4}}\sin \mathrm{(}{\mathit{\phi }}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{)}\mathit{,}\\ \mathrm{0}& \mathrm{=}& {\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{1}}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{4}}\sin \mathrm{(}\mathrm{3}{\mathit{\phi }}_{\mathrm{1}}\mathrm{+}{\mathit{\phi }}_{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{3}}{\mathit{C}}_{\mathrm{2}\mathit{,}\mathrm{4}}\sin \mathrm{\left(}{\mathit{\phi }}_{\mathrm{2}}\mathrm{\right)}\\ & & \mathrm{+}\mathrm{2}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{1}}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{4}}\sin \mathrm{(}{\mathit{\phi }}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{2}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{2}}{\mathit{C}}_{\mathrm{3}\mathit{,}\mathrm{4}}\sin \mathrm{\left(}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{\right)}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{4}\mathit{,}\mathrm{2}}{\mathit{C}}_{\mathrm{4}\mathit{,}\mathrm{3}}\sin \mathrm{\left(}{\mathit{\phi }}_{\mathrm{2}}\mathrm{\right)}\mathit{.}\end{array}$(C.8)There are trivial solutions at 0 and π, but other (asymmetric) solutions might exist. The existence of asymmetric solutions is due to the influence of first-order resonances between non-consecutive pairs. In the case where only consecutive pairs are involved in resonances, Eq. (C.7) simplifies, and it can be shown that only symmetric solutions exist. The number of solutions of a highly non-linear set of equations such as Eq. (C.8) is not easily predicted. Moreover, among those solutions, some correspond to elliptical (stable) fixed points, and the others to hyperbolic (unstable) fixed points. In the case of Kepler-223, I solve for the position of fixed points numerically, and only consider the stable fixed points. I find six possible stable solutions.

Once a solution is found for φ_i (i ≤ n−2), the angles σ_i, and the eccentricity ratios can easily be deduced from Eq. (C.6) ${\mathit{e}}_{\mathit{i}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{\mathit{\sigma }}_{\mathit{i}}}\mathrm{\approx }{\mathit{e}}_{\mathrm{0}}\sum _{\mathit{j}\ne \mathit{i}}\frac{{\mathit{m}}_{\mathit{j}}}{{\mathit{m}}_{\mathrm{0}}}{\mathit{C}}_{\mathit{i,j}}^{\mathrm{\prime }}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{p}}_{\mathit{i,j}}\mathrm{·}{\phi }}\mathit{,}$(C.9)with $\begin{array}{ccc}{\mathit{C}}_{\mathit{i,j}}^{\mathrm{\prime }}& \mathrm{=}& \mathrm{-}\frac{{\mathit{a}}_{\mathit{n,}\mathrm{0}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{{\mathit{a}}_{\max \mathrm{\left(}\mathit{i,j}\mathrm{\right)}\mathit{,}\mathrm{0}}\sqrt{{\mathit{a}}_{\mathit{i,}\mathrm{0}}}}\left({\mathit{C}}_{\mathit{i,j}}^{\mathrm{dir}}\mathrm{-}\frac{{\mathit{C}}_{\mathit{i,j}}^{\mathrm{ind}}}{\sqrt{{\mathit{\alpha }}_{\mathit{i,j}}}}\right)\mathit{,}\\ {\mathit{e}}_{\mathrm{0}}& \mathrm{=}& \mathrm{-}{\mathit{n}}_{\mathit{n,}\mathrm{0}}{\left(\frac{\mathit{\partial }{\mathrm{\mathscr{H}}}_{\mathrm{0}}}{\mathit{\partial ϵ}}\right)}^{-1}\mathrm{·}\end{array}$(C.10)For a 3:4:6:8 resonant chain such as Kepler-223, I obtain $\begin{array}{ccc}\frac{{\mathit{e}}_{\mathrm{1}}}{{\mathit{e}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{-}\mathrm{3}{\mathit{\lambda }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{)}}& \mathrm{\approx }& \mathrm{6.2526}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\mathrm{+}\mathrm{1.9999}\frac{{\mathit{m}}_{\mathrm{3}}}{{\mathit{m}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathrm{2}{\mathit{\phi }}_{\mathrm{1}}}\mathit{,}\\ \frac{{\mathit{e}}_{\mathrm{2}}}{{\mathit{e}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{-}\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{2}}\mathrm{)}}& \mathrm{\approx }& \mathrm{-}\mathrm{6.5665}\frac{{\mathit{m}}_{\mathrm{1}}}{{\mathit{m}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{3}{\mathit{\phi }}_{\mathrm{1}}}\mathrm{+}\mathrm{3.0911}\frac{{\mathit{m}}_{\mathrm{3}}}{{\mathit{m}}_{\mathrm{0}}}\\ & & \mathrm{+}\mathrm{1.4999}\frac{{\mathit{m}}_{\mathrm{4}}}{{\mathit{m}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{\phi }}_{\mathrm{2}}}\mathit{,}\\ \frac{{\mathit{e}}_{\mathrm{3}}}{{\mathit{e}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathrm{4}{\mathit{\lambda }}_{\mathrm{4}}\mathrm{-}\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{3}}\mathrm{)}}& \mathrm{\approx }& \mathrm{-}\mathrm{0.5712}\frac{{\mathit{m}}_{\mathrm{1}}}{{\mathit{m}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}{\mathit{\phi }}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}\mathrm{)}}\mathrm{-}\mathrm{3.3120}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{2}{\mathit{\phi }}_{\mathrm{2}}}\\ & & \mathrm{+}\mathrm{3.1263}\frac{{\mathit{m}}_{\mathrm{4}}}{{\mathit{m}}_{\mathrm{0}}}\mathit{,}\\ \frac{{\mathit{e}}_{\mathrm{4}}}{{\mathit{e}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathrm{4}{\mathit{\lambda }}_{\mathrm{4}}\mathrm{-}\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{4}}\mathrm{)}}& \mathrm{\approx }& \mathrm{-}\mathrm{0.4284}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{\mathit{\phi }}_{\mathrm{2}}}\mathrm{-}\mathrm{3.2833}\frac{{\mathit{m}}_{\mathrm{3}}}{{\mathit{m}}_{\mathrm{0}}}\mathrm{·}\end{array}$(C.11)

All Tables

All Figures

Fig. 1 Location of the (six) fixed points for the Kepler-223 resonant chain (4 planets), as a function of the parameter δ/ Γ (see Sect. 2). I plot the eccentricities of the planets (left), the difference of longitudes of periastron between two successive planets, as well as the two-planet resonant angle σ4 (center), and the Laplace angles (right). Each fixed point is represented by the same color in all the plots. The gray bands (on the left and right columns) show the observed parameters of the planets (taken from Mills et al. 2016). For the eccentricities (left), the gray bands correspond to the 1σ uncertainties. For the Laplace angles, they correspond to the observed libration amplitudes (which are well constrained by TTVs, see Mills et al. 2016). The observed positions of the planets (and especially the Laplace angles), are consistent with a libration around the red fixed point (see also Fig. 2).

In the text

	Fig. 2 Comparison of the observed libration amplitudes of the Laplace angles in the Kepler-223 system (gray rectangle, values taken from Mills et al. 2016), with the positions of the (six) fixed points as determined from the analytical model. I use the same colors as in Fig. 1. The observations are consistent with a libration around the red fixed point.
In the text

Fig. 3 Example of a simulation that successfully reproduces the observed configuration of the Kepler-223 system. The scenario of capture is of type BCA (see Sect. 3.2). I plot the period ratio between each consecutive pair of planets (top), the planets’ eccentricities (middle), and the two Laplace angles (bottom). The two horizontal black lines in the bottom plot represent the equilibrium values of the Laplace angles expected from the analytical model (see Sect. 3.1).

In the text