Analytical model of multiplanetary resonant chains and constraints on migration scenarios
^{1} Observatoire de l’Université de Genève, 51 chemin des Maillettes, 1290 Sauverny, Switzerland
email: jeanbaptiste.delisle@unige.ch
^{2} ASD, IMCCE, Observatoire de Paris – PSL Research University, UPMC Univ. Paris 6, CNRS, 77 avenue DenfertRochereau, 75014 Paris, France
Received: 23 March 2017
Accepted: 22 June 2017
Resonant chains are groups of planets for which each pair is in resonance, with an orbital period ratio locked at a rational value (2/1, 3/2, etc.). Such chains naturally form as a result of convergent migration of the planets in the protoplanetary disk. In this article, I present an analytical model of resonant chains of any number of planets. Using this model, I show that a system captured in a resonant chain can librate around several possible equilibrium configurations. The probability of capture around each equilibrium depends on how the chain formed, and especially on the order in which the planets have been captured in the chain. Therefore, for an observed resonant chain, knowing around which equilibrium the chain is librating allows for constraints to be put on the formation and migration scenario of the system. I apply this reasoning to the four planets orbiting Kepler223 in a 3:4:6:8 resonant chain. I show that the system is observed around one of the six equilibria predicted by the analytical model. Using Nbody integrations, I show that the most favorable scenario to reproduce the observed configuration is to first capture the two intermediate planets, then the outermost, and finally the innermost.
Key words: celestial mechanics / planets and satellites: general
© ESO, 2017
1. Introduction
Mean motion resonances (MMR) between two planets are a natural outcome of the convergent migration of planets in a gasdisk (e.g., Weidenschilling & Davis 1985). The planets initially form farther away from each other, and planetdisk 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 multiplanetary 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 Kepler223 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 semianalytical model of threeplanet resonances taking into account only the interactions between consecutive planets in the chain, with a particular focus on the Kepler60 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, fourplanet (or more) resonances are not considered. Moreover, for some threeplanet resonances, the interactions between nonconsecutive planets cannot be neglected. For instance, in a 3:4:6 resonant chain, each planet is locked in a firstorder 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 Kepler223. 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 (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(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 twoplanet resonances (Delisle et al. 2012, 2014). The Hamiltonian of the system takes the form (Laskar 1991) (3)where is the gravitational constant, m_{i} is the mass of body i, a_{i} is the semimajor axis, r_{i} the position vector, and the canonically conjugated momentum of planet i (in astrocentric coordinates, see Laskar 1991). The first sum on the righthand side of Eq. (3) is the Keplerian part of the Hamiltonian (planetstar interactions), while the second sum is the perturbative part (planetplanet 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 twoplanet 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 Kepler223 (chain of 4 planets, 6 d.o.f.), and the problem is, in most cases, nonintegrable. 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 twoplanet 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 where is the angular momentum of planet i, β_{i} = m_{i}m_{0}/ (m_{0} + m_{i}), . At low eccentricities the deficit of angular momentum D_{i} is proportional to . The Hamiltonian (Eq. (3)) can be expressed using these actionangle coordinates (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 and ), and a Fourrier series of the angles, where the coefficients are functions of Λ_{i} and Λ_{j} (i.e., of the semimajor axes). For a system that is close to the resonance or resonant, the semimajor axes remain close to the nominal resonant values (Kepler’s third law) (7)I introduce (8)where (9)and expand the Keplerian part at degree 2, and the perturbative part at degree 0 in ΔΛ_{i}(10)where n_{i,0} is the nominal mean motion of planet i, such that (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 Kepler223). I introduce new canonically conjugated angles and actions as follows (see Eqs. (A.1) and (A.2)) (12)G is the total angular momentum, and is a conserved quantity; Its canonically conjugated angle (φ_{4}) does not appear in the Hamiltonian. The angle φ_{3} 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 φ_{3} in the Fourrier expansion of the Hamiltonian^{1}. The angle φ_{1} is the argument of the Laplace resonance between the three innermost planets. The angle φ_{2} 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 twoplanet 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 σ_{3} and σ_{4}. The arguments of the 3/2 resonance between planets 2 and 3 are σ_{2}−2φ_{2} = 3λ_{3}−2λ_{2}−ϖ_{2} and σ_{3}−2φ_{2} = 3λ_{3}−2λ_{2}−ϖ_{3}. The arguments of the 4/3 resonance between planets 1 and 2 are σ_{1}−2φ_{2}−3φ_{1} = 4λ_{2}−3λ_{1}−ϖ_{1} and σ_{2}−2φ_{2}−3φ_{1} = 4λ_{2}−3λ_{1}−ϖ_{2}. Arguments of resonances between nonconsecutive pairs can also be expressed in the same way. For instance, for the 2/1 resonance between planets 1 and 3, the arguments are σ_{1}−2φ_{2}−φ_{1} = 2λ_{3}−λ_{1}−ϖ_{1} and σ_{3}−2φ_{2}−φ_{1} = 2λ_{3}−λ_{1}−ϖ_{3}. 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 (), and obtain an expression of the form (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 semimajor 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 firstorder 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 (14)where (15)which measures the distance of the system to the exact resonance, where(16)is the total deficit of angular momentum, and (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 semimajor 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 Kepler223 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 twoplanet 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). 

Open with DEXTER 
For a given value of δ/ Γ, the fixed points are found by solving the following system of equations (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 nonzero real part.
3. Application to Kepler223
In this section, I apply my model to the four planets orbiting in a 3:4:6:8 resonant chain around Kepler223. 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 Kepler223, 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 φ_{1} and φ_{2} very well, while the eccentricities are only roughly determined. The twoplanet 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 φ_{1} and φ_{2}, 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, twoplanet 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 Kepler223 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. 

Open with DEXTER 
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 Kepler223 (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 Nbody 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^{6}, 3 × 10^{5}, 2 × 10^{5}, and 1.5 × 10^{5} 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 λ_{1} = 0, and the other planets’ semimajor axes and mean longitudes are randomly drawn. The mean longitudes are drawn from a uniform distribution between 0 and 2π. The semimajor axis ratio a_{2}/a_{1} is uniformly drawn in the range [1.23,1.3], a_{3}/a_{2} in the range [1.32,1.38], and a_{4}/a_{3} 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^{4} 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 firstorder 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 Kepler223. 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 twoplanet resonances, and then the two pairs join to form the fourplanet 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 Kepler223 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 Kepler223 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). 

Open with DEXTER 
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 twoplanet 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 Kepler223 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 Nbody 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 semimajor axes and initial mean longitudes of the planets. The statistics would probably slightly change with a different setup 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 Kepler223).

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 Kepler223, 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 Kepler223 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
 Cresswell, P., & Nelson, R. P. 2006, A&A, 450, 833 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Delisle, J.B., Laskar, J., & Correia, A. C. M. 2014, A&A, 566, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J.B., Correia, A. C. M., & Laskar, J. 2015, A&A, 579, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goździewski, K., Migaszewski, C., Panichi, F., & Szuszkiewicz, E. 2016, MNRAS, 455, L104 [NASA ADS] [CrossRef] [Google Scholar]
 Henrard, J. 1983, Icarus, 53, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1991, in Predictability, Stability, and Chaos in NBody Dynamical Systems, eds. S. Roeser, & U. Bastian, 93 [Google Scholar]
 Laskar, J. 2000, Phys. Rev. Lett., 84, 3240 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Laskar, J., & Robutel, P. 1995, Celes. Mech. Dyn. Astron., 62, 193 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Libert, A.S., & Tsiganis, K. 2011, Celes. Mech. Dyn. Astron., 111, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B. 2015, Int. J. Astrobiol., 14, 291 [CrossRef] [Google Scholar]
 Papaloizou, J. C. B. 2016, Celes. Mech. Dyn. Astron., 126, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Terquem, C. 2010, MNRAS, 405, 573 [NASA ADS] [Google Scholar]
 Steffen, J. H., Fabrycky, D. C., Agol, E., et al. 2013, MNRAS, 428, 1077 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J., & Davis, D. R. 1985, Icarus, 62, 16 [NASA ADS] [CrossRef] [Google Scholar]
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 (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 (A.2)where (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 twoplanet resonant angles between the two outermost planets, and any twoplanet (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 (A.4)Then, the Fourier expansion of the perturbative part of the averaged Hamiltonian exhibits angles of the form (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φ_{1} + ..., then φ_{1} = 0 is not equivalent to φ_{1} = 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. (A.6)where (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 (A.8)where (A.9)is the total deficit of angular momentum. I introduce (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 (A.11)which measures the distance of the system to the exact resonance, such that (A.12)with ΔL_{i} = L_{i}− ∑ _{j ≤ i}c_{i}l_{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 (B.1)where i<j, α_{i,j} = a_{i,0}/a_{j,0}, and (B.2)for a 2/1 resonance between i and j (planets pairs 1, 3 and 2, 4), (B.3)for a 3/2 resonance (planets pair 2, 3), and (B.4)for a 4/3 resonance (planets pairs 1, 2 and 3, 4). The coefficients and 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 (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)) from which I deduce The parameter δ/ Γ only appears in these equations through the value of . 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 Kepler223, I obtain (see Eqs. (13) and (C.7)) (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 firstorder resonances between nonconsecutive 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 nonlinear 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 Kepler223, 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) (C.9)with (C.10)For a 3:4:6:8 resonant chain such as Kepler223, I obtain (C.11)
All Tables
All Figures
Fig. 1 Location of the (six) fixed points for the Kepler223 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 twoplanet 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). 

Open with DEXTER  
In the text 
Fig. 2 Comparison of the observed libration amplitudes of the Laplace angles in the Kepler223 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. 

Open with DEXTER  
In the text 
Fig. 3 Example of a simulation that successfully reproduces the observed configuration of the Kepler223 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). 

Open with DEXTER  
In the text 