The onset of instability in resonant chains

There is evidence that most chains of mean motion resonances of type $k$:$k-1$ among exoplanets become unstable once the dissipative action from the gas is removed from the system, particularly for large $N$ (the number of planets) and $k$ (indicating how compact the chain is). We present a novel dynamical mechanism that can explain the origin of these instabilities and thus the dearth of resonant systems in the exoplanet sample. It relies on the emergence of secondary resonances between a fraction of the synodic frequency $2 \pi (1/P_1-1/P_2)$ and the libration frequencies in the mean motion resonance. These secondary resonances excite the amplitudes of libration of the mean motion resonances thus leading to an instability. We detail the emergence of these secondary resonances by carrying out an explicit perturbative scheme to second order in the planetary masses and isolating the harmonic terms that are associated with them. Focusing on the case of three planets in the 3:2 -- 3:2 mean motion resonance as an example, a simple but general analytical model of one of these resonances is obtained which describes the initial phase of the activation of one such secondary resonance. The dynamics of the excited system is also briefly described. This scheme shows how one can obtain analytical insight into the emergence of these resonances, and into the dynamics that they trigger. Finally, a generalisation of this dynamical mechanism is obtained for arbitrary $N$ and $k$. This leads to an explanation of previous numerical experiments on the stability of resonant chains, showing why the critical planetary mass allowed for stability decreases with increasing $N$ and $k$.


Introduction
The formation of planetary systems is one of the key questions of planetary science, however it remains to this date observationally poorly constrained. One can nonetheless contemplate fully formed planetary systems, of which we have examples galore thanks to exoplanet-hunting missions such as HARPS and Nasa's Kepler surveys, and consider what physical and dynamical mechanisms can produce them. Despite our limited knowledge on the real nature of exoplanetary systems, there are some clear trends in the exoplanet sample, which thus impose constraints on formation scenarios.
A very common type of exoplanet, which was unknown in our Solar System, is what we now call Super-Earths or Mini-Neptunes. These are planets having a mass of about 1 to 20M ⊕ (Earth's masses) and are found on relatively short orbital periods, of less than about 200 days. They are estimated to orbit a third to a half of all Sun-like stars (Mayor et al. 2011;Howard et al. 2012;Fressin et al. 2013;Petigura et al. 2013), and multiplanetary systems are not rare (of the order of a few hundred). Given that several of them host H/He gaseous atmospheres that cannot be explained by production of volatiles after the formation of the planet (e.g. Rogers 2015;Zeng et al. 2019), these planets are believed to form within the lifetime of their protoplanetary disc, and therefore interact dynamically with it. This type of interaction is called type-I migration: on the one hand the eccentricities of the orbits are damped by the disc, on the other hand (and on longer timescales) the orbit's size changes over time, usually shrinking, so that the planet is seen to migrate inward towards the star. At the inner edge of the disc, carved by the magnetic activity of the host star itself, another torque is activated which halts inward migration (a so-called planetary trap, e.g. Masset et al. 2006). In systems with multiple planets, when the inner planet has reached the inner edge of the disc, the second planet is still migrating inward, so the two planets are approaching each other. A preferred outcome of this convergent migration is the formation of compact chains of mean motion resonances, where the period ratio of neighbouring planets is close to a ratio of simple integers (Terquem and Papaloizou 2007;Cresswell and Nelson 2008;Morbidelli et al. 2008;Ogihara et al. 2015;Izidoro et al. 2017Izidoro et al. , 2019Pichierri et al. 2018).
Since these are transiting planets, their orbital period is known with extremely good precision; the period ratio distribution is therefore one of the best constrained distribution for exoplanets. Within the observed Super-Earth population, we do observe relatively long, coplanar resonant chains of planets, such as Trappist-1 (Gillon et al. 2016(Gillon et al. , 2017Luger et al. 2017) and Kepler-223 (Mills et al. 2016). However, an initially puzzling realisation is that the overall distribution of the period ratios is marked by systems that show little preference for near-integer period ratios, hosting planets with much wider orbital separations than those characterising resonant chains (e.g. Winn and Fabrycky 2015). This appears at first in striking contradiction with the type-I migration scenario for Super-Earths and Mini-Neptunes, which naturally produces resonant chains. This paradox is however only apparent, as pointed out by Izidoro et al. (2017). Their analysis showed that many orbital properties of observed Kepler systems (including orbital spacing and multiplicity distribution) are very well reproduced if a large fraction of resonant systems eventually become unstable in the Gy evo-Article number, page 1 of 21 arXiv:2004.07789v1 [astro-ph.EP] 16 Apr 2020 A&A proofs: manuscript no. PM2020InstabilityOfResonantChains_ArxivAA lution following the dissipation of the disc, with instability rates of ∼ 95%. The remaining stable systems naturally represent the observed resonant systems, such as Trappist-1, etc. In the original Izidoro et al. (2017) paper, only a limited fraction of resonant systems constructed via type-I migration went unstable within reasonable systems' lifetimes after the removal of the disc. However, in Izidoro et al. (2019) these high rates of postdisc phase instabilities needed to explain the Kepler data are actually recovered, especially in the simulations where the formed systems are more massive and more compact. They therefore conclude that the final number of planets in the chain, the compactness of the system and the planets' masses are crucial parameters that differentiate between systems that remain stable after disc removal (for total integration times of 50 -300 My) and system that suffer dynamical instabilities (collisions or ejections).
The results of Izidoro et al. (2017Izidoro et al. ( , 2019) motivate a careful dynamical analysis on the threshold of stability in mean motion resonant chains, and in this paper we focus on the dynamical mechanisms leading to the instability even in absence of external perturbations 1 . On this subject, an important numerical study was performed by Matsumoto et al. (2012). There, the authors studied numerically the stability of resonant multi-planetary systems for high-integer first-order mean motion resonances. They built the desired resonant configuration by simulating the convergent type-I migration phase in a protoplanetary disc of gas; then they slowly depleted the disc. They observed that there is a critical number of planets N crit above which the resonant systems go naturally unstable, with a crossing time comparable to that of non-resonant systems, and studied how this number changes with the planetary masses (m pl /M * , where M * is the stellar mass) and compactness of the chain (index k of the k:k − 1 resonance). More specifically, they demonstrated numerically that the critical number N crit which guarantees stability decreases with increasing compactness of the chain (increasing k) and increasing planetary mass m pl . The dynamical reason of the instability, however, was not discussed, nor the exact scaling law that links N crit , m pl and k.
The main goal of this paper is to investigate both analytically and numerically the dynamical mechanisms at the origin of the onset of instability in resonant chains, in order to explain the result of Matsumoto et al. (2012) and the large instability fraction of resonant chains observed in the (Izidoro et al. 2017(Izidoro et al. , 2019 simulations. More precisely we focus on the stability of resonant configurations with small amplitude of libration around a resonant equilibrium point. These configurations are the resonant states less susceptible to instabilities (Pichierri et al. 2018), and therefore represent a natural testing ground to assess the limits of stability of resonant chains. Because we intend to work analytically, and since the planetary Hamiltonian is not a continuous function of the number of planets N, it is convenient to rephrase the findings of Matsumoto et al. (2012) with the following equivalent statement: given the number N of planets and the compactness of the system (the resonant index k), there is a limit mass (m pl /M * ) crit for stability, which decreases with increasing N and k. Thus, in this paper we address the question of why resonant chains at an initial state of low amplitude of libra-tion become unstable if the planets are too massive, for different values of N and k. This work is the continuation of our previous paper Pichierri et al. (2018), in which we considered the stability of two deeply resonant planets as a function of the planetary mass.
In order to fix ideas, as in the case of two resonant planets, we consider systems of planets of the same mass, m i ≡ m pl , ∀i = 1, . . . , N. This is a useful simplification which allows one to grasp the main points having to work with only one parameter. We note also that individual Kepler systems seem to show a homogeneity in planetary masses (Weiss et al. 2018;Millholland et al. 2017), so this simplification does not constitute a major inconvenience. We will also consider coplanar orbits for simplicity. Indeed, if the chains that we intend to study are the result of capture in mean motion resonances during the disc phase, any significant mutual inclinations would have been damped out by the disc. Moreover, the few confirmed truly resonant systems (such as Trappist-1 or Kepler-223) show very small mutual inclinations. This suggests that resonant chains form in a relatively planar orbital configuration.
The remainder of this paper is organised as follows. In Section 2 we detail the setup for our numerical investigations, similar to the one used in Pichierri et al. (2018). In the (N + 1)-body simulations with N = 3 resonant planets, a new dynamical phenomenon is observed which was not present in the case N = 2, that triggers the instability of the resonant chains. In Section 3 we give a phenomenological description of this dynamical feature, and how it can explain the dependence of the limit mass for stability with the number N of the planets and the index k of the resonance, thus elucidating the numerical findings of Matsumoto et al. (2012); Izidoro et al. (2019). In Section 4, we give a detailed analytical description of this dynamical phenomenon in the exemplifying case N = 3, k = 3, and in Section 5 we generalise the analytical scheme to arbitrary N and k. Our conclusions are presented in Section 6. Finally, in Appendix A we summarise the main aspects of the numerical setup which allows to capture planets into mean motion resonance at different desired eccentricities.

Numerical maps of stability of resonant planets
In this section we describe the setup of our numerical investigation of resonant chains. Motivated by the results of Matsumoto et al. (2012), we investigate the stability of planets in chains of first order mean motion resonances in terms of the critical planetary mass m crit allowed for stability. Specifically, we want to understand why m crit decreases with the number of the planets N and the index of the resonance k along the chain.
The setup of our numerical experiments is the same as in our previous paper Pichierri et al. (2018) on two resonant planets, and we review it here briefly for ease of reading but refer to the first paper for the details. The underlying idea is similar to that of Matsumoto et al. (2012) (see also for example Ramos et al. 2017;Deck and Batygin 2015;Xu et al. 2018): planets are captured into mean motion resonance by running (N + 1)-body simulations with added dissipative forces that mimic disc-planet interactions of the type-I migration regime (relevant for Super-Earths and Mini-Neptunes). However, unlike Matsumoto et al. (2012), we do not attempt resonant capture experiments with different masses. The reason is that for relatively large planetary masses, close to the instability limit, the capture itself can become quite chaotic which may lead to large amplitudes of libration. Then, it becomes difficult to compare the long-term stability of these systems with large amplitude of libration with those with smaller masses that settle near the resonant equilibrium point. Instead, for a theoretical understanding of stability of a resonant chain as a function of planetary mass only, it is preferable to capture all the planets in resonance at low libration amplitudes at small masses and then, after gas removal, slowly increase the planetary masses until an instability is achieved. We stress that this growth in mass should not be interpreted as a physical process. It is just a numerical artifice to explore resonant dynamics as a function of the planetary mass and achieve an analytic understanding of the instability process.
Our numerical experiments to probe the stability of resonant planets thus consist of two phases. First, the desired number of planets is captured deeply in the desired resonant chain at low planetary mass, and we consider planets of the same mass for simplicity. We implement a planetary trap at the inner edge in order to ensure convergent migration which is needed for the planets to capture (e.g. Masset et al. 2006). Then the disc is slowly dissipated away, leaving the system in a state of small libration around a resonant equilibrium point (Pichierri et al. 2018, see also Appendix A) and only the pure conservative dynamics remains. In the second phase, the value of m pl is slowly increased at each time-step, maintaining the small amplitude of oscillation around the resonance, until un instability is reached (usually, the instability results in planetary collisions); again, this mass increase is purely fictitious and serves the only purpose to study for each value of m pl the stability of resonant configurations with the same level of excitation of the resonant degrees of freedom.
Built around the same numerical setup, we review below a few important aspects discussed in Pichierri et al. (2018) on the case of two resonant planets, as they will turn out to be relevant for the case of three and more planets as well. We then discuss the application of the numerical simulations on the stability of three resonant planets in Subsection 2.2.

Review on the case of two resonant planets
There are some points that should be revisited from our previous paper Pichierri et al. (2018) on the stability of mean motion resonances in two-planets systems. We summarise them below, and refer to Pichierri et al. (2018) for a full discussion.
The first point is that, when two planets are in a first order mean motion resonance k:k − 1, there are two frequencies associated with the libration of the system around the resonant equilibrium point, which we call resonant frequencies and indicate with ω res,i , i = 1, 2. These frequencies dictate the evolution of the system over long timescales, and are essentially associated to the evolution of the two resonant angles ψ i = kλ 2 − (k − 1)λ 1 − i , i = 1, 2, which indeed have slow variation under the assumption that the system is in the k:k − 1 mean motion resonance. Instead, on shorter timescales, the evolution is dominated by the non-resonant combination δλ 1,2 = λ 1 − λ 2 of the mean longitudes λ i , which is a fast-evolving angle; this angle is called synodic angle, and its frequency is called the synodic frequency ω syn . Sinceλ i = n i , n 1 /n 2 k/(k − 1) by the resonance condition, and n i is linked to the semi-major axis by Kepler's third law n i = GM * /a 3 , we have that the synodic frequency ω syn =δλ 1,2 = n 1 /k = GM * /a 3 /k. Thus, the synodic frequency is independent of the planetary mass and only depends on the nominal separation of the planets to the star. Instead, the resonant frequencies grow with the planetary mass m pl : for example, in a simple pendulum approximation of the mean motion resonant dynamics, the resonant frequencies are expected to grow as √ m pl (see Subsect. 4.3). Thus, for small enough plane-tary masses, the synodic frequency is much higher than the resonant frequencies, so that the two contributions happen on totally different timescales and are perfectly decoupled: then, the fast synodic degree of freedom can be averaged out and only the purely resonant evolution (the combination of both resonant frequencies) matters over a long time. However, at large enough planetary masses, the resonant frequencies might become comparable with the synodic frequency. When the ratio between the synodic frequency and resonant frequencies is close to an integer ratio, a secondary resonance is encountered: this means that the purely resonant degrees of freedom can now exchange energy with the synodic degree of freedom. Therefore, these secondary resonances between the synodic and resonant frequencies could in principle destabilise a resonant pair of planets. In Pichierri et al. (2018) we found that, in the case of two planets in first order mean motion resonance, these secondary resonances are active at such high planetary masses that the system actually becomes unstable at smaller values of m pl because of close encounters between the planets. Therefore, we concluded that these secondary resonances are not responsible for instability in a system of two resonant planets. The second point thus concerns the instability caused by close encounters in the case of resonant planets. This type of planetary instability is a well understood phenomenon, so that we can discriminate the orbital configurations that are stable with respect to close encounters (also called Hill-stable) and those that are not (e.g. Gladman 1993;Marchal and Bozis 1982;Petit et al. 2018). Following for example the approximation for initially circular and coplanar planets made in Gladman (1993), one has that (for a general, non-resonant system) if the orbital distance d = a 2 − a 1 satisfies 3r H,1,2 3.46r H,1,2 , (2.1) then the system is Hill-stable (see also Obertas et al. 2017). 2 Here r H,1,2 is the mutual Hill radius of the two planets, defined as r H,1,2 = a 1 + a 2 2 We found in Pichierri et al. (2018) that resonant planets are more stable with respect to close encounters than non-resonant ones, in the sense that, to suffer mutual scattering, the planets need to approach to each other significantly closer than d crit ; however we did find that close encounters destabilise the systems at lower planetary masses than the aforementioned secondary resonances would. Moreover, we found that the larger the amplitude of oscillation associated with the resonant motion around the resonant equilibrium point, the closer to d crit is the minimal physical distance for instability (the same remains true with respect to the more general criterion found e.g. in Marchal and Bozis 1982;Petit et al. 2018). It will be important to keep these two points in mind even in the case of three and more resonant planets, as they will be relevant for understanding their stability. We investigate the case N ≥ 3 below.

Numerical stability maps for three resonant planets
The first step is to perform numerical experiments as explained at the beginning of Section 2. We refer to Pichierri et al. (2018) for a more in-depth discussion on the setup for capture into mean motion resonance (including an analytical understanding of this process which is consistent with the Hamiltonian formalism and adiabatic theory), the subsequent phase of fictitious mass growth and how it can be understood analytically. There is only one small difference to be pointed out in the capture phase of our simulations. In Pichierri et al. (2018) we could obtain any desired value of e 2 (equivalently, e 1 ) by changing the value of the eccentricity damping timescale τ e . By setting a large value for τ e , large planetary eccentricities could be obtained (cfr. Equation (A.8)). Here, because the planets capture in resonance in sequence (first planet 1 and 2, then planet 3) if τ e is large, e 1 and e 2 can grow significantly before planet 3 enters in resonance. This can force large secular eccentricity oscillations of planet 3, which may preclude its resonant capture (see e.g. Batygin 2015 on criteria for resonant capture). We give the details of the setup for capture in Appendix A and describe a numerical recipe to overcome this difficulty, which poses no problem at all in the context of the second phase where we actually investigate the stability of the chains as a function of planetary mass.
2.2.1. Numerical stability maps for N = 3 and k = 3 We show in Figure 1a the result of four simulations of the second phase of our numerical experiments for the case N = 3 and the 3:2 -3:2 resonant chain, starting from different initial eccentricities. On the horizontal axis we report the (increasing) planetary mass, while on the vertical axis we show the evolution of the eccentricity. The simulations are stopped when an instability occurs (a collision in all cases). This plot is to be compared with the similar Figures 9 and 10 in Pichierri et al. (2018) for the case of N = 2 and the same resonance index k = 3, and uses the same scale on both axes to allow for an easier comparison. The approximate location of the observed instability for two planets in the same resonance is represented in Figure 1a by a dotted line. Comparing the cases N = 2 and N = 3, there are two important observations to make. The first is that the instabilities occur at lower masses in the case N = 3 than in the case N = 2. This is in agreement with the results of Matsumoto et al. (2012). This anticipated instability, in terms of planetary mass, is unlikely to be due to too-close encounters between pairs of planets as it was the case N = 2. This is because a resonant chain repeats the same orbital geometry between adjacent planets of a two-planet resonance of the same order. Thus, if the critical mass m crit corresponding to the instability in the case N = 3 is smaller, the minimal approach distance between each pair of neighbouring planets is necessarily larger in terms of mutual Hill radii than that causing an instability for N = 2. There is no apparent reason for which the threshold distance for destabilising two-body encounters should significantly change with the number N of planets in the system. So, the instability is likely to have a different cause. Upon close examination of the (N + 1)-body integrations shown in Figure 1a, one notices that an interesting phenomenon is evident. For m pl /M * < 1.28 × 10 −3 , the amplitude of oscillation of the eccentricity grows linearly with the planets' mass. This is due to the increasing amplitude of the fast-frequency term associated to the synodic terms (the same effect was present in the case of two planets); instead, the amplitude of libration associated to the purely resonant dynamics is conserved adiabatically. Then, at m pl /M * 1.28 × 10 −3 (the dashed vertical line in the figure) there is a sudden excitation of the amplitude of eccentricity oscillations. Upon close inspection of the numerical output with high temporal resolution, we realise that this excitation is now due to an actual increase of the amplitude of libration inside the resonance, as will be clear below. After the excitation at m pl /M * 1.28 × 10 −3 , the systems temporarily remain in resonance, albeit with an increased libration amplitude of the resonant angles; soon after, while the planetary mass is still increasing, the systems finally become unstable as the planets experience close encounters, eventually leading to collisions. This is observed in all simulations.
We have seen in Pichierri et al. (2018) (see also Sect. 2.1) that, with increasing amplitude of libration, the planets need to farther away from each other (in terms of mutual Hill radius) to be stable. On the other hand, the larger is the libration amplitude in the resonance, the closer the planets approach each other during their evolution. Thus, in order to remain stable, the planetary masses have to be smaller, so that the mutual Hill radius r H,1,2 shrinks and their minimal physical distance in terms of r H,1,2 remains large. In other words, we concluded that more excited resonant states become unstable at smaller planetary masses. So, our interpretation for the anticipated instability in the N = 3 case is the following: first some dynamical process excites the libration amplitude; then the planets become encounter-unstable because the threshold distance for instability exceeds the actual minimal distance of approach between planet pairs. Thus, below we will look for the dynamical mechanism increasing the libration amplitude. It should be noticed that if such mechanism exists, it would also preclude capture in the resonance at small libration amplitude for the corresponding planetary mass, which is what was observed by Matsumoto et al. (2012).

Numerical and analytical investigation of the phenomenon
In the previous subsection we have underlined the importance of the observed increase in the amplitude of libration around the equilibrium point in the (N + 1)-body simulations, and its relevance for triggering the instability of resonant chains. In the following we aim at better understanding the dynamical origin of this growth of libration amplitude. Our approach is to find a simplified N-planets Hamiltonian model which captures the main features of the dynamics that are observed in the complete (N + 1)-body integrations. This is because the complete model contains a virtually infinite number of harmonics, making it extremely hard to proceed analytically or to obtain any insights from the observed evolution. If we are able to observe the same phenomenon in a simplified problem it will be easier to isolate its origin. Thus, in the following we start from a Hamiltonian planetary model that has only a minimal number of terms (harmonics) and we progressively add more terms until we observe in the integration of the considered Hamiltonian the same phenomenon that we have seen in the full numerical integration. The Hamiltonian models are integrated numerically, while slowly increasing the mass of the planets at each integration time step in accordance with the (N + 1)-body simulations in Subsection 2.2.1. Only when the numerical integrations show very good agreement with the full (N + 1)-body integrations, will we consider the corresponding Hamiltonian as a good approximation to the full one and work directly with the former. Before we get into the technicalities of our investigation, we plan out our methodology below. The dotted curve indicates the limit of stability for a system of two planets deep in the 3:2 mean motion resonance (Pichierri et al. 2018): this shows that three resonant planets go unstable at lower masses than two resonant planets, in accord with Matsumoto et al. (2012). As explained in the main text, the anticipated instability is unlikely caused by close encounters, which were causing the instability in the the two-planet case. Indeed, in the case of three resonant planets a new dynamical phenomenon appears which is not observed in simulations of two planets: the system experiences an excitation in amplitude of oscillation before going unstable. This excitation, starting at m pl /M * 1.28 × 10 −3 (vertical dashed line) is more clearly visible in panel (b), where the result of one such numerical simulation is shown in light green. In panel (b) this simulation is also compared with the integration of two simplified models (dark green and orange lines), with the same initial conditions as the numerical simulation of the complete equations of motion. In both simplified models, only terms up to first order in the eccentricities are considered (cfr. Subsect. 4.2). The orange line represents the evolution of the averaged equations of motion where all non-resonant terms have been dropped: the evolution is initially qualitatively similar to the complete simulation, however no excitation is observed (cfr. Subsect. 4.3). The dark green line represents the evolution of a model with both resonant and synodic interaction terms for each planet pair: although only terms up to order one in the eccentricities have been considered, we see that the excitation at m pl /M * 1.28 × 10 −3 is well reproduced in this simplified system (cfr. Subsect. 4.4).
The first reasonable choice for the numerical integrations is to consider the averaged equations of motion, expanded to some order in the eccentricity. This corresponds to dropping all nonresonant harmonics from the planetary Hamiltonian (cfr. Subsect. 4.1) and only keeping resonant harmonics up to some order in e. This results in a system governed by a Hamiltonian H := H kepl + H res ; this approach is presented in Subsection 4.3. By doing so, one realises that these terms cannot alone be responsible for the increase in amplitude of libration observed in the (N + 1)-body integrations. This fact is anticipated in Figure 1b, where we plot with a dark orange line the evolution of the system governed byH over one of the full (N + 1)-body integration with the same initial conditions; we see that at first the two simulations are qualitatively equivalent (the slight differences emerge solely from the expansion up to first order in the eccentricities made in the truncated modelH), but the averaged model does not reproduce the excitation observed in the (N + 1)body simulation at the location of the dashed vertical line. Actually, we will show that such excitation in the purely averaged model is not possible at any value of the planetary mass m pl . This is the first main result of this section: the purely resonant system H = H kepl + H res with initial conditions at vanishing amplitude around a resonant equilibrium point is (Lyapunov) stable for all planetary masses. The next step is therefore to include additional non-resonant terms, which were naturally present in the full Hamiltonian that governs the evolution of the (N +1)-body integrations. Maintaining for simplicity the expansion to first order in the eccentricity (which should be valid at least when all eccentricities are small enough), we then add synodic terms. In the case of three planets, these include the harmonics λ 1 − λ 2 and λ 2 − λ 3 , which we add in an additional interaction Hamiltonian H syn . As we show in Subsection 4.4, the introduction of these terms is responsible for the same phenomenon observed in Figure 1a. This fact is anticipated in Figure 1b, where we plot with a darker colour the evolution of the system governed by H * := H kepl + H res + H syn over one of the full (N + 1)-body integration with the same initial conditions, and we see that there is good qualitative agreement between the two evolutions. We also investigate the possibility of adding only one of the two synodic terms, but show that both are needed to reproduce the phenomenon at similar planetary masses, which is a result that we will also explain analytically (cfr. Subsect. 4.4.2). In the light of this, we will use the evolution yielded by the simplified model H * = H kepl + H res + H syn as a guide to understand the relevant dynamics contained in the full (N + 1)-body integrations. At the same time, working with a controlled number of interaction terms allows us to proceed analytically (see Subsect. 4.4.1) and to understand what is the dy-namical mechanism that gives rise to the increase in amplitude of libration around the resonant equilibrium point. Carrying out the calculation explicitly in the specific case of N = 3 planets and for the 3:2 -3:2 chain, we show in Subsect. 4.4.2 that this is due to a set of secondary resonances between a fraction of the synodic frequency (which remains relatively constant with increasing m pl ) and specific combinations of the libration frequencies around the equilibrium point (which increase with m pl , as we will show). Considering relevant canonical action-angle variables centred at the equilibrium, such secondary resonances have the effect of exciting the action to values farther and farther away from the origin. This is the second main result of this section: the synodic contribution introduces terms of order O(m 2 pl ) which include secondary resonances between a fraction the synodic frequency and the resonant libration frequencies, which are responsible for the excitation of the system and eventually for its instability. In Subsect. 4.4.3 we build a model for the secondary resonance that is encountered in the specific case N = 3 and k = 3, but the method can be easily generalised to the other secondary resonances that can in principle be encountered. Finally, we proceed to generalise this result to more populated and/or more compact resonant chains in Section 5.

The origin of instability in resonant chains
In the next section, we will begin a careful analysis of the dynamics for three planets in a chain of mean motion resonances based on the insights elucidated above, which were lead by numerical integrations such as those of Figure 1. In particular, we aim at gaining a deep understanding of the process which causes the sudden excitation in the systems shown in Figure 1. As anticipated at the end of the last subsection, this process involves secondary resonances between some fraction of the synodic frequency ω syn = d dt (λ 1 − λ 2 ) and the resonant frequencies ω res,l associated with the libration of the system around the resonant equilibrium point. Before we delve into the dynamical details of these secondary resonances, let us delineate in a more general and practical sense why they are relevant for the problem of the stability of resonant chains of N planets.
The idea is that, normally, the synodic evolution (with characteristic frequency ω syn ) and the purely resonant evolution (with characteristic frequency ω res,l ω syn ) happen on such different timescales that there can be no interaction between them (as we already recalled in Subsection 2.1). However, a secondary resonance between them effectively allows energy to be transferred between the synodic and resonant degrees of freedom, and can ultimately cause an excitation of the latter which in turn makes the chain unstable to close encounters between the planets. Now, in the case of two planets, the resonant frequencies were too small compared to ω syn and grew too slowly with m pl , so that secondary resonances were active at such high planetary masses that the system was already unstable to close encounters (cfr. Subsect. 2.1). Note that for the same planetary mass m pl and for the same k, the libration frequencies for two and three resonant planets are roughly similar for similar eccentricities. However, the key point is that in the case N ≥ 3 there is a fraction of the synodic frequency which appears in the Hamiltonian (in terms at second order in the planetary masses). In the case of three planets, this fraction is ω syn /k where k as usual is the index of the resonance 3 . Thus, in the case of three planets, in order to multiplying ω syn decreases with increasing N and with k, the conclusion is that the regime of secondary resonances between synodic and resonant degrees of freedom is encountered at lower masses for increasing k and/or increasing N, and therefore the critical mass (m pl /M * ) crit allowed for stability decreases with N and with k. This mechanism gives a dynamical explanation to the numerical findings of Matsumoto et al. (2012); Izidoro et al. (2019). In the rest of this paper, we give a detailed analytical description of the dynamical emergence of these secondary resonances.

Hamiltonian model
In this section we describe the analytical tools used to investigate the emergence of secondary resonances between synodic and resonant degrees of freedom. We begin introducing the general planetary Hamiltonian and the customary notation in Subsection 4.1, and we then consider the relevant harmonic terms in the Hamiltonian that interest us in Subsection 4.2. Then, in Subsections 4.3 and 4.4 respectively we consider the averaged model H and the model H * which includes synodic terms. There, we give an analytical descriptions of the main dynamical features of the simulations shown in Figure 1.

Planetary Hamiltonian
We start with the Hamiltonian H of N planets of masses m i , i = 1, . . . , N orbiting a star of mass M * . We let u i be the inertial barycentric cartesian coordinate of each planet, andũ i = m iui the conjugated momentum. We write H in canonical heliocentric variables (p i , r i ), i = 1, . . . , N, defined from the inertial barycenat second order in m pl is the following. The Hamiltonian of three planets contains both the δλ 1,2 = λ 1 − λ 2 and δλ 2,3 = λ 2 − λ 3 harmonics. If both planet pairs are in the k:k − 1 mean motion resonance, one can write δλ 2,3 as (k − 1)δλ 1,2 /k plus some correction harmonic terms that only depend on the resonant angles (cfr. (4.14) with constant index k along the chain); this can be easily understood by noting thatδλ 2,3 should be comparable to (k − 1)δλ 1,2 /k in a k:k − 1 chain. Then, the two angles δλ 1,2 and (k − 1)δλ 1,2 /k get combined at second order in m pl which yields a harmonic containing δλ 1,2 /k plus purely resonant harmonics (cfr. (4.38)).
Pichierri, Morbidelli: The onset of instability in resonant chains tric canonical variables (u,ũ i ) as Poincaré 1892;Laskar 1990). Doing so, the Hamiltonian can be split as In other words, the Hamiltonian appears as a sum of two terms. One term is the sum of the Keplerian unperturbed Hamiltonians for each planet H kepl,i , describing the planet-star interactions. The other is the perturbing Hamiltonian H pert which describes all planet-planet interactions; H pert itself is split into direct terms, − 1≤i< j≤N Gm i m j / r i − r j , and indirect terms, 1≤i< j≤N p i · p j /M * , which come from having considered canonical heliocentric rather than barycentric variables.
where m pl is the typical mass of the planets) so it can be seen as a small perturbation to the integrable Keplerian Hamiltonian. For each planet, the canonical modified Delaunay variables can be introduced, which are action-angle variables for the reference Keplerian problems H kepl,i . We will consider only coplanar motion for the planets, so we only have two pairs of action-angle variables (Λ i , λ i ) and (Γ i , γ i ). Their definition in terms of the orbital elements is (e.g. Morbidelli 2002) As usual, for each planet a i is the semi-major axis, e i is the eccentricity, λ i is the mean longitude, i is the mean anomaly i is the longitude of the pericentre and µ i = m i M * /(M * + m i ) m i is the reduced mass; the index i = 1, . . . , N refers to the i-th planet, with planets ordered with increasing semi-major axis. We note that, as in Pichierri et al. (2018), the orbital elements are defined starting from heliocentric positions and barycentric velocities (4.1) (they are the so-called formal osculating elements, Morbidelli 2002).
In the modified Delaunay variables (4.3) the Keplerian part rewrites while no simple expression exists for H pert , which is usually expanded in Fourier series of the angles. In this expansion, there are only combinations of λ i and γ i which satisfy the d'Alembert characteristics, and only harmonic terms combining angles from two planets. We won't go into the details of how this expansion is performed in general, which can be found in many works (e.g. Laskar and Robutel 1995;Murray & Dermott 1999), and we will only concentrate on the specific terms that interest us below.

Rescaled Hamiltonian and new set of canonical variables
In order to make the calculations and algebraic expressions less cumbersome, we start by performing the following simplifications. These are clearly general and are carried out here for any number N of planets, but we will give specific examples to the case of 3 planets to fix ideas.
Firstly, since the instabilities for N ≥ 3 planets occur at much lower values of m pl /M * than for 2 planets, we approximate the reduced mass µ = m pl M * M * +m pl ∼ m pl and M * + m pl ∼ M * . Then, we recall that all the planets have the same mass m pl , and we intend later on to make use of the tools of perturbation theory to study the dynamics of the resonant chains. It is therefore convenient to write the Hamiltonian in the form of a sum of an integrable part which does not depend on the small parameter m pl , plus a small perturbation proportional to m pl . The natural choice is to rescale all the actions (Λ, Γ) of the modified Delaunay variables by the planetary mass m pl , which yields where for simplicity we have maintained the same notation as for the non-rescaled variables. In order to maintain the canonicity of the Hamiltonian, H itself must be rescaled by m pl . With this choice the reduced N-planets Hamiltonian takes the form (again, as for the canonical variables we do not change the notation for the rescaled Hamiltonian) where H kepl is independent of m pl , and the (rescaled) perturbation is of order O(m pl ): For a pair of neighbouring planets labelled by the indices i and i + 1 which are near a k (i) : (k (i) − 1) mean motion resonance, the perturbing resonant contribution to first order in the eccentricity takes the form where the (rescaled) coefficients α are where f ( j,i) res (α (i) res ) are functions of the Laplace coefficients b ( j) s (Murray & Dermott 1999), themselves (weakly) depending on the semi-major axis ratios (they include both direct and indirect terms; indirect terms only appear in the 2:1 mean motion resonance). Here as usual α (i) res =ā i /ā i+1 = (k (i) − 1)/k (i) 2/3 is the nominal semi-major axis ratio corresponding to the resonance location in the Keplerian approximation, so the Laplace coefficients are the same for each pair of planets in a resonant chain repeating the k:k − 1 commensurability. Moreover, we have evaluated the Λ 2 i+1 at denominator at its nominal Keplerian valueΛ i+1 (because the terms in (4.9) are already of order O(e), e.g. Batygin and Morbidelli 2013). Doing so, the coefficients α (i) j are effectively constants for a given chain and a given nominal orbital separation, and they represent the strengths of the resonances.
The other terms in the perturbing function H pert that are of interest to us are the synodic terms for each neighbouring planet pair. At lowest order in the eccentricities and lowest harmonic order in λ i − λ i+1 , they take the form where the coefficients C i for the rescaled Hamiltonian are (4.11) and have the same scaling inΛ i+1 as the coefficients in (4.9) but a different dependence on the Laplace coefficients b ( j) s (e.g. Murray & Dermott 1999; the term − α (i) res −1/2 comes from the indirect term of the perturbing function). Notice that (4.10) is of order 0 in eccentricity. The term O(e) cannot exist, because it would not satisfy the d'Alembert rules. So, (4.10) is all we have for the terms dependent on the difference of the mean longitudes of neighbouring planets λ i −λ i+1 , but independent of the resonant angles, in an expansion up to O(e) of the original Hamiltonian. At order 1 in eccentricity, there are also terms coupling resonant and synodic angles (e.g. the terms k , for an arbitrary integer j). Because, in what follows, they would behave like those in H syn in (4.10) but are O(e) smaller, we neglect them for simplicity. Notice also that in (4.10) we can limit ourselves to the lowest multiples of λ i −λ i+1 because we are looking for the slowest possible synodic frequency, as explained in Section 3.
In the following we will want to consider the case of N planets, each pair being near a k (i) : (k (i) − 1) mean motion resonance, and thus introduce the resonant angles as canonical coordinates. However, at the same time, we will want to make use of the the non-resonant synodic angles λ i − λ i+1 , so it is preferable that one of them, say λ 1 − λ 2 , be also one of the canonical variables. The natural choice is to use as canonical positions the resonant angles ψ (i) is the longitude of conjunction for the i-th pair) and the apsidal differences δγ i,i+1 = γ i − γ i+1 for i = 1, . . . , N − 1, then define δλ 1,2 = λ 1 − λ 2 and finally keep an angle which will not appear explicitly in the Hamiltonian, such as γ N . These linear changes of variables for the positions are easily extended to a canonical transformation (the transformation on the actions is linear, with matrix equal to the transpose of the inverse of the matrix defining the transformation on the angles). For example, in the case N = 3 the new angles will be while the new conjugated actions are 4 the canonicity of this transformation can easily be checked using the Poisson bracket criterion. This canonical change of variable has the advantage of being easily generalisable to any number N of planets and of having the specific angular momentum L appearing as an explicit constant of motion, since its conjugated angle γ N = −γ N never appears explicitly in the transformed Hamiltonian (all the other angles satisfy the d'Alembert rules, while this one does not so it cannot appear in the Hamiltonian function, even the non-averaged one). We remark that L is now the specific angular momentum because the actions have been rescaled by the planetary mass; this also entails that when integrating the system (4.6) with increasing m pl , L will always remain constant. Moreover, the action ∆λ 1,2 conjugated to the angle δλ 1,2 is simply a factor away from the action K used in Pichierri et al. (2018) (see also e.g. ; this action has been called the "spacing parameter" (Michtchenko et al. 2008) and is a constant in the averaged model where all non-resonant contributions to H pert are dropped, yielding information on the nominal locationΛ of the resonance at hand.
We integrate this Hamiltonian with a numerical integrator while slowly increasing the planetary mass at each time step as detailed above. We use again as an example the case of the 3:2 -3:2 chain starting with an initial planetary mass m pl /M * = 10 −5 and we choose as initial condition that of Figure 1b. The resulting evolution of the canonical actionsp is shown in dark green in Figure 2, panels (a) to (d) (the evolution of the eccentricity has been already presented in Figure 1b). We observe that the four resonant degrees of freedom are never unstable even up to masses significantly higher than the critical mass (m pl /M * ) crit 1.28 × 10 −3 which is found in the numerical (N +1)-body simulations with the same initial conditions ( Figure  1b, light green evolution in Fig. 2).
We can present an analytical explanation for this. As in Pichierri et al. (2018), we find the stable resonant equilibrium points forH(x; L, ∆λ 1,2 , m pl ) in the variablesx, while keeping L and ∆λ 1,2 constants and for different values of m pl , yieldinḡ x eq (m pl ) =x eq (m pl ; L, ∆λ 1,2 ). Notice that at these low eccentricities we are interested in symmetric linearly stable equilibria only, so the equilibrium valuesq eq of the angles are simply ψ (1) 1,eq = 0, ψ (2) 1,eq = 0, δγ 1,2,eq = π, δγ 2,3,eq = π, (4.18) and we only need to solve for the equilibrium actionsp eq = (Ψ (1) 1,eq , Ψ (1) 2,eq , ∆γ 1,2,eq , ∆γ 2,3,eq ). In Figure 2, panels (a) to (d), we superimposed the analytically calculated equilibrium points (dashed purple lines) and the numerically-obtained evolution, showing excellent agreement, which implies that the numerical solution stays on the stable equilibrium at all times. Then, we diagonalise the system around the equilibrium pointx eq ; since it is a stable equilibrium point, all eigenvalues are purely imaginary and the diagonalisation procedure yields a Hamiltonian of the form in cartesian coordinatesx = T (ξ, η), with T a transformation matrix. Using canonical polar coordinates (I l , φ l ) l=1,...,4 with which appears as the sum of four decoupled harmonic oscillators plus higher order terms. The resulting four frequencies ω l , l = 1, . . . , 4 are shown in Figure 2e as a function of the planetary mass m pl , and we notice right away that they all have the same sign. This means that at vanishing amplitudes of libration the Hamiltonian has an extremum at the equilibrium point (a maximum) so that we can use the Hamiltonian itself as a Lyapunov function to deduce that the equilibrium point is Lyapunov stable for all planetary masses. This means also that if the initial amplitude of libration around the equilibrium point is small, it has to remain small at all times.
Article number, page 9 of 21    Figure 1b. We show in dark green in panels (a) to (d) the evolution of the actionsp = (Ψ (1) 1 , Ψ (1) 2 , ∆γ 1,2 , ∆γ 2,3 ) as the planetary mass m pl is slowly increasing, and we match it to the calculated equilibriā p eq = (Ψ (1) 1,eq , Ψ (1) 2,eq , ∆γ 1,2,eq , ∆γ 2,3,eq ) (purple dashed line); we also add the corresponding (N + 1)-body integration with the same initial condition (light green). A legend for panels (a) to (d) is given in panel (a). We see that the system remains stable well after the value of m pl /M * 1.28 × 10 −3 corresponding to the onset of excitation in Figure 1b. Panels (e) and (f) contain the analytical explanation of the observed stability: we plot with coloured lines all the frequencies of the four degrees of freedom and we notice that they have the same sign, therefore the Hamiltonian has a maximum at the equilibrium point and for low amplitude of librations the system remains Lyapunov-stable even if the frequencies grow in absolute value. In panel (e) we used the same eccentricities that correspond to the initial conditions of panels (a) to (d), e 0.01; in panel (f) we used higher initial eccentricities, e 0.1. We note that the scaling law for ω l (m pl ) changes depending on the eccentricity (see black solid and dashed lines).
Since it will be useful later on, we also consider here how the libration frequencies grow with m pl . This is shown in Figure  2 panels (e) and (f). We find numerically that ω 1,2 ∝ m 2/3 pl at low eccentricities (e 0.01, panel (e)) while ω 1,2 ∝ m 1/2 pl at higher eccentricities (e 0.1, panel (f)). Notice that for a pendulumtype Hamiltonian like H pend (Σ, σ) = aΣ 2 − m pl b cos σ (4.21) the libration frequency would be ∝ m 1/2 pl , so it might be interesting to ponder analytically why at low eccentricities we get a different scaling. The reason is that with changing mass we also change the corresponding equilibrium point, which means that the parameters a and b in the pendulum-like Hamiltonian above also depend on m pl , and the real scaling would therefore be abm pl . The way the equilibrium points adjust to changes in m pl here is by following lines of constant specific angular momentum (see above, and Pichierri et al. 2018). Thus, with changing mass we also change the eccentricity of the corresponding equilibrium point, i.e. b in (4.21), as m 1/3 pl . We finally remark that Batygin (2015) estimates for two planets the (highest) libration frequency, at small amplitude of librations around the resonant equilibrium point and for a value of the angular momentum at which the separatrix first appears. He finds that this frequency scales with ((m 1 + m 2 )/M * ) 2/3 : since the appearance of the separatrix happens at small eccentricities, this is consistent with our findings.

The synodic contribution
In the previous subsection we have shown that the purely resonant system is Lyapunov-stable for all planetary masses. The next natural step is therefore to introduce non-resonant contribution of the disturbing function. To lowest order in e, we introduce the two synodic terms (4.10) for the inner and outer pairs that had been averaged out before, resulting in We integrate this Hamiltonian for the 3:2 -3:2 chain with the same numerical scheme described before and the same initial conditions as in the previous section. This gives the evolution of the actions displayed in dark green in Figure 3, which is matched against the (N+1)-body integration with the same initial datum (lighter green) and the locations of the equilibria forH calculated in the previous section for different planetary masses (purple dashed lines). The comparison for the eccentricity evolutions, instead of the canonical variables, has been already presented in Figure 1b. We notice two important aspects of these  1 , Ψ (2) 1 , ∆γ 1,2 , ∆γ 2,3 , ∆λ 1,2 ) for the system H * = H kepl + H res + H syn with the same initial condition as the (N + 1)-body integration of Figure 1b (the evolution of these variables in the (N + 1)-body integration is also shown here in light green for reference). The system follows on average the purple dashed lines, which correspond to the equilibria p eq for the systemH = H kepl + H res . In orange we show the evolution of the averaged variables p calculated through analytical averaging of the fast synodic frequencies, Equation (4.36). Note that, for m pl /M * < 1.28 × 10 −3 , p has very little oscillation around the p eq curve, compared to the p evolution. Instead, for m pl /M * > 1.28 × 10 −3 , the amplitude of oscillations of p and p around p eq are almost the same. This reveals that, while the initial oscillation of p is entirely due to the synodic terms and is effectively removed by passing to the p variables, it is then dominated by an increased amplitude of libration in the resonance. The evolution of the angular momentum L is also shown in panel (f), and it is of course a conserved quantity; panel (f) also contains the legend for all panels in this figure.
plots. The first is that, initially, for all variables the evolution described by H * follows on average that described byH (compare Fig. 2 with Fig. 3, and the dashed purple lines). This can be easily understood realising that H * contains fast, non-resonant angles, which, up to first order in the small parameter m pl , have simply been averaged out inH; therefore, as long as the O(m 2 pl ) contributions are unimportant, the only difference between the two evolutions are the short-periodic, O(m pl ) oscillations due to the δλ 1,2 synodic angle. We will actually study this effect analytically below. However, as soon as the O(m 2 pl ) remainder introduces important contributions to the dynamics, as in the case of the emergence of a secondary resonance, the dynamics described by the averagedH approximation is not valid anymore. This is indeed what we see in Figure 3, where a phenomenon similar to the one observed in the (N + 1)-body integrations appears, and at roughly the same value of m pl /M * 1.28 × 10 −3 , which was not found inH. Notice that such a secondary resonance cannot be caused by an interaction of the resonant degrees of freedomx only, as we have shown that these are stable for all values of m pl . Therefore, these secondary resonances must come from an interaction between some (combination) of the four resonant degrees of freedom and the synodic degree of freedom (∆λ 1,2 , δλ 1,2 ). In the following, we use the analytical tools of the Lie series perturbation theory in order to pinpoint the relevant secondary resonances that arise at order 2 in the planetary mass m pl . We carry out the calculation for the case of three resonant planets in any resonant chain order to get the general picture, but we will focus on the case of k (1) = k (2) = k, and k = 3 when needed.

Eliminating the O(m pl ) synodic term
In the previous section we dropped H syn out by averaging the Hamiltonian. But simple averaging or dropping of harmonics is not a rigorous procedure and, as we have seen, can alter the real dynamics. Averaging is just the first step of more complex, rigorous, perturbation approach, as we describe here.
The first step is to find a canonical transformation that, to first order in m pl , eliminates the synodic contribution m pl H syn from H * . This will introduce O(m 2 pl ) terms that we want to calculate explicitly, since they contain harmonics mixingq and δλ 1,2 , potentially associated to secondary resonances.
In order to eliminate m pl H syn at O(m pl ), we need to find a generating function χ syn that solves the homological equation where ∇ is the gradient with respect to the canonical variables and J is the standard symplectic matrix J = 0 −I n I n 0 . (4.26) Clearly χ syn will be of order m pl so we can write χ syn = m pl χ syn . From the expression for H syn , Equation (4.15c), we see that χ syn Article number, page 11 of 21 A&A proofs: manuscript no. PM2020InstabilityOfResonantChains_ArxivAA will have the form χ syn = m pl C 1 η 1 sin(δλ 12 ) (4.27) where the divisors η 1 and η 2 are immediately found in terms of the frequencies (4.16) of the unperturbed Keplerian Hamiltonian and the combination of angles appearing in the harmonics in H syn , yielding (4.28) These divisors are not vanishing nor small, since clearly η 1 = n 1 − n 2 , η 2 = n 2 − n 3 (remember that n i is the mean motion frequency of planet i and that the harmonics in H syn in the modified Delaunay variables were simply λ 1 −λ 2 and λ 2 −λ 3 ) and the planets are evidently far from the 1:1 resonance. Therefore equation (4.24) can indeed be solved.
This is however only valid until a point in which the O(m 2 pl ) contribution, which is still present in (4.34), has resonant effects (which happens at m pl /M * 1.28×10 −3 in Fig. 3). Indeed, as it is typical in perturbation theory, these terms are expected to contain higher-order harmonics which were not present in the original Hamiltonian H * = H kepl +H res +H syn : then, if these newly introduced O(m 2 pl ) Hamiltonian terms contains angles which, for certain values of m pl , have a vanishing or small enough frequency, they could not be eliminated by a further perturbative step because of the problem of small divisors, and may thus change the dynamics considerably. We therefore proceed to analyse these terms below.

The O(m 2 pl ) contribution
In this subsection, we look closely at the O(m 2 pl ) terms in (4.34). We are specifically interested in the harmonics that they contain, to find explicitly which combinations of angles q = (ψ (1) 1 , ψ (2) 1 , δγ 1,2 , δγ 2,3 ) and δλ 1,2 can give rise to secondary resonances at values of the planetary masses close to those where the increase in amplitude of libration is observed in the numerical integrations. Since the synodic frequency of δλ 1,2 is much higher than the libration frequencies characteristic of the anglesq , the most interesting harmonics are the ones where the lowest fraction of δλ 1,2 appears next to a combinations ofq . This is because these are the harmonic terms that will be linked to the secondary resonances that appear at lowest resonant libration frequencies, that is, by Figure 2 panels (e) and (f), at lowest planetary mass. The following calculation is clearly general, but to simplify matters we will quickly specialise to the case of a chain of three planets with both pairs in the same resonance, k (1) = k (2) = k, as well as to the reference case k = 3 for which the numerical integrations in Figure 1 were performed.
We start with the main term {{H kepl ,χ syn },χ syn } of order m 2 pl in (4.34). Since H kepl does not contain any angles, all secondary resonance contributions must come from combinations of the harmonics contained inχ syn . Recall that we de-finedχ syn containing both synodic terms with harmonics λ 1 − λ 2 and λ 2 − λ 3 , which we wrote in Equation (4.14) in terms of the new variables q. Therefore, the harmonics that are included in {{H kepl ,χ syn },χ syn } are combinations of these synodic harmonics; more specifically, they come from the products of their cosines 6 . Using the standard trigonometric identity cos(a) cos(b) = 1 2 (cos(a − b) + cos(a + b)), the resulting harmonics are 2δλ 1,2 , (4.37) so the harmonic with the lowest fraction of δλ 1,2 is (k (2) − k (1) + 1)δλ 1,2 − ψ (1) 1 + ψ (2) 1 + δγ 1,2 /k (2) . Specialising now to the case of a chain with the same resonance index k (1) = k (2) = k, this simply gives With the aid of an algebraic manipulator one can compute the full expression of {{H kepl ,χ syn },χ syn } and select the desired harmonic term (we used the software package Wolfram Mathematica), thus obtaining its coefficient (actually, one can see that this term emerges solely from the term ∝ {{1/Λ 2 2 ,χ syn },χ syn }). We avoid writing here the full expression, which is rather cumbersome, moreover as in (4.34) we evaluate it at the reference values of the actions so the term multiplying the cosine becomes a numerical coefficient, and we write this term as H scnd.res,kepl = const × m 2 pl cos (δλ 1,2 − ψ (1) 1 + ψ (2) 1 + δγ 1,2 )/k . Since we want to compare the frequency of δλ 1,2 /k to that of (−ψ (1) 1 + ψ (2) 1 + δγ 1,2 )/k, we need to consider the resonant HamiltonianH in the x variables, expand the "barred" variablesx around the equilibrium point characterised by the equilibrium actionsp eq and the equilibrium anglesq eq (Equation (4.18)) as in Subsect. 4.3, and then introduce the transformationx → (I l , φ l ) l=1,...,4 to the action-angle variables (I, φ), which transformsH into the sum of decoupled harmonic oscillators plus higher order terms, Equation (4.20). It is also useful to translate the value of ∆λ 1,2 around its initial reference value ∆λ 1,2 introducing ∆λ 1,2 = ∆λ 1,2 + δ∆λ 1,2 which is clearly a canonical transformation. Therefore, we write H scnd.res,kepl in terms of the variables (I, δ∆λ 1,2 , φ, δλ 1,2 ). The Hamiltonian H scnd.res,kepl will now contain harmonic terms of type sin cos (δλ 1,2 /k + h· φ), (4.40) where h· φ is an integer combination with coefficients h 1 , . . . , h 4 ∈ Z of the angles φ 1 , . . . , φ 4 , which we can calculate explicitly. Therefore, whenever d dt δλ 1,2 /k = − d dt (h· φ) a 6 This can be easily understood noting that if χ = sin(q 1 ) + sin(q 2 ) and f depends only on the actions p i then {{ f, χ}, χ} = ∂ 2 f ∂p 2 1 cos 2 q 1 + 2 ∂ 2 f ∂p 1 ∂p 2 cos q 1 cos q 2 + ∂ 2 f ∂p 2 secondary resonance is crossed. We can rewrite this expression asδλ 1,2 /k + h· ω = 0. Since the Hamiltonian has d'Alembert characteristics in each pair (I l , φ l ), and the values of the actions I are initially (that is, before their excitation) small, the strongest secondary resonances will come from lowest integer combinations h· φ, that is, where most h l are zero. We also note that sincė δλ 1,2 > 0 and the frequencies ω l are all negative, a secondary resonance term can only appear when h· ω < 0, which together with the requirement that |h| be small is tantamount to requiring that all non-zero integers h l are positive. Since we calculated ω(m pl ) in Subsect. 4.3, we can calculate for each h the relative frequency δ λ 1,2 /k + h· ω (m pl ) as a function of m pl , and check if any of these vanish for some value of m pl , which corresponds to crossing a secondary resonance.
We carried out the calculation with the aid of the Mathematica software in the reference case k = 3 and a 1 0.1, which corresponds to the evolution shown in Figure 1b (and also Figures 2 and 3). We found that H scnd.res,kepl contains, among many others, the following terms The nature of these harmonics is clearly general, while the numerical coefficients are specific to the reference case k = 3 and a 1 0.1 mentioned above. We then calculated for each of the harmonics in (4.41) their frequency δ λ 1,2 /k + h· ω (m pl ) as a function of the mass. The results are presented in Figure 4.
We immediately remark that in the case of the harmonic δλ 1,2 /3 + 2φ 2 , the crossing of the secondary resonance happens precisely at the value of planetary mass m pl /M * 1.28 × 10 −3 where the numerical integrations showed the increase in amplitude of libration (see Figures 1, 2, 3). This is evidence that this phenomenon was indeed caused by the crossing of this secondary resonance.
Before we continue with an analytical description of the dynamics caused by this resonance, we should however go back and discuss a few technical details.
Firstly, if we had used in H syn only one synodic term, not all of the the harmonics in (4.37) would appear 7 . In particular, the harmonic δλ 1,2 /3 + 2φ 2 would not appear, so that the observed dynamical effects linked to the crossing of secondary resonances at m pl /M * 1.28 × 10 −3 are not expected. Indeed, we (1 + cos(2q 1 )).
Clearly we do not obtain the needed (k (2) − k (1) + 1)δλ 1,2 − ψ (1) 1 + ψ (2) 1 + δγ 1,2 /k (2) in (4.37) neither when q 1 = λ 1 − λ 2 = δλ 1,2 nor when q 1 = λ 2 − λ 3 = 1 k (2) (k (1) − 1)δλ 1,2 + ψ (1) 1 − ψ (2) 1 − δγ 1,2 .  +ω 1 +ω 4 Fig. 4: Frequencies of the angles δλ 1,2 /3 + h· φ as a function of the planetary mass in the case of the 3:2 -3:2 mean motion resonance chain with a 1 0.1 (the situation depicted in Figure 1b). Notice that the synodic frequencyδλ 1,2 varies only slightly due to the change in the equilibrium pointx eq for the averaged HamiltonianH, which is followed by the full system H until the second order effects become significant (cfr. Equation (4.34)). The main change comes from the resonant frequencies ω, whose dependence on the planetary mass is depicted in Figure 2e. The result is that the frequencies δ λ 1,2 /k + h· ω (m pl ) vanish within a small range of values of the planetary mass m pl , meaning that a capture into a secondary resonance becomes possible. By comparing with Figure 1b, we see that δλ 1,2 /3+2φ 2 has vanishing frequency at the same value of m pl /M * 1.28 × 10 −3 at which the excitation of the system occurs.  Figure 1b (lightest green) and two numerical simulations with the same initial conditions where only one of the two synodic terms λ 1 − λ 2 and λ 2 − λ 3 appears (two darker shades of green). These two semi-synodic simulations initially appear identical, which is easily understood from the fact that they follow in average the evolution ofH = H kepl + H res which is the same for the two (cfr. Equation (4.32)). The important point is that in both cases, when only one synodic angle is considered, the system is not excited at value of m pl /M * 1.28 × 10 −3 , where it is excited in the (N + 1)-body simulation as well as in the numerical simulation which includes both synodic terms, see Figure 1b. This shows that both synodic terms must be included in order to have a good quantitative agreement with the (N + 1)-body simulations.
With these clarifications, we can proceed with the model of the secondary resonance linked to the angle δλ 1,2 /3+2φ 2 , which, as we discussed above, has vanishing frequency exactly at the value of m pl when the increase in the amplitude of libration is observed in Figure 1b. This realisation is further supported by Figure 6. There, we plot the evolution of the actions I l , l = 1, . . . , 4 along the simulation, with the planetary mass m pl on the horizontal axis. We see that initially only one action is excited, namely I 2 , and after that the nonlinearities inherent in the system cause an exchange of energy between the degrees of freedom. This also suggests that the model that we are about to construct, which is valid only for small I's, breaks down whenever one of the actions is excited. This however presents no impediment in the description of the first phase, when the secondary resonance is encountered. One question that we wish to answer for example is whether or not there is or can be a capture in this secondary resonance or rather a jump across resonance. The integrable, low order model that we construct below can indeed answer this question.

4.4.3.
Model of the secondary resonance for δλ 1,2 /3 + 2φ 2 In the following we detail how we can construct a model for the resonance associated with the angle δλ 1,2 /3 + 2φ 2 since, as we saw before, it is the one that causes the observed increase in amplitude of libration. A similar approach can be implemented for the other resonances in (4.41).
We start by performing a canonical transformation which selects δλ 1,2 /3 + 2φ 2 as an angle. Notice that, since φ 2 appears with a coefficient 2 and so √ 2I 2 appears as a power two in (4.41b), we have a secondary resonance of order 2; hence it is useful to define the resonant angle θ as 2θ = δλ 1,2 /3 + 2φ 2 in order to maintain the d'Alembert characteristics so that the Hamiltonian will not be singular at the origin. The resulting transformation is Θ = I 2 , θ = δλ 1,2 /6 + φ 2 , I * r = I r , r = 1, 3, 4, φ * r = φ r , r = 1, 3, 4, (4.43) δ∆λ * 1,2 = δ∆λ 1,2 − I 2 /6 δλ * 1,2 = δλ 1,2 , whose canonicity follows immediately from the preservation of the Poisson brackets. We can already notice that Θ = I 2 appears as the conjugated action to the angle θ associated to the secondary resonance: this explains why in Figure 6 it is I 2 which is initially excited. The other variables do not feel the resonance, except δ∆λ 1,2 which must change according to the change in I 2 in order to maintain δ∆λ * 1,2 constant; however since I 2 gets divided by 6 this change is minute, but nevertheless clearly visible in Figure 3e. The pair (Θ, θ) is the pair of resonant variables for this specific secondary resonance, while the others will have a faster evolution, which can be "averaged" away, in order to yield a 1-d.o.f. system that we write H (φ * ,δλ * 1,2 ) (Θ, θ; I * , δ∆λ * 1,2 ). The notation • (φ * ,δλ * 1,2 ) means that we eliminated perturbatively to lowest order the non-secondary-resonant contributions from the angles (φ * , δλ * 1,2 ), and we stressed that the variables (I * , δ∆λ * 1,2 ) will only play the role of parameters for H (φ * ,δλ * 1,2 ) . Ultimately, the functional form of H (φ * ,δλ * 1,2 ) (Θ, θ) will be that of a Andoyer Hamiltonian, that is the coefficient c will be of order m 2 pl , while δ and β will be of order m pl . Since the system is initially close to the resonant equilibrium point, Θ is small and we can drop the O(Θ 3 ) terms. However, as we will see below, the parameter β (the second derivative at Θ = 0) plays a crucial role in determining if there can be capture into the secondary resonance or not, so we must keep track of all O(Θ 2 ) terms of the θ-independent part, that is, the first two terms in (4.44). The main contribution to the θ-independent part comes from theH term (the O(m pl ) term in (4.34)), while c √ 2Θ 2 cos(2θ) comes from (4.41b) and is O(m 2 pl ). Concerning the first part deriving fromH, we should stress that even if δ∆λ 1,2 appeared as a constant of motion when this Hamiltonian was treated alone, when the O(m 2 pl ) is taken into account the transformation (4.43) transforms δ∆λ 1,2 into δ∆λ * 1,2 + Θ/6, where δ∆λ * 1,2 is the new constant of motion. Therefore we must keep δ∆λ 1,2 as a variable inH and apply (4.43) to it.
With these considerations in mind we can obtain analytical insights on the dynamics, at least as long as the actions remain small (recall that before any secondary resonance is encountered, the system is very close to the equilibrium point at vanishing amplitude of libration). It is interesting for example to explore analytically if there can be a capture in this secondary resonance or not. Capture into resonance is possible (but not guaranteed) only ifδβ < 0. Intuitively, this is because near the origin one haṡ θ δ+βΘ and the resonance condition imposes that (on average) this quantity vanishes; therefore, at the centre of the resonance Θ = −δ/β, which only makes sense when β and δ have opposite signs. Thus, since β remains relatively constant (see below),  (4.20)). Then, I 2 increases steadily, symptom of an interaction with a secondary resonance that involves Θ = I 2 as a resonant action; this is confirmed by the canonical change of coordinates (4.43). Soon after I 2 is large enough, the degrees of freedom start interacting and exchanging energy, due to the non-linear effects.
only whenδβ is negative does the resonance centre appear from the origin and move at higher values of Θ, while ifδβ is positive the resonance centre approaches the origin from far away, the orbit is invested by a separatrix and then the resonance disappears leaving behind an excited orbit. We already know from Figure 4 that, as the planetary mass increases,θ goes from positive values to negative values, that is, thatδ < 0: this means a capture into this secondary resonance is possible only if β > 0.
To obtain the sign of β in (4.44) we need to compute its value explicitly. We do this in steps as follows. First, we fix a value of m pl right before the observed increase of amplitude of libration, m pl /M * 1.28 × 10 −3 , we calculate the equilibrium pointx eq =x eq (m pl ) and we apply the canonical diagonalisation procedure as explained in Subsect. 4.3. This yields four pairs of cartesian canonical variables (ξ, η) which replace thex: x = T (ξ, η), with T the diagonalasing matrix. Second, as in Subsect. 4.3, we introduce canonical polar coordinates (I l , φ l ) l=1,...,4 by ξ l = √ 2I l cos φ l , η l = √ 2I l sin φ l . The HamiltonianH will then depend on the variables (I 1 , . . . , I 4 , δ∆λ 1,2 , φ 1 , . . . , φ 4 , δλ 1,2 ). Third, we writeH in the variables (4.43);H contains a term in Θ 2 independent of the angles, but its coefficient is not β. To obtain the value of β we need to perform a fourth step, and calculate H (φ * ,δλ * 1,2 ) , that is the perturbative elimination inH of all the non-secondary-resonant contributions from the angles (φ * , δλ * 1,2 ), up to order 2 in Θ. This is because, as detailed below, the elimination of these harmonics can generate terms in Θ 2 independent of the angles, that need to be added to the original term to obtain β. To this end, we takeH(I 1 , . . . , I 4 , δ∆λ 1,2 , φ 1 , . . . , φ 4 , δλ 1,2 ) and expand it to order 2 with respect to the actions. Since these terms satisfy the d'Alembert characteristics in (I, φ), we only obtain terms like where α ∈ N 4 0 , |α| = α 1 + · · · + α 4 is restricted to |α|/2 ≤ 2, and c α (δ∆λ 1,2 ) is a coefficient which depends on δ∆λ 1,2 only. These coefficients are expanded around δ∆λ 1,2 = 0 to an optimal order which can be obtained in the following manner. Note that, from δ∆λ 1,2 = δ∆λ * 1,2 + Θ/6 (Equation (4.43)) each term of order d in δ∆λ 1,2 contributes a term of order d in Θ, so for each term of order |α|/2 in I we must obtain c α (δ∆λ 1,2 ) only up to order 2 − |α|/2 in δ∆λ 1,2 (where • is the floor function) to achieve the desired second order with respect to all the actions. We can then organise all terms with respect to the order of expansion in I and δ∆λ 1,2 , and write for each addendH j s/2 = O(I s/2 ) × O(δ∆λ 1,2 j ). To get a sense of what these terms look like, we write the terms only up to order s = 2 in √ I: the subsequent terms of higher order in √ I follow this structure but the possible combinations of the angles get substantially more numerous and we avoid writing them all here in the interest of brevity. Among them, there are of course also the terms ∝ I 2 l appearing without angles, as well as the term ∝ δ∆λ 1,2 2 , which contribute directly to the Θ 2 term in (4.44). We note that the first termsH 0 1/2 ∝ √ 2I l cos(φ l ) (corresponding to |α| = 1 and constant in δ∆λ 1,2 ) are actually zero by definition of equilibrium point calculated at the reference value δ∆λ 1,2 = 0. Also, the coefficients inH 0 1 in front of the I l 's are just the frequencies ω l , since these are the ones calculated in (4.20). Therefore, in this case the role of the integrable part of the Hamiltonian for a perturbation theory step is naturally played byH 0 1 = 4 l=1 ω l I l . To understand how the perturbative elimination of the nonresonant harmonics involving (φ * , δλ * 1,2 ) can generate terms in Θ 2 independent of the angles, consider that if f n = O(I n ) and 1 , χ m } +H j n,pert = 0, we must naturally use a χ n = O(I n ). This introduces new terms {{H 0 1 , χ n }, χ n } and {H j n,pert , χ n } which are O(I 2n−1 ) (given thatH 0 1 is O(I)). Thus, terms of order 2 can be generated for example if n = 3/2. We actually need to calculate explicitly only those terms that yield a Θ 2 independent of the angles, as the others would be eliminated further. Such terms derive fromH 0 3/2 (which governs the nonlinear interactions between the four resonant degrees of freedom around the equilibrium point, which Figure 6 proves to be strong) andH 1 1/2 (which describes the fact that the equilibrium pointp eq ofH shifts as δ∆λ 1,2 changes under the effects of the O(m 2 pl ) terms). We implemented this procedure with the aid of the algebraic manipulator Mathematica. In our reference case k = 3 and a 1 0.1 at a mass m pl /M * 1.28 × 10 −3 right before the development of the excitation of the resonant degrees of freedom ( Figure 6) this yields The fact that δ is positive and small is consistent with the fact that we put ourselves right before the development of the excitation (cfr. Figure 4). The fact that β ∼ 100 is positive yields an analytical confirmation that there can be capture into this secondary resonance. After we have obtained H (φ * ,δλ * 1,2 ) (cfr. Equation (4.47)), we can easily complete the determination of the model (4.44) for this secondary resonance. To do this, with the help of the algebraic manipulator Mathematica we use the canonical transformation (4.43) applied to the term (4.41b) which contains the resonant harmonic 2θ, and we obtain A phase is introduced which does not change the dynamics and could easily be eliminated by a simple rotation. We can now compare the evolution predicted by this model to the numerical integration of H * = H kepl + H res + H syn . The evolution of the action Θ = I 2 is already shown in Figure 6. We plot in Figure 7 the evolution of the angle 2θ (which produces a numerical evolution that is graphically more legible than that of θ). One can see that the angle starts librating at the same value of m pl where the conjugated action Θ = I 2 starts increasing in Figure 6: this shows that there is a passage across the resonance. Note that in such dynamics, the orbit finds itself close to the separatrix after the passage through the resonance, the adiabatic principle is not applicable and the orbit can end up in the inner circulation region (in any case, when the higher order interaction terms between the actions become too strong, a simple description of the dynamics becomes hopeless). In order to get a better sense of the dynamical interaction with this secondary resonance, we can fix different values for δ in (4.48) and look at the corresponding phase diagrams. Notice that changing δ essentially corresponds to changing m pl ; we also checked that for different planetary masses near m pl /M * 1.28 × 10 −3 the coefficients β and c do not change considerably, so we keep them fixed to obtain a qualitatively correct description of the dynamical portraits. Figure 8 shows the level plots of the Hamiltonian (4.48), for different values of δ (i.e. of the frequency of δλ 1,2 /3 + 2φ 2 at Θ = 0), in the variables (X = √ 2Θ cos(2θ), Y = √ 2Θ sin(2θ)); we also overplot the evolution of (X, Y) obtained from the simulation (a combination of Figures 6 and 7), truncated at the value of the planetary mass corresponding to the same δ used to plot the phase diagrams. Initially, there is only one stable centre at the origin (panel (a)) and the orbit circulates anti-clockwise around it with constant amplitude. Then, we see that a resonant island bifurcates from the origin in the bottom-right quadrant of the phase diagram (panel (b)), which is followed by the dynamical evolution. Almost immediately after, a second bifurcation occurs at roughly the same δ, so the inner circulation region starts to grow around the origin and catches up with the orbit (panels (c) and (d)). After crossing the inner separatrix, the dynamical evolution drops off the resonant island, falls inside the inner circulation region and and the angle 2θ starts to circulate in clockwise fashion (panels (e) and (f)). This missed capture into resonance is one of the two probabilistic fates for a second order resonance wheṅ δβ < 0 and when the two bifurcations occur at close values of δ. However, in this specific case we checked that this evolution is actually the result of more complicated interactions among the variables (I l , φ l ) themselves, as well as another secondary resonance involving the variables (I l , φ l ) and (δ∆λ 1,2 , δλ 1,2 ). First, from Figure 6, one can see that after the initial increase of I 2 , I 4 starts increasing also, after which there are wide oscillations of I 2 and I 4 in opposite phase. This is symptomatic of the effect of the term √ I 2 I 4 cos(φ 2 − 2φ 4 ), which is quasi-resonant because ω 2 2ω 4 (Figure 2e). To prove this, we plot in Figure 9 the action I 4 + 2I 2 , which is the constant of motion relative to this harmonic term: we see that the aforementioned coupled oscillations undergone by I 2 and I 4 are completely eliminated. On the other hand, for m pl /M * > 1.297 × 10 −3 we see a much longer period large oscillation, which diverges towards the end of the integration. We interpret this as evidence of a transition of the system into the secondary resonance with argument δλ 1,2 /3 + φ 2 + 2φ 4 (which also has a small frequency, since δλ 1,2 /3 + 2φ 2 and −φ 2 + 2φ 4 are both slow angles). The reason is that (up to a constant) I 4 + 2I 2 can also be seen as a conjugated action of δλ 1,2 /3 + φ 2 + 2φ 4 through the canonical change of variables (I 4 + 2I 2 )/4, δλ 1,2 /3 + φ 2 + 2φ 4 , (I 4 − 2I 2 )/4, − δλ 1,2 /3 − φ 2 + 2φ 4 , (4.49) δ∆λ 1,2 − I 2 /3, δλ 1,2 .
We see that after the angle θ = δλ 1,2 /6 + φ 2 leaves the first resonance (at mass m pl /M * 1.29 × 10 −3 , see Figure 7) the action 2I 2 + I 4 keeps growing, which indicates a transition to this new resonance involving δλ 1,2 /3 + φ 2 + 2φ 4 . This analysis shows that the evolution presented above is very rich, and does not allow any simple description of it. In any case, Figure 8 does not leave any doubt that the initial growth of I 2 is due to the interaction with the secondary resonance associated to the angle δλ 1,2 /3 + 2φ 2 , and that the simple model we have derived yields an effective understanding of the evolution, at least at a qualitative level.
The takeaway is the following. We showed that the numerical integration of the system H * = H kepl + H res + H syn presents an evolution that is similar to that obtained in the full (N + 1)body simulations where the resonant degrees of freedom get excited; we checked that the purely resonant system instead does not undergo the same evolution, and gave an analytical explanation to this fact. We then showed analytically that a set of secondary resonances are present in the H * system, which involve a fraction of the synodic frequency and combination of the resonant frequencies, and which appear at order two in the planetary mass. Then, we found the specific secondary resonance that is encountered in the numerical integration of H * ; we built an integrable model for this resonance valid as long as the actions remain small, and confirmed analytically that there can be capture into this resonance. Finally, we verified that the numerical evolution we obtained in the numerical integration corresponds to a temporary capture into the considered secondary resonance,  Figure 6: the action conjugated to θ is Θ = I 2 (Equation (4.43)), and we see that when θ start librating I 2 increases, indicating that the system has captured into this secondary resonance. As in Figure 6, after the actions get excited the integrable approximation to the dynamics is not valid anymore.
The colour of the dots in this figure only serve as a legend for the value of the planetary mass: we use the same colour-coding in Figure 8, where we take snapshots of the evolution of the pair (Θ, θ) at different values of m pl .  followed by a rich and fascinating series of interactions with additional secondary resonances.

Mass-limit for stability as a function of number of planets and resonance index
We now have all the information needed to derive the general result anticipated at the end of Section 3, namely the dependence of the maximal planetary mass ensuring stability as a function of N and k (i.e. the planet number and resonant index). We sketch below how the results found in the previous section can be generalised to the case of N ≥ 3 equal-mass planets in a given k:k−1 mean motion resonance chain. Following the development presented in Section 4, but for an arbitrary case of N planets, we start introducing the Hamiltonian H * = H kepl + H res + H syn , where H res contains the resonant interactions between all N − 1 pairs of neighbouring planets, and H syn contains terms of type cos(λ i − λ i+1 ), i = 1, . . . , N − 1; for both H res and H syn we will consider interaction terms up to  Fig. 9: Evolution of the action I 4 + 2I 2 , which is a constant of motion relative to the harmonic term φ 2 − 2φ 4 , as well as the resonant action conjugated to the slow angle δλ 1,2 /3 + φ 2 + 2φ 4 , cfr. Eq. (4.49).
order one in the eccentricities, as we did in the previous section. We then make use of a canonical transformation analogous to (4.12, 4.13): we introduce the resonant angles ψ (i) 1 = kλ i+1 − (k − 1)λ i + γ i and the angles δγ i,i+1 = γ i − γ i+1 for each planet pair, i = 1, . . . , N −1 (this gives 2(N −1) resonant degrees of freedom), then we have the synodic angle δλ 1,2 = λ 1 −λ 2 , and finally a γ N = −γ N which will not appear in the Hamiltonian (its conjugated action will again be the total orbital angular momentum).
We can immediately generalise the result of Subsection 4.3 and state that the 2(N − 1) purely resonant degrees of freedom are Lyapunov stable at low amplitude of libration around the resonant equilibrium point for any number of planets N. This is because, when adding an outer resonant pair to the system with the same resonance index k, the Hamiltonian simply repeats itself, since to first order in the eccentricities we are only considering the mutual planetary perturbations due to immediately neighbouring planets and the structure of each term is the same, namely (4.8). So, each planet is either the inner or outer planet, or a middle planet as in the case already considered of a threeplanets system. Therefore, all resonant libration frequencies will always have the same sign, and the reasoning of Subsection 4.3 stands.
As in the case of three resonant planets, we thus conclude that the instability must be due to an interaction between the synodic degree of freedom and the purely resonant degrees of freedom. Then, it is natural to investigate when a regime of secondary resonances analogous to (4.38) and (4.41) can be encountered. To answer this question, we proceeded analytically following the steps of Subsect. 4.4.2. We introduce a generating Hamiltonian χ syn which eliminates the synodic contribution H syn , so χ syn in Delaunay variables will have harmonic terms sin(λ i − λ i+1 ), i = 1, . . . , N − 1. Transforming H * = H kepl + H res + H syn with the Lie series generated by χ syn eliminates H syn to first order in m pl (the planetary mass for all planets), but introduces new terms to order 2 in m pl (among which the most important is {{H kepl , χ syn }, χ syn }, like in the case N = 3); these newly introduced terms will contain a fraction of the synodic angle δλ 1,2 , and we are interested in the smallest fraction of δλ 1,2 that appears. Like in the case of three planets, the term (m 2 pl /2){{H kepl , χ syn }, χ syn } combines together all synodic angles λ i − λ i+1 . Notice that in the new coordinates ψ (i) 1 , δγ i,i+1 and δλ 1,2 , each λ i − λ i+1 can be written as λ i − λ i+1 = k−1 k i−1 (λ 1 − λ 2 ) = k−1 k i−1 δλ 1,2 plus terms including ψ ( j) 1 and δγ j, j+1 . However we do not need to keep track of the ψ ( j) 1 's and δγ j, j+1 's since we are only interested in the way the angle δλ 1,2 appears in the O(m 2 pl ) terms. The smallest fraction of δλ 1,2 will be generated by combining the synodic angles relative to the two outermost pairs λ N−2 − λ N−1 and λ N−1 − λ N , since already they contain the smallest fraction of δλ 1,2 . Multiplying them together (using cos(a) cos(b) = 1 2 (cos(a − b) + cos(a + b))) yields a harmonic term of type where the + . . . terms represents a combination of ψ ( j) 1 's and δγ j, j+1 's, in which again we are not interested. Therefore, the lowest synodic frequency that appears in the O(m 2 ) term is This is the fraction of the synodic frequency which can resonate with the libration frequencies ω l of the resonant degrees of freedom. Since ω l increase with m pl (as ω ∼ m 1/2 pl or m 2/3 pl according to the eccentricities), there will be a critical mass after which a regime of secondary resonances is encountered, which can excite the system and cause its instability. Since the factor 1 k k−1 k N−3 multiplyingδλ 1,2 decreases with increasing N and k, the conclusion is that the regime of secondary resonances between synodic and resonant degrees of freedom is encountered at lower masses for increasing k and/or increasing N, and therefore the critical mass (m pl /M * ) crit allowed for stability decreases with increasing N and k. This gives an analytical explanation to the numerical findings of Matsumoto et al. (2012).

Conclusions
In this paper, we have considered the stability of chains of mean motion resonances, in relation to the observed exoplanet population. Previous works have demonstrated that the paucity of resonances in the exoplanets sample is not in contradiction with the scenario of capture into mean motion resonance during planet migration in the disc phase, if post-disc instability rates are as high as 90% (Izidoro et al. 2017(Izidoro et al. , 2019. This motivates a detailed study on the stability of these chains. Previous numerical investigations pointed out that there is a critical planetary mass above which the instability time of resonant systems is comparable to that of non-resonant ones, and that this limit mass decreases with increasing number of planets and/or increasing index k of the resonance (Matsumoto et al. 2012). The dynamical origin of these instabilities was however not discussed. In this paper we thus investigated analytically and numerically the origin of these instabilities. From the numerical perspective, we used numerical experiments where we first put low-mass planets deep in resonance (at low level of excitation of the resonant modes) and secondly we slowly (and fictitiously) increased the planetary mass to follow the low-amplitude regime until the onset of instability. We confirmed that the instability for three resonant planets occurs at smaller masses than in the two-planet case, and we identified a novel dynamical mechanism which excites the amplitude of libration of the resonant degrees of freedom. The excited systems can then become unstable by suffering close encounters and collisions.
Therefore, we investigated analytically this phenomenon, using a simplified Hamiltonian which reproduced well the observed excitation of the system. Carrying out the calculation explicitly in the case k = 3, we showed that the observed excitation is due to a set of secondary resonances between a combination of the resonant libration frequencies and a fraction of the synodic frequency. We identified the specific secondary resonance that caused the effect observed in the numerical integrations, and built a simple integrable model for this resonance which captures qualitatively the dynamics until the excitation of the system is too severe, showing for example that there can be a capture into this specific resonance. This technique can be generalised to the other secondary resonances.
We therefore proposed that in the numerical simulations the systems become unstable due to a crossing of this type of secondary resonances, which excites the planets' orbits and leads to a phase of close encounters and collisions. This gives a critical mass at which a regime of secondary resonances is encountered, and after which the system can be destabilised. This scheme can then be generalised to an arbitrary number of planets N and/or an arbitrary index of the first-order mean motion resonance k of the chain. One can easily calculate for different N's and k's the lowest fraction of the synodic frequency that can resonate with the resonant frequencies, and see that it decreases with increasing values of N and k (Eq. (5.2)). Consequently, the regime of secondary resonances between synodic and resonant degrees of freedom is encountered when the resonant libration frequencies are slower. Because the libration frequencies grow with the planetary mass, this implies that the instability of the resonant chain occurs at lower masses for increasing k and/or increasing N, and therefore the critical mass allowed for stability decreases with N and with k. This gives an analytical explanation to the numerical findings of Matsumoto et al. (2012).
The takeaway is that we now have a dynamical understanding of the origin of the instabilities observed in the numerical experiments of Matsumoto et al. (2012); Izidoro et al. (2019), which captures the trend in the dependence of the critical mass allowed for stability on the index of the resonance k and the number of planets N. Having understood this mechanism, we will be able to perform a more focused and quantitative analysis on the threshold of stability of resonant chains with different N, k and m pl , and produce an explicit criterion for the stability against secondary resonances of the type described here. This will be the subject of future work. evaluated at the location of the planet. The eccentricity damping is modeled aṡ where τ e is given, in the limit of vanishing eccentricities, by Moreover, we smoothly reverse the sign of the torque at the desired location of the inner edge of the disc in order to stop inward migration, simulating the effects of a cavity.
When two planets are embedded in the disc, the inner planet stops migrating at the edge of the disc and the outer planet continues to migrate inward until it approaches a mean motion resonance with the first, see (Pichierri et al. 2018) for the details. By balancing the eccentricity excitation due to the resonance and the eccentricity damping provided from the disc, one finds that the equilibrium eccentricities at the equilibrium captured state are given by (R 3/2 − 1) τ mig,2 − Re 2 1 τ e,1 − e 2 2 τ e,2 = 0, (A.8) where R = a 2 /a 1 and τ mig,2 is migration rate of the outer planet (cfr. Eq. (A.20) in Pichierri et al. 2018). Since in resonance e 1 ∝ e 2 with a factor that only depends on the masses of the planets, one has that the final equilibrium eccentricity at the capture state is of the order of ∼ H/r, the aspect ratio of the disc. Thus, by changing the aspect ratio of the disc (or equivalently the eccentricity damping timescale) one can reach a resonant configuration with virtually any desired value of the eccentricities. Modifying the disc structure is not an issue here, since the the sole role of the first phase of the numerical experiments described in Section 2 is to obtain a deeply resonant configuration with a desired eccentricity, in order to subsequently study the stability of the obtained configuration as a function of the planetary mass in the second phase.
The case of capture of three (or more) planets in resonance at different desired eccentricities is similar, with only one minor difference. Because the planets capture in resonance in sequence (first planet 1 and 2, then planet 3, etc.) if τ e is large, e 1 and e 2 can grow significantly before planet 3 enters in resonance. This can force large secular eccentricity oscillations of planet 3, which may preclude its resonant capture (see e.g. Batygin (2015) on criteria for resonant capture). We therefore use the following numerical recipe that allows us to capture all N planets at the desired resonance and at any reasonably eccentric configurations. We first capture all planets at small eccentricities, that is with small τ e . Then we slowly increase the value of τ e while the planets are already locked in resonance: since the strength of the resonant interaction stays the same while K = τ a /τ e decreases, Equation (A.8) shows that the planets will adjust to the change in τ e by becoming more eccentric. By doing so adiabatically we obtain, at the fixed initial planetary mass, resonant chains with the same small amplitude of libration around the resonant equilibrium point with different equilibrium eccentricities. This method does not follow the real dynamical evolution of planetary systems, however we reiterate that the role of this first phase is simply to put the planets deeply into a desired resonant chain with a desired eccentricity (i.e. angular momentum) in order to subsequently study the stability of the obtained configuration as a function of the planetary mass.