# The onset of instability in resonant chains

Gabriele Pichierri, Alessandro Morbidelli

► To cite this version:

Gabriele Pichierri, Alessandro Morbidelli. The onset of instability in resonant chains. Monthly Notices of the Royal Astronomical Society, 2020, 494 (4), pp.4950-4968. 10.1093/mnras/sta1102. hal-02996588

HAL Id: hal-02996588

https://hal.science/hal-02996588

Submitted on 4 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The onset of instability in resonant chains

Gabriele Pichierri and Alessandro Morbidelli<fnref n="2" />

<figure><img src="image_3.png" /></figure>

Accepted 2020 April 20. Received 2020 April 20; in original form 2020 March 13

## ABSTRACT

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. Finally, a generalization 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$.

**Key words:** methods: analytical–methods: numerical–celestial mechanics–planets and satellites: dynamical evolution and stability.

## 1 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 none the less 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–20M_⊕$ (Earth's masses) and are found on relatively short orbital periods, of less than about 200 d. 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, Howard & Marcy 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, and, on the other hand (and on longer time-scales), 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 & Papaloizou 2007; Cresswell & Nelson 2008; Morbidelli et al. 2008; Ogihara, Morbidelli & Guillot 2015; Izidoro et al. 2017, 2019; Pichierri, Morbidelli & Crida 2018).

<footnote n="1">¹Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="2">²Laboratoire Lagrange, Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, Boulevard de l'Observatoire, F-06304 Nice Cedex 4, France</footnote>

<footnote>* E-mail: pichierri@mpia.de</footnote>
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, 2017; Luger et al. 2017) and Kepler-223 (Mills et al. 2016). However, an initially puzzling realization 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 characterizing resonant chains (e.g. Winn & 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 evolution following the dissipation of the disc, with instability rates of ~95 per cent. The remaining stable systems naturally represent the observed resonant systems, such as Trappist-1, Kepler-223, etc. In Izidoro et al.'s (2017) original 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 post-disc 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 Myr) and system that suffer dynamical instabilities (collisions or ejections).

The results of Izidoro et al. (2017, 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 the absence of external perturbations.<fnref n="1" /> On this subject, an important numerical study was performed by Matsumoto, Nagasawa & Ida (2012), where the authors studied numerically the stability of resonant multiplanetary 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 Izidoro et al.'s (2017, 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 libration 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 \equiv m_{pl}$, $\forall i = 1, \dots, N$. This is a useful simplification that 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 (Millholland, Wang & Laughlin 2017; Weiss et al. 2018), 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 organized 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) and 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 generalize the analytical scheme to arbitrary $N$ and $k$. Our conclusions are presented in Section 6. Finally, in Appendix A (available as online supplementary material), we summarize the main aspects of the numerical setup, which allows to capture planets into mean motion resonance at different desired eccentricities.

## 2 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.

<footnote n="1">¹External perturbations have also been invoked to increase the fraction of unstable systems, such as the turbulence in the disc (which prevents capture in deep resonance; Batygin & Adams 2017) or the scattering of leftover planetesimals from the planetary region (Chatterjee, Krantzler & Ford 2016).</footnote>
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; also see e.g. Deck & Batygin 2015; Ramos et al. 2017; Xu, Lai & Morbidelli 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 a 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 an 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 Section 2.2.

## 2.1 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 summarize 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 $\omega_{res,i}$, $i=1, 2$. These frequencies dictate the evolution of the system over long time-scales, and are essentially associated with the evolution of the two resonant angles $\psi_i = k\lambda_2 - (k-1)\lambda_1 - \omega_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 time-scales, the evolution is dominated by the non-resonant combination $\delta\lambda_{1,2} = \lambda_1 - \lambda_2$ of the mean longitudes $\lambda_i$, which is a fast-evolving angle; this angle is called the synodic angle, and its frequency is called the synodic frequency $\omega_{syn}$. Since $\lambda_i = n_i$, $n_1/n_2 \simeq k/(k-1)$ by the resonance condition, and $n_i$ is linked to the semi-major axis by Kepler's third law $n_i = \sqrt{GM_*/a^3}$, we have that the synodic frequency $\omega_{syn} = \delta\lambda_{1,2} = n_1/k = \sqrt{GM_*/a^3}/k$. Thus, the synodic frequency is independent of the planetary mass and depends only 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 $\sqrt{m_{pl}}$ (see Section 4.3). Thus, for small enough planetary masses, the synodic frequency is much higher than the resonant frequencies, so that the two contributions happen on totally different time-scales 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, destabilize 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. Marchal & Bozis 1982; Gladman 1993; Petit, Laskar & Boué 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

$$d \ge d_{\text{crit}} = 2\sqrt{3}r_{\text{H},1,2} \simeq 3.46r_{\text{H},1,2}, \quad (1)$$

then the system is Hill-stable (see also Obertas, Van Laerhoven & Tamayo 2017).<fnref n="2" /> Here $r_{\text{H},1,2}$ is the mutual Hill radius of the two planets, defined as

$$r_{\text{H},1,2} = \frac{a_1 + a_2}{2} \left( \frac{m_1 + m_2}{3M_*} \right)^{1/3}. \quad (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

<footnote n="2">²One should note that the resonant condition is a condition on the angles, which prevents the closest approach along the two planets' orbits to happen. This means that resonant systems are expected to be more protected than non-resonant ones with respect to close encounters. However, in Pichierri et al. (2018), we used the actual minimal approach distance $d$ rather than the orbital distance $d_{\text{orb}}$ in (1) to measure the limits of stability against close encounters, and to compare the result with the general case of non-resonant systems.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Numerical investigation of the stability of three planets deep in the 3:2–3:2 mean motion resonance chain, as a function of the planetary mass $m_{pl}$, equal for all planets. In panel (a), four numerical simulations (the coloured markers) are performed starting from low-mass planets ($m_{pl} = 10^{-5}M_*$) and slowly increasing the planetary mass until an instability occurs (a collision in all cases). 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 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_* \simeq 1.28 \times 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 circles. In panel (b), this simulation is also compared with the integration of two simplified models (dark green diamonds and orange triangles), 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 (cf. Section 4.2). The orange triangles represent 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 (cf. Section 4.3). The dark green diamonds represent 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_* \simeq 1.28 \times 10^{-3}$ is well reproduced in this simplified system (cf. Section 4.4).</figcaption></figure>

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 destabilize 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 & 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 \ge 3$ below.

## 2.2 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 time-scale $\tau_e$. By setting a large value for $\tau_e$, large planetary eccentricities could be obtained (cf. equation A8). Here, because the planets capture in resonance in sequence (first planets 1 and 2, then planet 3) if $\tau_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 Fig. 1(a) 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 figs 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 Fig. 1(a) by a dotted line. Comparing the cases $N = 2$ and 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 destabilizing 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 Fig. 1(a), one notices that an interesting phenomenon is evident. For $m_{pl}/M_* < 1.28 \times 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 with the synodic terms (the same effect was present in the case of two planets); instead, the amplitude of libration associated with the purely resonant dynamics is conserved adiabatically. Then, at $m_{pl}/M_* \simeq 1.28 \times 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 realize 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_* \simeq 1.28 \times 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 Section 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 the small libration amplitude for the corresponding planetary mass, which is what was observed by Matsumoto et al. (2012).

### 2.2.2 Numerical and analytical investigation of the phenomenon

In the previous sub-section, 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 that 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 Section 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 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 non-resonant harmonics from the planetary Hamiltonian (cf. Section 4.1) and only keeping resonant harmonics up to some order in $e$. This results in a system governed by a Hamiltonian $\tilde{H} := H_{\text{kepl}} + H_{\text{res}}$; this approach is presented in Section 4.3. By doing so, one realizes 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 Fig. 1(b), where we plot with dark orange triangles the evolution of the system governed by $\tilde{H}$ over one of the full $(N + 1)$-body integrations 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 model $\tilde{H}$), 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 $\tilde{H} = H_{\text{kepl}} + H_{\text{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 $\lambda_1 - \lambda_2$ and $\lambda_2 - \lambda_3$, which we add in an additional interaction Hamiltonian $H_{\text{syn}}$. As we show in Section 4.4, the introduction of these terms is responsible for the same phenomenon observed in Fig. 1(a). This fact is anticipated in Fig. 1(b), where we plot with darker-coloured diamonds the evolution of the system governed by $H^* := H_{\text{kepl}} + H_{\text{res}} + H_{\text{syn}}$ over one of the full $(N + 1)$-body integrations 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 (cf. Section 4.4.2). In the light of this, we will use the evolution yielded by the simplified model $H^* = H_{\text{kepl}} + H_{\text{res}} + H_{\text{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 Section 4.4.1) and to understand what is the dynamical 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 Section 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_{pl}^2)$, 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 Section 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 generalized to the other secondary resonances that can, in principle be encountered. Finally, we proceed to generalize this result to more populated and/or more compact resonant chains in Section 5.

## 3 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 Fig. 1. In particular, we aim at gaining a deep understanding of the process, which causes the sudden excitation in the systems shown in Fig. 1. As anticipated at the end of the last subsection, this process involves secondary resonances between some fraction of the synodic frequency $\omega_{syn} = \frac{d}{dt}(\lambda_1 - \lambda_2)$ and the resonant frequencies $\omega_{res,i}$ 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 $\omega_{syn}$) and the purely resonant evolution (with characteristic frequency $\omega_{res,i} \ll \omega_{syn}$) happen on such different time-scales that there can be no interaction between them (as we already recalled in Section 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 $\omega_{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 (cf. Section 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 \ge 3$, there is a fraction of the synodic frequency that appears in the Hamiltonian (in terms at second order in the planetary masses). In the case of three planets, this fraction is $\omega_{syn}/k$ where $k$ as usual is the index of the resonance.<fnref n="3" /> Thus, in the case of three planets, in order to reach a secondary resonance

involving synodic and resonant degrees of freedom, the resonant frequencies do not have to be as large, that is, the planetary masses do not have to be as large as in the two-planet case. This is why for $N \ge 3$, these secondary resonances can be relevant while they were not in the case $N = 2$.

To extend this principle to the general case $N \ge 3$, one can easily calculate that the smallest fraction of the synodic frequency that appears in the case of $N$ planets in a $k:k-1$ resonant chain is $\frac{1}{k}(\frac{k-1}{k})^{N-3}\omega_{syn}$ (cf. equation 52). Again, this frequency can resonate with the resonant frequencies $\omega_{res}$, and, just as before, $\frac{1}{k}(\frac{k-1}{k})^{N-3}\omega_{syn} \approx \frac{1}{k^2}(\frac{k-1}{k})^{N-3}n_1$ is independent of the planetary masses, and, for fixed orbital separation (fixed $n_1$), decreases with $k$ and $N$. Finally, the resonant frequencies $\omega_{res}$ still increase with $m_{pl}$ (and with $k$), more or less independently of the number of planets. Thus, there will be a critical mass after which a regime of secondary resonances is encountered, which can excite the system and cause its subsequent instability by close encounters. Since the factor $\frac{1}{k}(\frac{k-1}{k})^{N-3}$ multiplying $\omega_{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 $k$. This mechanism gives a dynamical explanation to the numerical findings of Matsumoto et al. (2012) and Izidoro et al. (2019). In the rest of this paper, we give a detailed analytical description of the dynamical emergence of these secondary resonances.

## 4 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 Section 4.1, and we then consider the relevant harmonic terms in the Hamiltonian that interest us in Section 4.2. Then, in Sections 4.3 and 4.4, respectively, we consider the averaged model $\mathcal{H}$ and the model $\mathcal{H}^*$ that includes synodic terms. There, we give an analytical descriptions of the main dynamical features of the simulations shown in Fig. 1.

### 4.1 Planetary Hamiltonian

We start with the Hamiltonian $\mathcal{H}$ of $N$ planets of masses $m_i$, $i = 1, \dots, N$ orbiting a star of mass $M_*$. We let $\mathbf{u}_i$ be the inertial barycentric Cartesian coordinate of each planet, and $\tilde{\mathbf{u}}_i = m_i\tilde{\mathbf{u}}_i$ the conjugated momentum. We write $\mathcal{H}$ in canonical heliocentric variables $(\mathbf{p}_i, \mathbf{r}_i)$, $i = 1, \dots, N$, defined from the inertial barycentric canonical variables $(\mathbf{u}, \tilde{\mathbf{u}}_i)$ as

$$ \begin{aligned} \mathbf{p}_0 &= \sum_{i=0}^{N} \tilde{\mathbf{u}}_i, & \mathbf{r}_0 &= \mathbf{u}_0, \\ \mathbf{p}_i &= \tilde{\mathbf{u}}_i, & \mathbf{r}_i &= \mathbf{u}_i - \mathbf{u}_0, & i &= 1, \dots, N \end{aligned} \quad (3) $$

(e.g. Poincaré 1892; Laskar 1990). Doing so, the Hamiltonian can be split as

$$ \begin{aligned} \mathcal{H}(\mathbf{p}, \mathbf{r}) &= \mathcal{H}_{\text{kepl}} + \mathcal{H}_{\text{pert}}, \\ \mathcal{H}_{\text{kepl}} &= \sum_{i=1}^{N} \left( \frac{\|\mathbf{p}_i\|^2}{2\mu_i} - \frac{G(M_* + m_i)\mu_i}{\|\mathbf{r}_i\|} \right) = \sum_{i=1}^{N} \mathcal{H}_{\text{kepl},i}, \\ \mathcal{H}_{\text{pert}} &= \sum_{1 \le i < j \le N} \left( \frac{\mathbf{p}_i \mathbf{p}_j}{M_*} - \frac{G m_i m_j}{\|\mathbf{r}_i - \mathbf{r}_j\|} \right). \end{aligned} \quad (4) $$

In other words, the Hamiltonian appears as a sum of two terms. One term is the sum of the Keplerian unperturbed

<footnote n="3">³A simple explanation for why this fraction of the synodic frequency naturally pops up in the equations of motion (that is in the Hamiltonian) at second order in $m_{pl}$ is the following. The Hamiltonian of three planets contains both the $\delta\lambda_{1,2} = \lambda_1 - \lambda_2$ and $\delta\lambda_{2,3} = \lambda_2 - \lambda_3$ harmonics. If both planet pairs are in the $k:k-1$ mean motion resonance, one can write $\delta\lambda_{2,3}$ as $(k-1)\delta\lambda_{1,2}/k$ plus some correction harmonic terms that depend only on the resonant angles (cf. 16) with constant index $k$ along the chain); this can be easily understood by noting that $\delta\lambda_{2,3}$ should be comparable to $(k-1)\delta\lambda_{1,2}/k$ in a $k:k-1$ chain. Then, the two angles $\delta\lambda_{1,2}$ and $(k-1)\delta\lambda_{1,2}/k$ get combined at second order in $m_{pl}$, which yields a harmonic containing $\delta\lambda_{1,2}/k$ plus purely resonant harmonics (cf. 39).</footnote>
Hamiltonians for each planet $\mathcal{H}_{\text{kepl},i}$, describing the planet-star interactions. The other is the perturbing Hamiltonian $\mathcal{H}_{\text{pert}}$ that describes all planet-planet interactions; $\mathcal{H}_{\text{pert}}$ itself is split into direct terms, $-\sum_{1 \le i < j \le N} G m_i m_j / \|\mathbf{r}_i - \mathbf{r}_j\|$, and indirect terms, $\sum_{1 \le i < j \le N} \mathbf{p}_i \mathbf{p}_j / M_*$, which come from having considered canonical heliocentric rather than barycentric variables. $\mathcal{H}_{\text{kepl}} = \sum_{i=1}^N \mathcal{H}_{\text{kepl},i}$ is integrable, while $\mathcal{H}_{\text{pert}}$ is of order $m_{\text{pl}}/M_*$ with respect to $\mathcal{H}_{\text{kepl}}$ (where $m_{\text{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 $\mathcal{H}_{\text{kepl},i}$. We will consider only coplanar motion for the planets, so we only have two pairs of action-angle variables ($\Lambda_i, \lambda_i$) and ($\Gamma_i, \gamma_i$). Their definition in terms of the orbital elements is (e.g. Morbidelli 2002)

$$ \Lambda_i = \mu_i \sqrt{G(M_* + m_i)a_i}, \quad \lambda_i = \ell_i + \bar{\omega}_i, \\ \Gamma_i = \Lambda_i (1 - \sqrt{1 - e_i^2}) \sim \Lambda_i e_i^2 / 2, \quad \gamma_i = -\bar{\omega}_i. \qquad (5) $$

As usual, for each planet, $a_i$ is the semi-major axis, $e_i$ is the eccentricity, $\lambda_i$ is the mean longitude, $\ell_i$ is the mean anomaly $\bar{\omega}_i$ is the longitude of the pericentre, and $\mu_i = m_i M_*/(M_* + m_i) \simeq m_i$ is the reduced mass; the index $i=1, \dots, 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 (3) (they are the so-called formal osculating elements; Morbidelli 2002).

In the modified Delaunay variables (5), the Keplerian part can be rewritten as

$$ \mathcal{H}_{\text{kepl}} = -G^2 \sum_{i=1}^{N} \frac{\mu_i^3 (M_* + m_i)^2}{2\Lambda_i^2}, \qquad (6) $$

while no simple expression exists for $\mathcal{H}_{\text{pert}}$, which is usually expanded in Fourier series of the angles. In this expansion, there are only combinations of $\lambda_i$ and $\gamma_i$ that satisfy the d'Alembert characteristics, and only harmonic terms combining angles from two planets. We will not go into the details of how this expansion is performed in general, which can be found in many works (e.g. Laskar & Robutel 1995; Murray & Dermott 1999), and we will only concentrate on the specific terms that interest us below.

## 4.2 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.

First, since the instabilities for $N \ge 3$ planets occur at much lower values of $m_{\text{pl}}/M_*$ than for two planets, we approximate the reduced mass $\mu = m_{\text{pl}} M_*/(M_* + m_{\text{pl}}) \sim m_{\text{pl}}$ and $M_* + m_{\text{pl}} \sim M_*$. Then, we recall that all the planets have the same mass $m_{\text{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_{\text{pl}}$, plus a small perturbation proportional to $m_{\text{pl}}$. The natural choice is to rescale all the actions ($\Lambda, \Gamma$) of the modified Delaunay variables by the planetary mass $m_{\text{pl}}$, which yields

$$ \Lambda = \sqrt{GM_*a}, \quad \Gamma = \Lambda(1 - \sqrt{1 - e^2}), \qquad (7) $$

where, for simplicity, we have maintained the same notation as for the non-rescaled variables. In order to maintain the canonicity of the Hamiltonian, $\mathcal{H}$ itself must be rescaled by $m_{\text{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)

$$ \mathcal{H} = \mathcal{H}_{\text{kepl}} + \mathcal{H}_{\text{pert}}, \quad \mathcal{H}_{\text{kepl}} = -\sum_{i=1}^{N} \frac{G^2 M_*^2}{2\Lambda_i^2}, \qquad (8) $$

where $\mathcal{H}_{\text{kepl}}$ is independent of $m_{\text{pl}}$, and the (rescaled) perturbation is of order $O(m_{\text{pl}})$:

$$ \mathcal{H}_{\text{pert}} = m_{\text{pl}} \mathcal{H}_{\text{pert}'} \qquad (9) $$

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

$$ \mathcal{H}_{\text{res}}^{(i)} = m_{\text{pl}} \left[ \alpha_1^{(i)} e_i \cos (k^{(i)} \lambda_{i+1} - (k^{(i)} - i)\lambda_i + \gamma_i) + \alpha_2^{(i)} e_{i+1} \cos (k^{(i)} \lambda_{i+1} - (k^{(i)} - i)\lambda_i + \gamma_{i+1}) \right], \qquad (10) $$

where the (rescaled) coefficients $\alpha$ are

$$ \alpha_j^{(i)} = - \frac{G^2 M_*}{\Lambda_{i+1}^2} f_{\text{res}}^{(j,i)}(\alpha_{\text{res}}^{(i)}), \qquad (11) $$

where $f_{\text{res}}^{(j,i)}(\alpha_{\text{res}}^{(i)})$ are functions of the Laplace coefficients $b_s^{(j)}$ (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, $\alpha_{\text{res}}^{(i)} = \bar{a}_i / \bar{a}_{i+1} = ((k^{(i)} - 1)/k^{(i)})^{2/3}$ is the nominal semi-major axial 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 $\Lambda_{i+1}^2$ at the denominator at its nominal Keplerian value $\bar{\Lambda}_{i+1}$ [because the terms in (11) are already of order $O(e)$, e.g. Batygin & Morbidelli 2013]. Doing so, the coefficients $\alpha_j^{(i)}$ 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 $\mathcal{H}_{\text{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 $\lambda_i - \lambda_{i+1}$, they take the form

$$ \mathcal{H}_{\text{syn}}^{(i)} = c_i \cos(\lambda_i - \lambda_{i+1}) = m_{\text{pl}} C_i \cos(\lambda_i - \lambda_{i+1}), \qquad (12) $$

where the coefficients $C_i$ for the rescaled Hamiltonian are

$$ C_i = - \frac{G^2 M_*}{\bar{\Lambda}_{i+1}^2} \times \left[ \frac{1}{2} b_{1/2}^{(i)} (\alpha_{\text{res}}^{(i)}) - (\alpha_{\text{res}}^{(i)})^{-1/2} \right], \qquad (13) $$

and have the same scaling in $\bar{\Lambda}_{i+1}$ as the coefficients in (11) but a different dependence on the Laplace coefficients $b_s^{(j)}$ (e.g. Murray & Dermott 1999; the term $-(\alpha_{\text{res}}^{(i)})^{-1/2}$ comes from the indirect term of the perturbing function). Notice that (12) is of order 0 in eccentricity. The term $O(e)$ cannot exist because it would not satisfy the d'Alembert rules. So, (12) is all we have for the terms dependent on the difference of the mean longitudes of neighbouring planets $\lambda_i - \lambda_{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^{(i)}\lambda_{i+1} - (k^{(i)} - i)\lambda_i + \gamma_i) + j(\lambda_i - \lambda_{i+1})$, for an arbitrary integer $j$). Because, in what follows, they would behave like those in $\mathcal{H}_{\text{syn}}$
in (12) but are O(e) smaller, we neglect them for simplicity. Notice also that in (12) we can limit ourselves to the lowest multiples of $\lambda_i - \lambda_{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 non-resonant synodic angles $\lambda_i - \lambda_{i+1}$, so it is preferable that one of them, say $\lambda_1 - \lambda_2$, be also one of the canonical variables. The natural choice is to use as canonical positions the resonant angles $\psi_1^{(i)} = \theta^{(i)} + \gamma_i$ (where $\theta^{(i)} = k^{(i)}\lambda_{i+1} - (k^{(i)} - 1)\lambda_i$ is the longitude of conjunction for the ith pair) and the apsidal differences $\delta\gamma_{i,i+1} = \gamma_i - \gamma_{i+1}$ for $i = 1, \dots, N-1$, then define $\delta\lambda_{1,2} = \lambda_1 - \lambda_2$ and finally keep an angle that will not appear explicitly in the Hamiltonian, such as $\gamma_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

$$ 
\begin{aligned}
\psi_1^{(1)} &= k^{(1)}\lambda_2 - (k^{(1)} - 1)\lambda_1 + \gamma_1, \\
\psi_1^{(2)} &= k^{(2)}\lambda_3 - (k^{(2)} - 1)\lambda_2 + \gamma_2, \\
\delta\gamma_{1,2} &= \gamma_1 - \gamma_2, \\
\delta\gamma_{2,3} &= \gamma_2 - \gamma_3, \\
\delta\lambda_{1,2} &= \lambda_1 - \lambda_2, \\
\gamma_3' &= -\gamma_3,
\end{aligned} 
\quad (14) 
$$

while the new conjugated actions are<fnref n="4" />

$$ 
\begin{aligned}
\psi_1^{(1)} &= \Lambda_1 + \Lambda_2 + \frac{k^{(2)} - 1}{k^{(2)}} \Lambda_3, \\
\psi_1^{(2)} &= \frac{1}{k^{(2)}} \Lambda_3, \\
\Delta\gamma_{1,2} &= -\left(\Lambda_1 + \Lambda_2 + \frac{k^{(2)} - 1}{k^{(2)}} \Lambda_3\right) + \Gamma_1, \\
\Delta\gamma_{2,3} &= -(\Lambda_1 + \Lambda_2 + \Lambda_3) + \Gamma_1 + \Gamma_2, \\
\Delta\lambda_{1,2} &= k^{(1)}\Lambda_1 + (k^{(1)} - 1)\Lambda_2 + \frac{(k^{(1)} - 1)(k^{(2)} - 1)}{k^{(2)}} \Lambda_3 = k^{(1)}\mathcal{K}, \\
\mathcal{L} &= (\Lambda_1 + \Lambda_2 + \Lambda_3) - (\Gamma_1 + \Gamma_2 + \Gamma_3).
\end{aligned} 
\quad (15) 
$$

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 generalizable 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 $\gamma'_N = -\gamma_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 (8) with increasing $m_{pl}$, L will always remain constant. Moreover, the action $\Delta\lambda_{1,2}$ conjugated to the angle $\delta\lambda_{1,2}$ is simply a factor away from the action K used in Pichierri et al. (2018; see also e.g. Batygin & Morbidelli 2013a); this action has been called the 'spacing parameter' (Michtchenko, Beaugé & Ferraz-Mello 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 $\bar{\Lambda}$ of the resonance at hand.

We note that in these variables in the case of three resonant planets in a $k^{(1)}$: ($k^{(1)} - 1$) - $k^{(2)}$: ($k^{(2)} - 1$) chain, the synodic harmonics for the two pairs of planets can be written as

$$ 
\begin{aligned}
\lambda_1 - \lambda_2 &= \delta\lambda_{1,2}, \\
\lambda_2 - \lambda_3 &= \frac{1}{k^{(2)}}((k^{(1)} - 1)\delta\lambda_{1,2} + \psi_1^{(1)} - \psi_1^{(2)} - \delta\gamma_{1,2}).
\end{aligned} 
\quad (16) 
$$

Then, in the new variables, the Keplerian Hamiltonian (8), the resonant contribution (10), and the synodic contribution (12) for pairs of neighbouring planets write

$$ 
\begin{aligned}
\mathcal{H}_{\text{kepl}} = & -\frac{G^2 M_*^2}{2(-\Delta\lambda_{1,2} + k^{(1)}\psi_1^{(1)} - k^{(2)}\psi_1^{(2)} + \psi_1^{(2)})^2} \\
& -\frac{G^2 M_*^2}{2(\Delta\lambda_{1,2} - k^{(1)}\psi_1^{(1)} + \psi_1^{(1)})^2} - \frac{G^2 M_*^2}{2(k^{(2)}\psi_1^{(2)})^2},
\end{aligned} 
\quad (17a) \\
\mathcal{H}_{\text{res}} = m_{\text{pl}}\mathcal{H}_{\text{res}'} \quad (17b) 
$$

$$ 
\begin{aligned}
\mathcal{H}_{\text{syn}} = m_{\text{pl}} & \left[ C_1 \cos(\delta\lambda_{1,2}) \right. \\
& \left. + C_2 \cos \left( \frac{1}{k^{(2)}} \left( (k^{(1)} - 1)\delta\lambda_{1,2} + \psi_1^{(1)} - \psi_1^{(2)} - \delta\gamma_{1,2} \right) \right) \right] \\
& = m_{\text{pl}} \mathcal{H}_{\text{syn}'}.
\end{aligned} 
\quad (17c) 
$$

$\mathcal{H}_{\text{kepl}}$ is independent of $m_{\text{pl}}$ and depends on the variables $\psi_1^{(1)}$, $\psi_1^{(2)}$, and $\Delta\lambda_{1,2}$ only; one can introduce the frequencies (analogous to the mean motions $n$)

$$ \eta_{\psi_1^{(1)}} := \frac{\partial \mathcal{H}}{\partial \psi_1^{(1)}}, \quad \eta_{\psi_1^{(2)}} := \frac{\partial \mathcal{H}}{\partial \psi_1^{(2)}}, \quad \eta_{\Delta\lambda_{1,2}} := \frac{\partial \mathcal{H}}{\Delta\lambda_{1,2}}. \quad (18) $$

$\mathcal{H}_{\text{res}}$ depends only on the angles through the harmonic terms $\cos \psi_1^{(1)}$, $\cos(\psi_1^{(1)} - \delta\gamma_{1,2})$, $\cos \psi_1^{(2)}$, and $\cos(\psi_1^{(2)} - \delta\gamma_{2,3})$, and each term has a coefficient depending on the actions (15) and the coefficients $\alpha_j^{(i)}$; the exact expression can easily be obtained by direct substitution. In (17b) and (17c), we use a prime (') to indicate that the Hamiltonian term has been rescaled by $m_{\text{pl}}$ itself, so it is O(0) in $m_{\text{pl}}$, and the dependence on $m_{\text{pl}}$ has been clearly expressed with a coefficient. In the following, we will also use the notation

$$ 
\begin{aligned}
\mathbf{x} = & (\psi_1^{(1)}, \psi_1^{(2)}, \Delta\gamma_{1,2}, \Delta\gamma_{2,3}, \Delta\lambda_{1,2}, \psi_1^{(1)}, \psi_1^{(2)}, \delta\gamma_{1,2}, \delta\gamma_{2,3}, \\
& \times \delta\lambda_{1,2})
\end{aligned} 
\quad (19) 
$$

for the canonical variables that enter in $\mathcal{H}$ (except the pair $(\mathcal{L}, \gamma'_3)$, since $\gamma'_3$ does not appear in $\mathcal{H}$ and $\mathcal{L}$ is a constant of motion); we write for the actions $\mathbf{p} = (\psi_1^{(1)}, \psi_1^{(2)}, \Delta\gamma_{1,2}, \Delta\gamma_{2,3}, \Delta\lambda_{1,2})$ and for the angles $\mathbf{q} = (\psi_1^{(1)}, \psi_1^{(2)}, \delta\gamma_{1,2}, \delta\gamma_{2,3}, \delta\lambda_{1,2})$.

<footnote n="4">⁴A note on notation can help clarify the meaning of the names of these variables. Reading the definitions for the angles, the upper indices (i) refer to which pair of planets is considered, so that $k^{(1)}$ is the index of the first-order mean motion resonance for the inner pair and $\psi^{(1)}$ refers to a resonant angle for that pair. Similarly, $\theta^{(i)} = k^{(i)}\lambda_{i+1} - (k^{(i)} - 1)\lambda_i$ is the longitude of conjunction of the ith pair. The subscript 1 in $\psi_1^{(i)}$ signifies the fact that for each pair, we choose to use the resonant angle that depends on the longitude of pericentre $\gamma$ of the innermost planet of the pair, so $\psi_1^{(i)} = \theta^{(i)} + \gamma_i$, while the other resonant angle would then be $\psi_2^{(i)} = \theta^{(i)} + \gamma_{i+1}$ (and it does not appear since we also use $\gamma_i - \gamma_{i+1}$ as canonical angles). The conventions for the other angles are evident. For the actions, we simply use an uppercase first letter to indicate to which angle each action is conjugated, except for the last action since it is just the orbital angular momentum, which we always indicate with $\mathcal{L}$, and it is always a constant of motion.</footnote>
## 4.3 Purely resonant dynamics

The purely resonant dynamics is the one governed by the Hamiltonian averaged over the fast angle $\delta\lambda_{1,2}$, i.e. $\tilde{H} = H_{\text{kepl}} + H_{\text{res}}$. $\tilde{H}$ is now rewritten in terms of the new canonical variables (14, 15) (cf. equation 17), and since the synodic terms have been removed by averaging, not only $\mathcal{L}$ but also $\Delta\lambda_{1,2}$ is a constant of motion, so that only the ‘barred’ variables $\bar{x} = (\bar{p}, \bar{q}) = (\psi_1^{(1)}, \psi_1^{(2)}, \Delta\gamma_{1,2}, \Delta\gamma_{2,3}, \psi_1^{(1)}, \psi_1^{(2)}, \delta\gamma_{1,2}, \delta\gamma_{2,3})$ evolve.<fnref n="5" /> These barred variables $\bar{x} = (\bar{p}, \bar{q})$ are simply a subset of the variables $\mathbf{x} = (\mathbf{p}, \mathbf{q})$ introduced above, and represent the purely resonant degrees of freedom (this notation is only introduced to separate these variables from the synodic canonical pair ($\Delta\lambda_{1,2}, \delta\lambda_{1,2}$); this will be a useful distinction later on).

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_{\text{pl}}/M_* = 10^{-5}$ and we choose as initial condition that of Fig. 1(b). The resulting evolution of the canonical actions $\bar{p}$ is shown in dark green triangles in Fig. 2, panels (a)–(d) (the evolution of the eccentricity has been already presented in Fig. 1b). We observe that the four resonant degrees of freedom are never unstable even up to masses significantly higher than the critical mass ($m_{\text{pl}}/M_*$)crit $\simeq 1.28 \times 10^{-3}$, which is found in the numerical ($N+1$)-body simulations with the same initial conditions (Fig. 1b, light green circles in Fig. 2).

We can present an analytical explanation for this. As in Pichierri et al. (2018), we find the stable resonant equilibrium points for $\tilde{H}(\bar{x}; \mathcal{L}, \Delta\lambda_{1,2}, m_{\text{pl}})$ in the variables $\bar{x}$, while keeping $\mathcal{L}$ and $\Delta\lambda_{1,2}$ constants and for different values of $m_{\text{pl}}$, yielding $\bar{x}_{\text{eq}}(m_{\text{pl}}) = \bar{x}_{\text{eq}}(m_{\text{pl}}; \mathcal{L}, \Delta\lambda_{1,2})$. Notice that at these low eccentricities, we are interested in symmetric linearly stable equilibria only, so the equilibrium values $\bar{q}_{\text{eq}}$ of the angles are simply

$$ \psi_{1,\text{eq}}^{(1)} = 0, \quad \psi_{1,\text{eq}}^{(2)} = 0, \quad \delta\gamma_{1,2,\text{eq}} = \pi, \quad \delta\gamma_{2,3,\text{eq}} = \pi, \qquad (20) $$

and we only need to solve for the equilibrium actions $\bar{p}_{\text{eq}} = (\psi_{1,\text{eq}}^{(1)}, \psi_{2,\text{eq}}^{(1)}, \Delta\gamma_{1,2,\text{eq}}, \Delta\gamma_{2,3,\text{eq}})$. In Fig. 2, panels (a)–(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 diagonalize the system around the equilibrium point $\bar{x}_{\text{eq}}$; since it is a stable equilibrium point, all eigenvalues are purely imaginary and the diagonalization procedure yields a Hamiltonian of the form

$$ \tilde{H}(\xi, \eta) = \sum_{l=1}^{4} \frac{\omega_l}{2} (\xi_l^2 + \eta_l^2) + \mathcal{O}(\|(\xi, \eta)^3\|) \qquad (21) $$

in Cartesian coordinates $\bar{x} = T(\xi, \eta)$, with $T$ being a transformation matrix. Using canonical polar coordinates $(I_l, \phi_l)_{l=1,...,4}$ with $(\xi_l = \sqrt{2}I_l \cos \phi_l, \eta_l = \sqrt{2}I_l \sin \phi_l)$, we get

$$ \tilde{H} = \sum_{l=1}^{4} \omega_l I_l + \mathcal{O}(\|\mathbf{I}^{3/2}\|), \qquad (22) $$

which appears as the sum of four decoupled harmonic oscillators plus higher order terms.

The resulting four frequencies $\omega_l$, $l = 1, ..., 4$ are shown in Fig. 2(c) as a function of the planetary mass $m_{\text{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.

Since it will be useful later on, we also consider here how the libration frequencies grow with $m_{\text{pl}}$. This is shown in Fig. 2, panels (e) and (f). We find numerically that $\omega_{1,2} \propto m_{\text{pl}}^{2/3}$ at low eccentricities ($e \simeq 0.01$, panel e) while $\omega_{1,2} \propto m_{\text{pl}}^{1/2}$ at higher eccentricities ($e \simeq 0.1$, panel f). Notice that for a pendulum-type Hamiltonian like

$$ \mathcal{H}_{\text{pend}}(\Sigma, \sigma) = a\Sigma^2 - m_{\text{pl}}b \cos \sigma, \qquad (23) $$

the libration frequency would be $\propto m_{\text{pl}}^{1/2}$, 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_{\text{pl}}$, and the real scaling would therefore be $\sqrt{abm_{\text{pl}}}$. The way the equilibrium points adjust to changes in $m_{\text{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 (23), as $m_{\text{pl}}^{1/3}$. 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.

## 4.4 The synodic contribution

In the previous sub-section, 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 the lowest order in $e$, we introduce the two synodic terms (12) for the inner and outer pairs that had been averaged out before, resulting in

$$ \mathcal{H}_{\text{syn}} = m_{\text{pl}} [C_1 \cos(\lambda_1 - \lambda_2) + C_2 \cos(\lambda_2 - \lambda_3)] = m_{\text{pl}} \mathcal{H}_{\text{syn}'} \qquad (24) $$

with coefficients given by (13). The full rescaled Hamiltonian written in the new variables (14, 15) is now

$$ \mathcal{H}^*(\mathbf{x}; \mathcal{L}, m_{\text{pl}}) = \mathcal{H}_{\text{kepl}} + \mathcal{H}_{\text{res}} + \mathcal{H}_{\text{syn}}. \qquad (25) $$

We have stressed that it depends parametrically on the constant of motion $\mathcal{L}$ and on the mass $m_{\text{pl}}$ through $\mathcal{H}_{\text{res}} = m_{\text{pl}}\mathcal{H}_{\text{res}'}$ and $\mathcal{H}_{\text{syn}} = m_{\text{pl}}\mathcal{H}_{\text{syn}'}$.

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 triangles in Fig. 3, which is matched against the ($N+1$)-body integration with the same initial datum (lighter green circles) and the locations of the equilibria for $\tilde{H}$ calculated in the previous section for different planetary

<footnote n="5">⁵From a technical point of view, these variables are only an approximation (up to order 0 in $m_{\text{pl}}$) to the actual canonical variables that eliminate the non-resonant contributions. These would be the primed variables introduced later on in (36).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The purely resonant evolution governed by $\tilde{\mathcal{H}}$ in the case of three planets in a 3:2–3:2 mean motion resonance chain, with the same initial conditions as in Fig. 1(b). We show in dark green triangles in panels (a)–(d) the evolution of the actions $\bar{p} = (\psi_1^{(1)}, \psi_2^{(1)}, \Delta\gamma_{1,2}, \Delta\gamma_{2,3})$ as the planetary mass $m_{\text{pl}}$ is slowly increasing, and we match it to the calculated equilibria $\bar{p}_{\text{eq}} = (\psi_{1,\text{eq}}^{(1)}, \psi_{2,\text{eq}}^{(1)}, \Delta\gamma_{1,2,\text{eq}}, \Delta\gamma_{2,3,\text{eq}})$ (purple dashed line); we also add the corresponding $(N+1)$-body integration with the same initial condition (light green circles). A legend for panels (a)–(d) is given in panel (a). We see that the system remains stable well after the value of $m_{\text{pl}}/M_* \simeq 1.28 \times 10^{-3}$ corresponding to the onset of excitation in Fig. 1(b). Panels (e) and (f) contain the analytical explanation of the observed stability: We plot with coloured continuous 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)–(d), $e \simeq 0.01$; in panel (f), we used higher initial eccentricities, $e \simeq 0.1$. We note that the scaling law for $\omega_l(m_{\text{pl}})$ changes depending on the eccentricity (see black dotted and dashed lines).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Panels (a)–(e) show in dark green triangles the evolution of the actions $\boldsymbol{p} = (\psi_1^{(1)}, \psi_1^{(2)}, \Delta\gamma_{1,2}, \Delta\gamma_{2,3}, \Delta\lambda_{1,2})$ for the system $\mathcal{H}^* = \mathcal{H}_{\text{kepl}} + \mathcal{H}_{\text{res}} + \mathcal{H}_{\text{syn}}$ with the same initial condition as the $(N+1)$-body integration of Fig. 1(b) [the evolution of these variables in the $(N+1)$-body integration is also shown here in light green circles for reference]. The system follows on average the purple dashed lines, which correspond to the equilibria $\boldsymbol{p}_{\text{eq}}$ for the system $\tilde{\mathcal{H}} = \mathcal{H}_{\text{kepl}} + \mathcal{H}_{\text{res}}$. With an orange continuous line, we show the evolution of the averaged variables $\boldsymbol{p}'$ calculated through analytical averaging of the fast synodic frequencies (equation 37). Note that, for $m_{\text{pl}}/M_* < 1.28 \times 10^{-3}$, $\boldsymbol{p}'$ has very little oscillation around the $\boldsymbol{p}_{\text{eq}}$ curve, compared to the $\boldsymbol{p}$ evolution. Instead, for $m_{\text{pl}}/M_* > 1.28 \times 10^{-3}$, the amplitude of oscillations of $\boldsymbol{p}'$ and $\boldsymbol{p}$ around $\boldsymbol{p}_{\text{eq}}$ are almost the same. This reveals that, while the initial oscillation of $\boldsymbol{p}$ is entirely due to the synodic terms and is effectively removed by passing to the $\boldsymbol{p}'$ variables, it is then dominated by an increased amplitude of libration in the resonance. The evolution of the angular momentum $\mathcal{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.</figcaption></figure>
masses (purple dashed lines). The comparison for the eccentricity evolutions, instead of the canonical variables, has been already presented in Fig. 1(b). We notice two important aspects of these plots. The first is that, initially, for all variables the evolution described by $\mathcal{H}^*$ follows on average that described by $\tilde{\mathcal{H}}$ (compare Fig. 2 with Fig. 3, and the dashed purple lines). This can be easily understood realizing that $\mathcal{H}^*$ contains fast, non-resonant angles, which, up to first order in the small parameter $m_{pl}$, have simply been averaged out in $\tilde{\mathcal{H}}$; therefore, as long as the $\mathcal{O}(m_{pl}^2)$ contributions are unimportant, the only difference between the two evolutions are the short-periodic, $\mathcal{O}(m_{pl})$ oscillations due to the $\delta\lambda_{1,2}$ synodic angle. We will actually study this effect analytically below. However, as soon as the $\mathcal{O}(m_{pl}^2)$ remainder introduces important contributions to the dynamics, as in the case of the emergence of a secondary resonance, the dynamics described by the averaged $\tilde{\mathcal{H}}$ approximation is not valid anymore. This is indeed what we see in Fig. 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_* \simeq 1.28 \times 10^{-3}$, which was not found in $\tilde{\mathcal{H}}$. Notice that such a secondary resonance cannot be caused by an interaction of the resonant degrees of freedom $\tilde{x}$ 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 ($\Delta\lambda_{1,2}, \delta\lambda_{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.

### 4.4.1 Eliminating the $\mathcal{O}(m_{pl})$ synodic term

In the previous section, we dropped $\mathcal{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}\mathcal{H}_{syn}'$ from $\mathcal{H}^*$. This will introduce $\mathcal{O}(m_{pl}^2)$ terms that we want to calculate explicitly, since they contain harmonics mixing $\tilde{q}$ and $\delta\lambda_{1,2}$, potentially associated with secondary resonances. In order to eliminate $m_{pl}\mathcal{H}_{syn}'$ at $\mathcal{O}(m_{pl})$, we need to find a generating function $\chi_{syn}$ that solves the homological equation

$$ \{\chi_{syn}, \mathcal{H}_{kepl}\} + \mathcal{H}_{syn} = 0. \qquad (26) $$

In the above equation, we used the Poisson bracket $\{\bullet, \bullet\}$ that operates on two dynamical variables $f_1$ and $f_2$, yielding a new dynamical variable $\{f_1, f_2\}$ defined by

$$ \{f_1, f_2\} := (\nabla f_1)^T \mathbf{J} (\nabla f_2), \qquad (27) $$

where $\nabla$ is the gradient with respect to the canonical variables and $\mathbf{J}$ is the standard symplectic matrix. Clearly, $\chi_{syn}$ will be of order $m_{pl}$ so we can write $\chi_{syn} = m_{pl}\chi_{syn}'$. From the expression for $\mathcal{H}_{syn}$ (equation 17c), we see that $\chi_{syn}$ will have the form

$$ \chi_{syn} = m_{pl} \left[ \frac{C_1}{\eta_1} \sin(\delta\lambda_{12}) + \frac{C_2}{\eta_2} \sin \left( \frac{1}{k^{(2)}} \left( (k^{(1)} - 1)\delta\lambda_{12} + \psi_1^{(1)} - \psi_1^{(2)} - \delta\gamma_{1,2} \right) \right) \right], \qquad (28) $$

where the divisors $\eta_1$ and $\eta_2$ are immediately found in terms of the frequencies (18) of the unperturbed Keplerian Hamiltonian and the combination of angles appearing in the harmonics in $\mathcal{H}_{syn}$, yielding

$$ \eta_1 = \eta_{\Delta\lambda_{1,2}}, \quad \eta_2 = \frac{1}{k^{(2)}} \left( (k^{(1)} - 1)\eta_{\Delta\lambda_{1,2}} + \eta_{\psi_1^{(1)}} - \eta_{\psi_1^{(2)}} \right). \qquad (29) $$

These divisors are not vanishing nor small, since clearly $\eta_1 = n_1 - n_2$, $\eta_2 = n_2 - n_3$ (remember that $n_i$ is the mean motion frequency of planet $i$ and that the harmonics in $\mathcal{H}_{syn}$ in the modified Delaunay variables were simply $\lambda_1 - \lambda_2$ and $\lambda_2 - \lambda_3$) and the planets are evidently far from the 1:1 resonance. Therefore, equation (26) can indeed be solved.

Having calculated $\chi_{syn}$, we can then write out how the Hamiltonian $\mathcal{H}^*$ transforms under the Lie series transformation $\exp(L_{\chi_{syn}})$ generated by $\chi_{syn}$. Here, $\exp(L_{\chi_{syn}})$ is given by

$$ \exp(L_{\chi_{syn}}) f = f + m_{pl} \{f, \chi_{syn}'\} + \frac{m_{pl}^2}{2} \{\{f, \chi_{syn}\}, \chi_{syn}'\} + \dots \qquad (30) $$

The new Hamiltonian $\mathcal{H}'$ is given by $\exp(L_{\chi_{syn}})\mathcal{H}^*$, and reads up to $\mathcal{O}(m_{pl}^2)$:

$$ \mathcal{H}' = [\mathcal{H}_{kepl}] + m_{pl}\{\mathcal{H}_{kepl}, \chi_{syn}'\} + \frac{m_{pl}^2}{2}\{\{\mathcal{H}_{kepl}, \chi_{syn}'\}, \chi_{syn}'\} + \dots \qquad (31a) $$

$$ [m_{pl}\mathcal{H}_{res}'] + m_{pl}^2\{\mathcal{H}_{res'}, \chi_{syn}'\} + \dots \qquad (31b) $$

$$ + m_{pl}\mathcal{H}_{syn}' + m_{pl}^2\{\mathcal{H}_{syn'}, \chi_{syn}'\} + \dots; \qquad (31c) $$

as it is typical in perturbation theory via Lie transform, this transformed Hamiltonian is written in terms of the new variables $x'$ by direct substitution of the old variables $x$ to the new, $x \to x'$. The change of variable is given by $x = \exp(L_{\chi_{syn}})x'$ (see equations 36 and 37). We note that the boxed terms are simply

$$ [\mathcal{H}_0 + m_{pl}\mathcal{H}_{res}']|_{x'} =: \tilde{\mathcal{H}}|_{x'}, \qquad (32) $$

that is, $\tilde{\mathcal{H}}$ written in the new variables $x'$ via direct substitution. Recall that $\tilde{\mathcal{H}}$ does not depend on $\delta\lambda_{1,2}$ and so $\Delta\lambda_{1,2}$ was a first integral; hence, only the 'averaged variables' $\tilde{x}' = (\psi_1^{(1)}, \psi_1^{(2)}, \Delta\gamma_{1,2}, \Delta\gamma_{2,3}, \psi_1^{(1)}, \psi_1^{(2)}, \delta\gamma_{1,2}, \delta\gamma_{2,3})$ remain as evolving variables (as in Section 4.3, we use a barred notation $\tilde{x}' = (\tilde{p}', \tilde{q}')$ for the purely resonant variables, the subset of $\tilde{x}'$ not including $(\Delta\lambda_{1,2}, \delta\lambda_{1,2})$). Concerning the remaining two $\mathcal{O}(m_{pl})$ terms in (31), these actually cancel out by construction, since $\chi_{syn}$ was chosen to satisfy (26). We can, therefore, write the transformed Hamiltonian as

$$ \mathcal{H}' = [\mathcal{H}_{kepl}|_{x'} + m_{pl}\mathcal{H}_{res'}|_{x'}] + \mathcal{O}(m_{pl}^2). \qquad (33) $$

This equation shows that $\tilde{\mathcal{H}}$, which we have studied in the previous section, approximates $\mathcal{H}'$ only to first order in $m_{pl}$. As long as the new term $\mathcal{O}(m_{pl}^2)$ does not contain resonant terms, the variables $\tilde{p}' = (\psi_1^{(1)}, \psi_1^{(2)}, \Delta\gamma_{1,2}, \Delta\gamma_{2,3})$ closely follow the equilibrium points $\tilde{p}_{eq}$ calculated from $\tilde{\mathcal{H}}$ in Section 4.3 (i.e. they have oscillations around $\tilde{p}_{eq}$ of order $\mathcal{O}(m_{pl}^2)$, while the oscillations of $\tilde{p}$ are $\mathcal{O}(m_{pl})$), while $\Delta\lambda_{1,2}$ undergoes oscillations of $\mathcal{O}(m_{pl}^2)$ around the initial value $\Delta\lambda_{1,2}$ (again a conserved quantity in the purely averaged model $\tilde{\mathcal{H}}$). This is what we observe in the numerical simulations;
see Fig. 3. We can therefore simplify the calculation by writing

$$
\chi_{\text{syn}'} = \underbrace{\chi_{\text{syn}}(\bar{p}' = \bar{p}_{\text{eq}}; \mathcal{L}, \Delta\lambda'_{1,2} = \overline{\Delta\lambda_{1,2}})}_{\bar{\chi}_{\text{syn}}} \\
+ \mathcal{O}(|(\bar{p}' - \bar{p}_{\text{eq}}, \Delta\lambda'_{1,2} - \overline{\Delta\lambda_{1,2}})|), \qquad (34)
$$

(where we called $\bar{\chi}'_{\text{syn}}$ the first term of the last equation), and drop-
ping the higher order terms, which correspond to small deviations
from $\bar{p}' = \bar{p}_{\text{eq}}$ and from the initial value $\Delta\bar{\lambda}_{1,2}$ of $\Delta\bar{\lambda}'_{1,2}$. With this
approximation we can eliminate the term $m_{\text{pl}}^2\{\mathcal{H}_{\text{syn}'}$, $\bar{\chi}'_{\text{syn}}\}$ in (31c)
because now $\frac{\partial \mathcal{H}_{\text{syn}'}}{\partial p} = \frac{\partial \bar{\chi}'_{\text{syn}}}{\partial p} = 0$, so $\{\mathcal{H}_{\text{syn}'}$, $\bar{\chi}'_{\text{syn}}\} = \frac{\partial \mathcal{H}_{\text{syn}'} \frac{\partial \bar{\chi}'_{\text{syn}}}{\partial q}}{\partial q} -$
$\frac{\partial \mathcal{H}_{\text{syn}'} \frac{\partial \bar{\chi}'_{\text{syn}}}{\partial q}}{\partial p} = 0$ (of course, by the same reasoning, also the higher
order terms of the Lie series for $\exp(L_{\bar{\chi}_{\text{syn}}})\mathcal{H}_{\text{syn}}$ cancel out). The
resulting Hamiltonian becomes

$$
\mathcal{H}' = \left[ \mathcal{H}_{\text{kepl}}|_{x'} + m_{\text{pl}} \mathcal{H}_{\text{res}}|_{x'} \right] \\
+ \left[ \frac{m_{\text{pl}}^2}{2} \left\{ \mathcal{H}_{\text{kepl}}, \bar{\chi}'_{\text{syn}} \right\}_{x'} + m_{\text{pl}}^2 \left\{ \mathcal{H}_{\text{res}}, \bar{\chi}'_{\text{syn}} \right\}_{x'} \right] + \dots \tag{35}
$$

We now explicit the transformation that to $\mathcal{O}(m_{\text{pl}})$ eliminates the
fast synodic evolution in the numerical integrations. This is given
by

$$
\begin{align}
\boldsymbol{p} &= \exp(L_{\bar{\chi}_{\text{syn}}}) \boldsymbol{p}' = \boldsymbol{p}' + m_{\text{pl}} \{\boldsymbol{p}', \bar{\chi}'_{\text{syn}}\} + \mathcal{O}(m_{\text{pl}}^2) \nonumber \\
\boldsymbol{q} &= \exp(L_{\bar{\chi}_{\text{syn}}}) \boldsymbol{q}' = \boldsymbol{q}'. \tag{36}
\end{align}
$$

Notice that the angles remain unchanged since $\bar{\chi}_{\text{syn}}$ is independent
of the actions, so $\frac{\partial \bar{\chi}_{\text{syn}}}{\partial p'} = 0$. The transformation for the actions
reads, to first order in $m_{\text{pl}}$:

$$
\begin{align}
\psi_1^{(1)} &= \psi_1^{(1)'} - m_{\text{pl}} \frac{1}{k^{(2)}} \frac{C_2}{\bar{\eta}_2} \cos(\delta\lambda_{2,3}), \\
\psi_1^{(2)} &= \psi_1^{(2)'} + m_{\text{pl}} \frac{1}{k^{(2)}} \frac{C_2}{\bar{\eta}_2} \cos(\delta\lambda_{2,3}), \\
\Delta\gamma_{1,2} &= \Delta\gamma_{1,2}^{\prime} + m_{\text{pl}} \frac{1}{k^{(2)}} \frac{C_2}{\bar{\eta}_2} \cos(\delta\lambda_{2,3}), \\
\Delta\gamma_{2,3} &= \Delta\gamma_{2,3}^{\prime}, \\
\Delta\lambda_{1,2} &= \Delta\lambda_{1,2}^{\prime} - m_{\text{pl}} \left[ \frac{C_1}{\bar{\eta}_1} \cos(\delta\lambda_{1,2}) + \frac{k^{(1)} - 1}{k^{(2)}} \frac{C_2}{\bar{\eta}_2} \cos(\delta\lambda_{2,3}) \right], \tag{37}
\end{align}
$$

where one has to replace $\delta\lambda_{2,3}$ with its expression in terms of the
variables (14), $\delta\lambda_{2,3} = \frac{1}{k^{(2)}}(k^{(1)} - 1)\delta\lambda_{1,2} + \psi_1^{(1)} - \psi_1^{(2)} - \delta\gamma_{1,2}$;
moreover, $\bar{\eta}_1$ and $\bar{\eta}_2$ are the frequencies (29) evaluated at the
reference values for the actions at each $m_{\text{pl}}$. We can invert these
expressions to obtain $p'$ from $(p, q)$, and the evolution of $p'$
represents that of $p$ where, to first order in $m_{\text{pl}}$, the short periodic
have been averaged out. The evolution of $p'$ is shown as an orange
continuous line in Fig. 3 in our reference $N=3$, $k=3$ example,
where we see that, initially, the averaged evolution follows closely
the analytical calculation of the equilibrium points of $\bar{\mathcal{H}}$ for different
planetary masses.

This is however only valid until a point in which the $\mathcal{O}(m_{\text{pl}}^2)$ con-
tribution, which is still present in (35), has resonant effects (which
happens at $m_{\text{pl}}/M_* \simeq 1.28 \times 10^{-3}$ in Fig. 3). Indeed, as it is typical
in perturbation theory, these terms are expected to contain higher
order harmonics that were not present in the original Hamiltonian
$\mathcal{H}^* = \mathcal{H}_{\text{kepl}} + \mathcal{H}_{\text{res}} + \mathcal{H}_{\text{syn}}$: then, if these newly introduced $\mathcal{O}(m_{\text{pl}}^2)$

Hamiltonian terms contains angles that, for certain values of $m_{\text{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.

4.4.2 The $\mathcal{O}(m_{\text{pl}}^2)$ contribution

In this sub-section, we look closely at the $\mathcal{O}(m_{\text{pl}}^2)$ terms in
(35). We are specifically interested in the harmonics that they
contain, to find explicitly which combinations of angles $\bar{q}' =$
$(\psi_1^{(1)'}, \psi_1^{(2)'}, \delta\gamma_{1,2}', \delta\gamma_{2,3}')$ and $\delta\lambda'_{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 $\delta\lambda'_{1,2}$ is much higher
than the libration frequencies characteristic of the angles $\bar{q}'$, the
most interesting harmonics are the ones where the lowest fraction
of $\delta\lambda'_{1,2}$ appears next to a combinations of $\bar{q}'$. 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 Fig. 2 panels (e) and (f), at the lowest planetary mass. The
following calculation is clearly general, but to simplify matters, we
will quickly specialize 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 Fig. 1
were performed.

We start with the main term $\{\mathcal{H}_{\text{kepl}}, \bar{\chi}'_{\text{syn}}\}$, $\bar{\chi}'_{\text{syn}}$ of order $m_{\text{pl}}^2$
in (35). Since $\mathcal{H}_{\text{kepl}}$ does not contain any angles, all secondary
resonance contributions must come from combinations of the
harmonics contained in $\bar{\chi}'_{\text{syn}}$. Recall that we defined $\bar{\chi}'_{\text{syn}}$ containing
both synodic terms with harmonics $\lambda_1 - \lambda_2$ and $\lambda_2 - \lambda_3$, which we
wrote in equation (16) in terms of the new variables $q$. Therefore, the
harmonics that are included in $\{\mathcal{H}_{\text{kepl}}, \bar{\chi}'_{\text{syn}}\}$, $\bar{\chi}'_{\text{syn}}$ are combinations
of these synodic harmonics; more specifically, they come from
the products of their cosines.<fnref n="6" /> Using the standard trigonometric
identity $\cos(a)\cos(b) = \frac{1}{2}(\cos(a-b) + \cos(a+b))$, the resulting
harmonics are

$$
\begin{align}
& 2\delta\lambda'_{1,2}, \nonumber \\
& ((k^{(2)} + k^{(1)} - 1)\delta\lambda'_{1,2} + \psi_1^{(1)'} - \psi_1^{(2)'} - \delta\gamma'_{1,2})/k^{(2)}, \tag{38} \\
& ((k^{(2)} - k^{(1)} + 1)\delta\lambda'_{1,2} - \psi_1^{(1)'} + \psi_1^{(2)'} + \delta\gamma'_{1,2})/k^{(2)}, \nonumber \\
& 2((k^{(1)} - 1)\delta\lambda'_{1,2} + \psi_1^{(1)'} - \psi_1^{(2)'} - \delta\gamma'_{1,2})/k^{(2)}, \nonumber
\end{align}
$$

so the harmonic with the lowest fraction of $\delta\lambda'_{1,2}$ is $((k^{(2)} - k^{(1)} + 1)\delta\lambda'_{1,2} - \psi_1^{(1)'} + \psi_1^{(2)'} + \delta\gamma'_{1,2})/k^{(2)}$. Specializing now to the case of a chain with the same resonance index $k^{(1)} = k^{(2)} = k$, this simply gives

$$
\frac{1}{k} (\delta\lambda'_{1,2} - \psi_1^{(1)'} + \psi_1^{(2)'} + \delta\gamma'_{1,2}). \qquad (39)
$$

With the aid of an algebraic manipulator, one can compute the full
expression of $\{\mathcal{H}_{\text{kepl}}, \bar{\chi}'_{\text{syn}}\}$, $\bar{\chi}'_{\text{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

<footnote n="6">⁶This can be easily understood noting that if $\chi = \sin(q_1) + \sin(q_2)$ and $f$
depends only on the actions $p_i$, then</footnote>

<footnote>$$
\{\{f, \chi\}, \chi\} = \frac{\partial^2 f}{\partial p_1^2} \cos^2 q_1 + 2 \frac{\partial^2 f}{\partial p_1 \partial p_2} \cos q_1 \cos q_2 + \frac{\partial^2 f}{\partial p_2^2} \cos^2 q_2.
$$</footnote>
solely from the term $\propto \{(1/\Lambda_2^2, \tilde{\chi}'_{\text{syn}}, \tilde{\chi}'_{\text{syn}})\}$. We avoid writing here the full expression, which is rather cumbersome; moreover, as in (35), 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

$$ \mathcal{H}_{\text{scnd.res.kepl}} = \text{const} \times m_{\text{pl}}^2 \cos \left( \frac{(\delta\lambda'_{1,2} - \psi_1^{(1)'} + \psi_1^{(2)'} + \delta\gamma'_{1,2})}{k} \right). \quad (40) $$

Since we want to compare the frequency of $\delta\lambda'_{1,2}/k$ to that of $(-\psi_1^{(1)'} + \psi_1^{(2)'} + \delta\gamma'_{1,2})/k$, we need to consider the resonant Hamiltonian $\tilde{H}$ in the $x'$ variables, expand the ‘barred’ variables $\tilde{x}'$ around the equilibrium point characterized by the equilibrium actions $\tilde{p}_{\text{eq}}$ and the equilibrium angles $\tilde{q}_{\text{eq}}$ (equation 20) as in Section 4.3, and then introduce the transformation $\tilde{x}' \to (I_l, \phi_l)_{l=1,...,4}$ to the action-angle variables $(I, \phi)$, which transforms $\tilde{H}$ into the sum of decoupled harmonic oscillators plus higher order terms (equation 22). It is also useful to translate the value of $\Delta\lambda'_{1,2}$ around its initial reference value $\overline{\Delta\lambda'_{1,2}}$ introducing $\Delta\lambda'_{1,2} = \overline{\Delta\lambda'_{1,2}} + \delta\Delta\lambda'_{1,2}$. Therefore, we write $\mathcal{H}_{\text{scnd.res.kepl}}$ in terms of the variables $(I, \delta\Delta\lambda'_{1,2}, \phi, \delta\lambda'_{1,2})$. The Hamiltonian $\mathcal{H}_{\text{scnd.res.kepl}}$ will now contain harmonic terms of type

$$ \left\{ \frac{\sin}{\cos} \right\} \left( \frac{\delta\lambda'_{1,2}}{k} + \mathbf{h}\boldsymbol{\phi} \right), \quad (41) $$

where $\mathbf{h}\boldsymbol{\phi}$ is an integer combination with coefficients $h_1, \dots, h_4 \in \mathbb{Z}$ of the angles $\phi_1, \dots, \phi_4$, which we can calculate explicitly. Therefore, whenever $\frac{d}{dt}(\delta\lambda'_{1,2}/k) = -\frac{d}{dt}(\mathbf{h}\boldsymbol{\phi})$, a secondary resonance is crossed. We can rewrite this expression as $\delta\lambda'_{1,2}/k + \mathbf{h}\boldsymbol{\omega} = 0$. Since the Hamiltonian has d'Alembert characteristics in each pair $(I_l, \phi_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 $\mathbf{h}\boldsymbol{\phi}$, that is, where most $h_l$ are zero. We also note that since $\delta\lambda'_{1,2} > 0$ and the frequencies $\omega_l$ are all negative, a secondary resonance term can only appear when $\mathbf{h}\boldsymbol{\omega} < 0$, which, together with the requirement that $|\mathbf{h}|$ be small, is tantamount to requiring that all non-zero integers $h_l$ are positive. Since we calculated $\omega(m_{\text{pl}})$ in Section 4.3, we can calculate for each $\mathbf{h}$ the relative frequency $(\delta\lambda'_{1,2}/k + \mathbf{h}\boldsymbol{\omega})(m_{\text{pl}})$ as a function of $m_{\text{pl}}$, and check if any of these vanish for some value of $m_{\text{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 \approx 0.1$, which corresponds to the evolution shown in Fig. 1(b) (and also Figs 2 and 3). We found that $\mathcal{H}_{\text{scnd.res.kepl}}$ contains, among many others, the following terms:

$$ 1.24 \times m_{\text{pl}}^2 \sqrt{2I_1} \left\{ \frac{\sin}{\cos} \right\} \left( \frac{\delta\lambda'_{1,2}}{3} + \phi_1 + \text{phase} \right), \quad (42a) $$

$$ 0.27 \times m_{\text{pl}}^2(2I_2) \left\{ \frac{\sin}{\cos} \right\} \left( \frac{\delta\lambda'_{1,2}}{3} + 2\phi_2 + \text{phase} \right), \quad (42b) $$

$$ 2.39 \times 10^{-3} \times m_{\text{pl}}^2 \sqrt{2I_1} \sqrt{2I_3} \left\{ \frac{\sin}{\cos} \right\} \left( \frac{\delta\lambda'_{1,2}}{3} + \phi_1 + \phi_3 + \text{phase} \right), \quad (42c) $$

$$ 1.6 \times m_{\text{pl}}^2 \sqrt{2I_1} \sqrt{2I_4} \left\{ \frac{\sin}{\cos} \right\} \left( \frac{\delta\lambda'_{1,2}}{3} + \phi_1 + \phi_4 + \text{phase} \right). \quad (42d) $$

The nature of these harmonics is clearly general, while the numerical

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Frequencies of the angles $\delta\lambda'_{1,2}/3 + \mathbf{h}\boldsymbol{\phi}$ as a function of the planetary mass in the case of the 3:2–3:2 mean motion resonance chain with $a_1 \approx 0.1$ (the situation depicted in Fig. 1b). Notice that the synodic frequency $\delta\lambda'_{1,2}$ varies only slightly due to the change in the equilibrium point $\bar{x}_{\text{eq}}$ for the averaged Hamiltonian $\tilde{H}$, which is followed by the full system $\tilde{H}'$ until the second-order effects become significant (cf. equation 35). The main change comes from the resonant frequencies $\omega$, whose dependence on the planetary mass is depicted in Fig. 2(c). The result is that the frequencies $(\delta\lambda'_{1,2}/k + \mathbf{h}\boldsymbol{\omega})(m_{\text{pl}})$ vanish within a small range of values of the planetary mass $m_{\text{pl}}$, meaning that a capture into a secondary resonance becomes possible. By comparing with Fig. 1(b), we see that $\delta\lambda'_{1,2}/3 + 2\phi_2$ has vanishing frequency at the same value of $m_{\text{pl}}/M_* \approx 1.28 \times 10^{-3}$ at which the excitation of the system occurs.</figcaption></figure>

coefficients are specific to the reference case $k=3$ and $a_1 \approx 0.1$ mentioned above. We then calculated for each of the harmonics in (42) their frequency $(\delta\lambda'_{1,2}/k + \mathbf{h}\boldsymbol{\omega})(m_{\text{pl}})$ as a function of the mass. The results are presented in Fig. 4.

We immediately remark that in the case of the harmonic $\delta\lambda'_{1,2}/3 + 2\phi_2$, the crossing of the secondary resonance happens precisely at the value of planetary mass $m_{\text{pl}}/M_* \approx 1.28 \times 10^{-3}$ where the numerical integrations showed the increase in amplitude of libration (see Figs 1–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.

First, if we had used in $\mathcal{H}_{\text{syn}}$ only one synodic term, not all of the harmonics in (38) would appear.⁷ In particular, the harmonic $\delta\lambda'_{1,2}/3 + 2\phi_2$ would not appear, so that the observed dynamical effects linked to the crossing of secondary resonances at $m_{\text{pl}}/M_* \approx 1.28 \times 10^{-3}$ are not expected. Indeed, we performed similar numerical integrations with only one of the synodic terms, $\lambda_1 - \lambda_2$, and separately $\lambda_2 - \lambda_3$, which are shown in Fig. 5, and there is no effect at the right value of $m_{\text{pl}}$. Secondary resonances do occur, but at larger values of $m_{\text{pl}}$, given that the generated harmonics have a larger coefficient for $\delta\lambda'_{1,2}$.

Secondly, so far, we have not considered the $\mathcal{O}(m_{\text{pl}}^2)$ term $\{\mathcal{H}_{\text{res.'}, \tilde{\chi}'_{\text{syn}}}\}$, which is also present in (35). However, with the same technique as above one can see that this term only yields harmonics

⁷To see this, as in footnote 6, we calculate for $\chi = \sin(q_1)$ and $f$, which depends on the actions only:

$$ \{f, \chi\}, \chi\} = \frac{\partial^2 f}{\partial p_1^2} \cos^2 q_1 = \frac{1}{2} \frac{\partial^2 f}{\partial p_1^2} (1 + \cos(2q_1)). $$

Clearly, we do not obtain the needed $(k^{(2)} - k^{(1)} + 1)\delta\lambda'_{1,2} - \psi_1^{(1)'} + \psi_1^{(2)'} + \delta\gamma'_{1,2})/k^{(2)}$ in (38) neither when $q_1 = \lambda_1 - \lambda_2 = \delta\lambda_{1,2}$ nor when $q_1 = \lambda_2 - \lambda_3 = \frac{1}{k^{(2)}}(k^{(1)} - 1)\delta\lambda_{1,2} + \psi_1^{(1)} - \psi_1^{(2)} - \delta\gamma_{1,2})$.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison between the full (N + 1)-body simulation from Fig. 1(b) (lightest green circles) and two numerical simulations with the same initial conditions where only one of the two synodic terms $\lambda_1 - \lambda_2$ and $\lambda_3$ appears (darker green triangles and diamonds). These two semi-synodic simulations initially appear identical, which is easily understood from the fact that they follow in average the evolution of $\tilde{H} = H_{\text{kepl}} + H_{\text{res}}$, which is the same for the two (cf. equation 33). The important point is that in both cases, when only one synodic angle is considered, the system is not excited at value of $m_{\text{pl}}/M_* \simeq 1.28 \times 10^{-3}$, where it is excited in the (N + 1)-body simulation as well as in the numerical simulation that includes both synodic terms; see Fig. 1(b). This shows that both synodic terms must be included in order to have a good quantitative agreement with the (N + 1)-body simulations.</figcaption></figure>

of type

$$ 
\begin{aligned}
\delta\lambda'_{1,2} \pm \psi_1^{(1)'} & \\
\delta\lambda'_{1,2} \pm \psi_1^{(2)'} & \\
\pm\delta\lambda'_{1,2} - \psi_1^{(1)'} + \delta\gamma'_{1,2} & \\
((k^{(1)}-1)\delta\lambda'_{1,2} - (\pm k^{(2)} - 1)\psi_1^{(1)'} - \psi_1^{(2)'} + (\pm k^{(2)}-1)\delta\gamma'_{1,2})/k^{(2)} & \\
((k^{(1)} - 1)\delta\lambda'_{1,2} + \psi_1^{(1)'} - (\pm k^{(2)} + 1)\psi_1^{(2)'} - \delta\gamma'_{1,2})/k^{(2)} & \\
((k^{(1)}-1)\delta\lambda'_{1,2} + \psi_1^{(1)'} - (\pm k^{(2)} + 1)\psi_1^{(2)'} - \delta\gamma'_{1,2} \pm k^{(2)}\delta\gamma'_{2,3})/k^{(2)} &
\end{aligned} 
\quad (43) $$

Whenever $k^{(1)} \ge 3$, as in our reference case $k^{(1)} = k^{(2)} = k = 3$, this does not contribute the needed harmonic (39) with $\delta\lambda_{1,2}$ appearing as a single $\delta\lambda_{1,2}/k$; it will only include multiples of $\delta\lambda_{1,2}/k$ and therefore to lowest order does not contribute to the secondary resonance harmonics in (42).

Finally, in (35), we used the simplification $\bar{p}' = \bar{p}_{\text{eq}}$, $\Delta\lambda_{1,2}' - \bar{\Delta}\lambda_{1,2}$ to define $\bar{\chi}_{\text{syn}}$ (cf. equation 34). However, the remaining terms of $\mathcal{O}(|\bar{p}' - \bar{p}_{\text{eq}}|, \Delta\lambda_{1,2}' - \bar{\Delta}\lambda_{1,2}|)$ do not contribute to the dynamics to lowest order. Indeed, concerning $m_{\text{pl}}\{\mathcal{H}_0, \chi\}|_{x'}$, this term only contains the two separate synodic harmonics already contained in $\chi$ and therefore does not yield terms linked to secondary resonances. Finally, the remaining terms in $\frac{m_{\text{pl}}^2}{2}[(\mathcal{H}_{\text{kepl}} \cdot \chi_{\text{syn}}'), \chi_{\text{syn}}]|_{x'}$ will only yield higher order terms in the actions $I$, so we can neglect them (recall that initially the values of the actions are small since we are close to the equilibrium point).

With these clarifications, we can proceed with the model of the secondary resonance linked to the angle $\delta\lambda'_{1,2}/3 + 2\phi_2$, which, as we discussed above, has vanishing frequency exactly at the value of $m_{\text{pl}}$ when the increase in the amplitude of libration is observed in Fig. 1(b). This realization is further supported by Fig. 6. There, we plot the evolution of the actions $I$, $l = 1, \dots, 4$ along the simulation, with the planetary mass $m_{\text{pl}}$ on the horizontal axis. We see that initially only one action is excited, namely $I_2$, and after that the non-linearities 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 $\delta\lambda'_{1,2}/3 + 2\phi_2$

In the following we detail how we can construct a model for the resonance associated with the angle $\delta\lambda'_{1,2}/3 + 2\phi_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 (42).

We start by performing a canonical transformation which selects $\delta\lambda'_{1,2}/3 + 2\phi_2$ as an angle. Notice that, since $\phi_2$ appears with a coefficient 2 and so $\sqrt{2}I_2$ appears as a power 2 in (42b), we have a secondary resonance of order 2; hence, it is useful to define the resonant angle $\theta$ as $2\theta = \delta\lambda'_{1,2}/3 + 2\phi_2$ in order to maintain the d'Alembert characteristics so that the Hamiltonian will not be singular at the origin. The resulting transformation is

$$ 
\begin{aligned}
\Theta &= I_2, & \theta &= \delta\lambda'_{1,2}/6 + \phi_2, \\
I_r^* &= I_r, & \phi_r^* &= \phi_r, & r &= 1, 3, 4, \\
\delta\Delta\lambda_{1,2}^* &= \delta\Delta\lambda'_{1,2} - I_2/6 & \delta\lambda_{1,2}^* &= \delta\lambda'_{1,2}, & (44)
\end{aligned} $$

whose canonicity follows immediately from the preservation of the Poisson brackets. We can already notice that $\Theta = I_2$ appears as the conjugated action to the angle $\theta$ associated with the secondary resonance: this explains why in Fig. 6 it is $I_2$, which is initially excited. The other variables do not feel the resonance, except $\delta\Delta\lambda'_{1,2}$, which must change according to the change in $I_2$ in order to maintain $\delta\Delta\lambda'_{1,2}$ constant; however, since $I_2$ gets divided by 6, this change is minute, but nevertheless clearly visible in Fig. 3e. The pair $(\Theta, \theta)$ 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-dof system that we write $\langle \mathcal{H}' \rangle_{(\phi^*, \delta\lambda^*_{1,2})}(\Theta, \theta; I^*, \delta\Delta\lambda^*_{1,2})$. The notation $\langle \cdot \rangle_{(\phi^*, \delta\lambda^*_{1,2})}$ means that we eliminated perturbatively to lowest order the non-secondary-resonant contributions from the angles $(\phi^*, \delta\lambda^*_{1,2})$, and we stressed that the variables $(I^*, \delta\Delta\lambda^*_{1,2})$ will only play the role of parameters for $\langle \mathcal{H}' \rangle_{(\phi^*, \delta\lambda^*_{1,2})}(\Theta, \theta)$ will be that of a Andoyer Hamiltonian, that is

$$ \langle \mathcal{H}' \rangle_{(\phi^*, \delta\lambda^*_{1,2})}(\Theta, \theta) = \delta\Theta + \frac{\beta}{2}\Theta^2 + \mathcal{O}(\Theta^3) + c(\sqrt{2\Theta})^2 \cos(2\Theta); \quad (45) $$

the coefficient $c$ will be of order $m_{\text{pl}}^2$, while $\delta$ and $\beta$ will be of order $m_{\text{pl}}$. Since the system is initially close to the resonant equilibrium point, $\Theta$ is small and we can drop the $\mathcal{O}(\Theta^3)$ terms. However, as we will see below, the parameter $\beta$ (the second derivative at $\Theta = 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 $\mathcal{O}(\Theta^2)$ terms of the $\theta$-independent part, that is, the first two terms in (45). The main contribution to the $\theta$-independent part comes from the $\tilde{H}$ term [the $\mathcal{O}(m_{\text{pl}})$ term in (35)], while $c\sqrt{2\Theta^2}\cos(2\Theta)$ comes from (42b) and is $\mathcal{O}(m_{\text{pl}}^2)$. Concerning the first part deriving from $\tilde{H}$, we should stress that even if $\delta\Delta\lambda'_{1,2}$ appeared as a constant of motion when this Hamiltonian was treated alone, when the $\mathcal{O}(m_{\text{pl}}^2)$ is taken into account, the transformation (44) transforms $\delta\Delta\lambda'_{1,2}$
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Evolution around $m_{pl}/M_* \simeq 1.28 \times 10^{-3}$ of the actions $I_l, l = 1, \dots, 4$ along the reference numerical simulation shown in Fig. 1(b). We see that the actions are initially relatively constant, and the system is well approximated by a Hamiltonian of the form $\sum_{l=1}^4 \omega_l I_l$ (cf. equation 22). Then, $I_2$ increases steadily, which is a symptom of an interaction with a secondary resonance that involves $\Theta = I_2$ as a resonant action; this is confirmed by the canonical change of coordinates (44). Soon after $I_2$ is large enough, the degrees of freedom start interacting and exchanging energy, due to the non-linear effects.</figcaption></figure>

into $\delta\Delta\lambda_{1,2}^* + \Theta/6$, where $\delta\Delta\lambda_{1,2}^*$ is the new constant of motion. Therefore, we must keep $\delta\Delta\lambda_{1,2}'$ as a variable in $\tilde{H}$ and apply (44) 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 $\delta\beta < 0$. Intuitively, this is because near the origin one has $\dot{\theta} \simeq \delta + \beta\Theta$ and the resonance condition imposes that (on average) this quantity vanishes; therefore, at the centre of the resonance $\Theta = -\delta/\beta$, which only makes sense when $\beta$ and $\delta$ have opposite signs. Thus, since $\beta$ remains relatively constant (see below), only when $\delta\beta$ is negative does the resonance centre appear from the origin and move at higher values of $\Theta$, while if $\delta\beta$ 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 Fig. 4 that, as the planetary mass increases, $\dot{\theta}$ goes from positive values to negative values, that is, $\dot{\delta} < 0$: this means a capture into this secondary resonance is possible only if $\beta > 0$.

To obtain the sign of $\beta$ in (45), 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_* \simeq 1.28 \times 10^{-3}$, we calculate the equilibrium point $\tilde{x}_{eq} = \tilde{x}_{eq}(m_{pl})$ and we apply the canonical diagonalization procedure as explained in Section 4.3. This yields four pairs of Cartesian canonical variables $(\xi, \eta)$ which replace the $\tilde{x}: \tilde{x} = T(\xi, \eta)$, with $T$ being the diagonalizing matrix. Secondly, as in Section 4.3, we introduce canonical polar coordinates $(I_l, \phi_l)_{l=1, \dots, 4}$ by $(\xi_l = \sqrt{2}I_l \cos \phi_l, \eta_l = \sqrt{2}I_l \sin \phi_l)$. The Hamiltonian $\tilde{H}$ will then depend on the variables $(I_1, \dots, I_4, \delta\Delta\lambda_{1,2}', \phi_1, \dots, \phi_4, \delta\lambda_{1,2}')$. Thirdly, we write $\tilde{H}$ in the variables (44); $\tilde{H}$ contains a term in $\Theta^2$ independent of the angles, but its coefficient is not $\beta$. To obtain the value of $\beta$ we need to perform a fourth step, and calculate $(\tilde{H})_{(\phi^*, \delta\lambda_{1,2}^*)}$, which is the perturbative elimination in $\tilde{H}$ of all the non-secondary-resonant contributions from the angles $(\phi^*, \delta\lambda_{1,2}^*)$, up to order 2 in $\Theta$. This is because, as detailed below, the elimination of these harmonics can generate terms in $\Theta^2$ independent of the angles, which need to be added to the original term to obtain $\beta$. To this end, we take $\tilde{H}(I_1, \dots, I_4, \delta\Delta\lambda_{1,2}', \phi_1, \dots, \phi_4, \delta\lambda_{1,2}')$ and expand it to order 2 with respect to the actions. Since these terms satisfy the d'Alembert

characteristics in $(I, \phi)$, we only obtain terms like

$$c_{\alpha}(\delta\Delta\lambda'_{1,2}) \times \sqrt{2I}^{\alpha} \cos(m\phi), \\ m_j = -\alpha_j, -\alpha_j + 2, \dots, \alpha_j - 2, \alpha_j, |\alpha| = 1, 2, 3, 4, \quad (46)$$

where $\alpha \in \mathbb{N}_0^4$, $|\alpha| = \alpha_1 + \dots + \alpha_4$ is restricted to $|\alpha|/2 \le 2$, and $c_{\alpha}(\delta\Delta\lambda'_{1,2})$ is a coefficient that depends on $\delta\Delta\lambda'_{1,2}$ only. These coefficients are expanded around $\delta\Delta\lambda'_{1,2} = 0$ to an optimal order which can be obtained in the following manner. Note that, from $\delta\Delta\lambda'_{1,2} = \delta\Delta\lambda_{1,2}^* + \Theta/6$ (equation 44) each term of order $d$ in $\delta\Delta\lambda_{1,2}'$ contributes a term of order $d$ in $\Theta$, so for each term of order $|\alpha|/2$ in $I$, we must obtain $c_{\alpha}(\delta\Delta\lambda'_{1,2})$ only up to order $\lfloor 2 - |\alpha|/2 \rfloor$ in $\delta\Delta\lambda_{1,2}'$ (where $\lfloor \cdot \rfloor$ is the floor function) to achieve the desired second order with respect to all the actions. We can then organize all terms with respect to the order of expansion in $I$ and $\delta\Delta\lambda_{1,2}',$ and write for each addend $\tilde{H}_{s/2}^j = \mathcal{O}(I^{s/2}) \times \mathcal{O}(\delta\Delta\lambda_{1,2}'^j)$. To get a sense of what these terms look like, we write the terms only up to order $s = 2$ in $\sqrt{I}$:

$$\tilde{H}_{1/2}^0 \propto \sqrt{2I} \cos(\phi_l), \quad \tilde{H}_{1/2}^1 \propto \delta\Delta\lambda_{1,2}' \sqrt{2I} \cos(\phi_l), \\ \tilde{H}_1^0 \propto I_l, \quad \tilde{H}_1^1 \propto \delta\Delta\lambda_{1,2}' \sqrt{2I_1} \sqrt{2I_2} \cos(\phi_{l_1} \pm \phi_{l_2}). \quad (47)$$

The subsequent terms of higher order in $\sqrt{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 $\propto I_l^2$ appearing without angles, as well as the term $\propto \delta\Delta\lambda_{1,2}'^2$, which contribute directly to the $\Theta^2$ term in (45). We note that the first terms $\tilde{H}_{1/2}^0 \propto \sqrt{2I} \cos(\phi_l)$ (corresponding to $|\alpha| = 1$ and constant in $\delta\Delta\lambda_{1,2}'$) are actually zero by the definition of equilibrium point calculated at the reference value $\delta\Delta\lambda_{1,2}' = 0$. Also, the coefficients in $\tilde{H}_1^0$ in front of the $I_s$ are just the frequencies $\omega_l$, since these are the ones calculated in (22). Therefore, in this case, the role of the integrable part of the Hamiltonian for a perturbation theory step is naturally played by $\tilde{H}_1^0 = \sum_{l=1}^4 \omega_l I_l$.

To understand how the perturbative elimination of the non-resonant harmonics involving $(\phi^*, \delta\lambda_{1,2}^*)$ can generate terms in $\Theta^2$ independent of the angles, consider that if $f_n = \mathcal{O}(I^n)$ and $\chi_m = \mathcal{O}(I^m)$, then $\{f_n, \chi_m\} = \mathcal{O}(I^{n+m-1})$, and so on for all the other terms of the Lie series: $\{\{f_n, \chi_m\}, \chi_m\} = \mathcal{O}(I^{n+2m-2})$, etc. When we eliminate a $\tilde{H}_{n,\text{pert}}^j = \mathcal{O}(I^j)$ term by solving the homological equation $\{\tilde{H}_1^0, \chi_m\} + \tilde{H}_{n,\text{pert}}^j = 0$, we must naturally use a $\chi_n = \mathcal{O}(I^n)$. This introduces new terms $\{\tilde{H}_1^0, \chi_n\}, \chi_n$ and $\{\tilde{H}_{n,\text{pert}}^j, \chi_n\}$, which are $\mathcal{O}(I^{2n-1})$ (given that $\tilde{H}_1^0 = \mathcal{O}(I)$). Thus, terms of order 2 can be
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Evolution around $m_{pl}/M_* \simeq 1.28 \times 10^{-3}$ of the angle $\theta$, the resonant angle for the secondary resonance encountered (we plot $2\theta$ instead of $\theta$, as explained in the main text), in panel (a). This figure should be compared to Fig. 6: the action conjugated to $\theta$ is $\Theta = I_2$ (equation 44), and we see that when $\theta$ start librating, $I_2$ increases, indicating that the system has captured into this secondary resonance. As in Fig. 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 Fig. 8, where we take snapshots of the evolution of the pair $(\Theta, \theta)$ at different values of $m_{pl}$.</figcaption></figure>

generated, for example, if $n = 3/2$. We actually need to calculate explicitly only those terms that yield a $\Theta^2$ independent of the angles, as the others would be eliminated further. Such terms derive from $\tilde{H}_{3/2}^0$ (which governs the non-linear interactions between the four resonant degrees of freedom around the equilibrium point that Fig. 6 proves to be strong) and $\tilde{H}_{1/2}^1$ (which describes the fact that the equilibrium point $\tilde{p}_{eq}$ of $\tilde{H}$ shifts as $\delta\Delta\lambda'_{1,2}$ changes under the effects of the $O(m_{pl}^2)$ terms).

We implemented this procedure with the aid of the algebraic manipulator MATHEMATICA. In our reference case $k = 3$ and $a_1 \simeq 0.1$ at a mass $m_{pl}/M_* \simeq 1.28 \times 10^{-3}$, right before the development of the excitation of the resonant degrees of freedom (Fig. 6), this yields

$$ \langle \tilde{H} \rangle_{(\phi^*, \delta\lambda'_{1,2}^*)}(\Theta) = \delta\Theta + \frac{\beta}{2}\Theta^2, \quad \delta \simeq 7.74 \times 10^{-3}, \quad \beta \simeq 101. \qquad (48) $$

The fact that $\delta$ is positive and small is consistent with the fact that we put ourselves right before the development of the excitation (cf. Fig. 4). The fact that $\beta \sim 100$ is positive yields an analytical confirmation that there can be capture into this secondary resonance.

After we have obtained $\langle \tilde{H} \rangle_{(\phi^*, \delta\lambda'_{1,2}^*)}$ (cf. equation 48), we can easily complete the determination of the model (45) for this secondary resonance. To do this, with the help of the algebraic manipulator MATHEMATICA, we use the canonical transformation (44) applied to the term (42b) that contains the resonant harmonic $2\theta$, and we obtain

$$ \langle \mathcal{H}' \rangle_{(\phi^*, \delta\lambda'_{1,2}^*)}(\Theta, \theta) = \delta\Theta + \frac{\beta}{2}\Theta^2 + c(\sqrt{2\Theta})^2 \cos(2\theta + 2\pi/6), \\ \delta \simeq 7.74 \times 10^{-3}, \quad \beta \simeq 101, \quad c \simeq -7.8 \times 10^{-4}. \qquad (49) $$

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 $\mathcal{H}^* = \mathcal{H}_{\text{kepl}} + \mathcal{H}_{\text{res}} + \mathcal{H}_{\text{syn}}$. The evolution of the action $\Theta = I_2$ is already shown in Fig. 6. We plot in Fig. 7 the evolution of the angle $2\theta$ (which produces a numerical evolution that is graphically more legible than that of $\theta$). One can see that the angle starts librating at the same value of $m_{pl}$ where the conjugated

action $\Theta = I_2$ starts increasing in Fig. 6: this shows that there is a passage across the resonance. Note that in such dynamic, 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 $\delta$ in (49) and look at the corresponding phase diagrams. Notice that changing $\delta$ essentially corresponds to changing $m_{pl}$; we also checked that for different planetary masses near $m_{pl}/M_* \simeq 1.28 \times 10^{-3}$ the coefficients $\beta$ and $c$ do not change considerably, so we keep them fixed to obtain a qualitatively correct description of the dynamical portraits.

Fig. 8 shows the level plots of the Hamiltonian (49), for different values of $\delta$ (i.e. of the frequency of $\delta\lambda'_{1,2}/3 + 2\phi_2$ at $\Theta = 0$), in the variables $(X = \sqrt{2\Theta}\cos(2\theta), Y = \sqrt{2\Theta}\sin(2\theta))$; we also overplot the evolution of $(X, Y)$ obtained from the simulation (a combination of Figs 6 and 7), truncated at the value of the planetary mass corresponding to the same $\delta$ 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-hand 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 $\delta$, 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 the angle $2\theta$ 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 when $\delta\beta < 0$ and when the two bifurcations occur at close values of $\delta$. However, in this specific case, we checked that this evolution is actually the result of more complicated interactions among the variables $(I_i, \phi_i)$ themselves, as well as another secondary resonance involving the variables $(I_i, \phi_i)$ and $(\delta\Delta\lambda'_{1,2}, \delta\lambda'_{1,2})$. First, from Fig. 6, one can see that after the initial increase of $I_2, I_4$ also starts increasing, after which there are wide oscillations of $I_2$ and $I_4$ in opposite phase. This is symptomatic of the effect of the term $\sqrt{I_2}I_4\cos(\phi_2 - 2\phi_4)$, which is quasi-resonant because $\omega_2 \simeq 2\omega_4$ (Fig. 2c). To prove this, we plot in Fig. 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 \times 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 $\delta\lambda'_{1,2}/3 + \phi_2 + 2\phi_4$ (which also has a small frequency, since $\delta\lambda'_{1,2}/3 + 2\phi_2$ and $-\phi_2 + 2\phi_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 $\delta\lambda'_{1,2}/3 + \phi_2 + 2\phi_4$ through the canonical change of variables :

$$ \begin{aligned} & (I_4 + 2I_2)/4, \qquad \delta\lambda'_{1,2}/3 + \phi_2 + 2\phi_4, \\ & (I_4 - 2I_2)/4, \qquad -\delta\lambda'_{1,2}/3 - \phi_2 + 2\phi_4, \\ & \delta\Delta\lambda'_{1,2} - I_2/3, \qquad \delta\lambda'_{1,2}. \end{aligned} \qquad (50) $$

We see that after the angle $\theta = \delta\lambda'_{1,2}/6 + \phi_2$ leaves the first resonance (at mass $m_{pl}/M_* \simeq 1.29 \times 10^{-3}$; see Fig. 7), the action $2I_2 + I_4$ keeps growing, which indicates a transition to this new resonance involving $\delta\lambda'_{1,2}/3 + \phi_2 + 2\phi_4$. This analysis
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Contour plots of the Hamiltonian of the integrable model (49) for the secondary resonance involving $2\theta = \delta\lambda'_{1,2}/3 + 2\phi_2$ and $\Theta = I_2$, shown in panels (a)–(f) in canonical Cartesian coordinates at different values of the parameter $\delta$. The change in $\delta$ represents the change in the planetary mass $m_{\text{pl}}$ implemented along the integration. The dots represent the evolution of the system along the numerical integration, and their colour indicates the value of $m_{\text{pl}}$ using the same colour scheme as in Fig. 7.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Evolution of the action $I_4 + 2I_2$, which is a constant of motion relative to the harmonic term $\phi_2 - 2\phi_4$, as well as the resonant action conjugated to the slow angle $\delta\lambda'_{1,2}/3 + \phi_2 + 2\phi_4$ (cf. equation 50).</figcaption></figure>

shows that the evolution presented above is very rich, and does not allow any simple description of it. In any case, Fig. 8 does not leave any doubt that the initial growth of $I_2$ is due to the interaction with the secondary resonance associated with the angle $\delta\lambda'_{1,2}/3 + 2\phi_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 $\mathcal{H}^* = \mathcal{H}_{\text{kepl}} + \mathcal{H}_{\text{res}} + \mathcal{H}_{\text{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 $\mathcal{H}^*$ system, which involve a fraction of the synodic frequency and combination of the resonant frequencies, and which appear at order 2 in the planetary mass. Then, we found the specific secondary resonance that is encountered in the numerical integration of $\mathcal{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, followed by a rich and fascinating series of interactions with additional secondary resonances.

## 5 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 generalized 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 $\mathcal{H}^* = \mathcal{H}_{\text{kepl}} + \mathcal{H}_{\text{res}} + \mathcal{H}_{\text{syn}}$, where $\mathcal{H}_{\text{res}}$ contains the resonant interactions between all N − 1 pairs of neighbouring planets, and $\mathcal{H}_{\text{syn}}$ contains terms of type $\cos(\lambda_i - \lambda_{i+1})$, i = 1, ..., N − 1; for both $\mathcal{H}_{\text{res}}$ and $\mathcal{H}_{\text{syn}}$, we will consider interaction terms up to order 1 in the eccentricities, as we did in the previous section. We then make use of a canonical transformation analogous to (14, 15): We introduce the resonant angles $\psi_1^{(i)} = k\lambda_{i+1} - (k-1)\lambda_i + \gamma_i$ and the angles $\delta\gamma_{i,i+1} = \gamma_i - \gamma_{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 $\delta\lambda_{1,2} = \lambda_1 - \lambda_2$, and, finally, a $\gamma'_N = -\gamma_N$, which will not appear in the Hamiltonian (its conjugated action will again be the total orbital angular momentum).

We can immediately generalize the result of Section 4.3 and state that the $2(N-1)$ purely resonant degrees of freedom are Lyapunov-stable at a 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 (10). So, each planet is either the inner or outer planet, or a middle planet as in the case already considered of a three-planets system. Therefore, all resonant libration frequencies will always have the same sign, and the reasoning of Section 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 (39) and (42) can be encountered. To answer this question, we proceeded analytically following the steps of Section 4.4.2. We introduce a generating Hamiltonian $\chi_{\text{syn}}$, which eliminates the synodic contribution $H_{\text{syn}}$, so $\chi_{\text{syn}}$ in Delaunay variables will have harmonic terms $\sin(\lambda_i - \lambda_{i+1})$, $i = 1, \dots, N-1$. Transforming $H^* = H_{\text{kepl}} + H_{\text{res}} + H_{\text{syn}}$ with the Lie series generated by $\chi_{\text{syn}}$ eliminates $H_{\text{syn}}$ to first order in $m_{\text{pl}}$ (the planetary mass for all planets), but introduces new terms to order 2 in $m_{\text{pl}}$ (among which the most important is $\{H_{\text{kepl}}, \chi_{\text{syn}}\}, \chi_{\text{syn}}\}$, like in the case $N=3$); these newly introduced terms will contain a fraction of the synodic angle $\delta\lambda_{1,2}$, and we are interested in the smallest fraction of $\delta\lambda_{1,2}$ that appears. Like in the case of three planets, the term $(m_{\text{pl}}^2/2)\{H_{\text{kepl}}, \chi_{\text{syn}}\}, \chi_{\text{syn}}\}$ combines together all synodic angles $\lambda_i - \lambda_{i+1}$. Notice that in the new coordinates $\psi_1^{(j)}$, $\delta\gamma_{i,i+1}$, and $\delta\lambda_{1,2}$, each $\lambda_i - \lambda_{i+1}$ can be written as $\lambda_i - \lambda_{i+1} = (\frac{k-1}{k})^{i-1}(\lambda_1 - \lambda_2) = (\frac{k-1}{k})^{i-1}\delta\lambda_{1,2}$ plus terms including $\psi_1^{(j)}$ and $\delta\gamma_{j,j+1}$. However, we do not need to keep track of the $\psi_1^{(j)}$s and $\delta\gamma_{j,j+1}^s$ since we are only interested in the way the angle $\delta\lambda_{1,2}$ appears in the $O(m_{\text{pl}}^2)$ terms. The smallest fraction of $\delta\lambda_{1,2}$ will be generated by combining the synodic angles relative to the two outermost pairs $\lambda_{N-2} - \lambda_{N-1}$ and $\lambda_{N-1} - \lambda_N$, since already they contain the smallest fraction of $\delta\lambda_{1,2}$. Multiplying them together (using $\cos(a)\cos(b) = \frac{1}{2}(\cos(a-b) + \cos(a+b))$) yields a harmonic term of type

$$ \cos\left(\left(\left(\frac{k-1}{k}\right)^{N-3} - \left(\frac{k-1}{k}\right)^{N-2}\right)\delta\lambda_{1,2} + \dots\right), \quad (51) $$

where the $+\dots$ terms represents a combination of $\psi_1^{(j)}$s and $\delta\gamma_{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

$$ \frac{1}{k} \left( \frac{k-1}{k} \right)^{N-3} \delta\lambda_{1,2} \approx \frac{1}{k^2} \left( \frac{k-1}{k} \right)^{N-3} n_1. \quad (52) $$

This is the fraction of the synodic frequency which can resonate with the libration frequencies $\omega_l$ of the resonant degrees of freedom. Since $\omega_l$ increase with $m_{\text{pl}}$ (as $\omega \sim m_{\text{pl}}^{1/2}$ or $m_{\text{pl}}^{2/3}$ 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 $\frac{1}{k}(\frac{k-1}{k})^{N-3}$ multiplying $\delta\lambda_{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_{\text{pl}}/M_*)_{\text{crit}}$ allowed for stability decreases with increasing $N$ and $k$. This gives an analytical explanation to the numerical findings of Matsumoto et al. (2012).

## 6 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 per cent (Izidoro et al. 2017, 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 a 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 that 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 that 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 generalized 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 destabilized. This scheme can then be generalized 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$ (equation 52). 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) and 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*<sub>pl</sub>, and produce an explicit criterion for the stability against secondary resonances of the type described here. This will be the subject of future work.

## ACKNOWLEDGEMENTS

G. P. thanks the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support.

## REFERENCES

* Batygin K., 2015, MNRAS, 451, 2589

* Batygin K., Adams F. C., 2017, AJ, 153, 120

* Batygin K., Morbidelli A., 2013b, A&A, 556, A28

* Chatterjee S., Krantzler S. O., Ford E. B., 2016, Proc. IAU Symp. 29A, Astronomy in Focus. Kluwer, Dordrecht, p. 30

* Cresswell P., Nelson R. P., 2008, A&A, 482, 677

* Deck K. M., Batygin K., 2015, ApJ, 810, 119

* Fressin F. et al., 2013, ApJ, 766, 81

* Gillon M. et al., 2016, Nature, 533, 221

* Gillon M. et al., 2017, Nature, 542, 456

* Gladman B., 1993, Icarus, 106, 247

* Howard A. W. et al., 2012, ApJS, 201, 15

* Izidoro A., Ogihara M., Raymond S. N., Morbidelli A., Pierens A., Bitsch B., Cossou C., Hersant F., 2017, MNRAS, 470, 1750

* Izidoro A., Bitsch B., Raymond S. N., Johansen A., Morbidelli A., Lambrechts M., Jacobson S. A., 2019, Formation of planetary systems by pebble accretion and migration: Hot super-Earth systems from breaking compact resonant chains, preprint (arXiv:1902.08772)

* Laskar J., 1990, Icarus, 88, 266

* Laskar J., Robutel P., 1995, Celest. Mech. Dyn. Astron., 62, 193

* Luger R. et al., 2017, Nat Astron., 1, 0129

* Marchal C., Bozis G., 1982, Celest. Mech., 26, 311

* Masset F. S., Morbidelli A., Crida A., Ferreira J., 2006, ApJ, 642, 478

* Matsumoto Y., Nagasawa M., Ida S., 2012, Icarus, 221, 624

* Mayor M. et al., 2011, The HARPS search for southern extra-solar planets XXXIV. Occurrence, mass distribution and orbital properties of super-Earths and Neptune-mass planets, preprint (arXiv:1109.2497)

* Michtchenko T. A., Beaugé C., Ferraz-Mello S., 2008, MNRAS, 387, 747

* Millholland S., Wang S., Laughlin G., 2017, ApJ, 849, L33

* Mills S. M., Fabrycky D. C., Migaszewski C., Ford E. B., Petigura E., Isaacson H., 2016, Nature, 533, 509

* Morbidelli A., 2002, Modern celestial mechanics: aspects of solar system dynamics, Taylor & Francis, London

* Morbidelli A., Crida A., Masset F., Nelson R. P., 2008, A&A, 478, 929

* Murray C. D., Dermott S. F., 1999, Solar system dynamics, Cambridge University Press, Cambridge, UK

* Obertas A., Van Laerhoven C., Tamayo D., 2017, Icarus, 293, 52

* Ogihara M., Morbidelli A., Guillot T., 2015, A&A, 578, A36

* Petigura E. A., Howard A. W., Marcy G. W., 2013, Proc. Natl. Acad. Sci., 110, 19273

* Petit A. C., Laskar J., Boué G., 2018, A&A, 617, A93

* Pichierri G., Morbidelli A., Crida A., 2018, Celest. Mech. Dyn. Astron., 130, 54

* Poincaré H., 1892, Les methodes nouvelles de la mecanique celeste, Gauthier-Villars et fils, Paris

* Ramos X. S., Charalambous C., Benítez-Llambay P., Beaugé C., 2017, A&A, 602, A101

* Rogers L. A., 2015, ApJ, 801, 41

* Terquem C., Papaloizou J. C. B., 2007, ApJ, 654, 1110

* Weiss L. M. et al., 2018, AJ, 155, 48

* Winn J. N., Fabrycky D. C., 2015, ARA&A, 53, 409

* Xu W., Lai D., Morbidelli A., 2018, MNRAS, 481, 1538

* Zeng L. et al., 2019, Proc. Natl. Acad. Sci., 116, 9723

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

**APPENDIX A:** Simulating capture into mean motion resonance via convergent migration

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

This paper has been typeset from a *TEX/LATEX* file prepared by the author.