An integrable model for first-order three-planet mean motion resonances

Recent works on three-planet mean motion resonances (MMRs) have highlighted their importance for understanding the details of the dynamics of planet formation and evolution. While the dynamics of two-planet MMRs are well understood and approximately described by a one degree of freedom Hamiltonian, little is known of the exact dynamics of three-bodies resonances besides the cases of zeroth-order MMRs or when one of the body is a test particle. In this work, I propose the first general integrable model for first-order three-planet mean motion resonances. I show that one can generalize the strategy proposed in the two-planet case to obtain a one degree of freedom Hamiltonian. The dynamics of these resonances are governed by the second fundamental model of resonance. The model is valid for any mass ratio between the planets and for every first-order resonance. I show the agreement of the analytical model with numerical simulations. As examples of application I show how this model could improve our understanding of the capture into MMRs as well as their role on the stability of planetary systems.

The dynamics of two-planet MMRs are now well understood. The dynamics of first-order MMRs are given by the so-called second fundamental model of resonance (Henrard and Lemaitre 1983). The resonant dynamics of two massive planets on eccentric orbits are reduced to this one degree of freedom integrable Hamiltonian thanks to a sequence of transformations first introduced by Sessin and Ferraz-Mello (1984) and Henrard et al. (1986). A generalization of this model for two-planet MMRs of arbitrary order has been proposed by Hadden (2019).
On the other hand, the study of three-bodies resonances has long been restricted to the asteroids dynamics in the Solar System (Nesvorný and Morbidelli 1998;Cachucho et al. 2010). The question has gained a renewed interest thanks to the discovery of close-in exoplanet systems trapped into zeroth-order order resonances such as Trappist-1 (Gillon et al. 2017;Agol et al. 2021) or more recently TOI-178 (Leleu et al. 2021). Three-planet resonances are also thought to be the main driver of the early instability of tightly packed systems (Quillen 2011;Petit et al. 2020).
The thorough numerical studies by Guzzo (2005) or Charalambous et al. (2018) has shown the key role of three-planet resonances in shaping the phase space. Particularly, Charalambous et al. (2018) highlighted their role during resonance capture, where first-order three-planet MMRs act as a guide towards two-planet MMR during multiple planet migration. Semi-analytical studies (Gallardo et al. 2016) have also been performed and given a qualitative understanding of the threeplanet resonance dynamics. A general integrable model for zeroth-order three-planet MMRs is directly obtained from a second-order averaging since no transformation of the coordinates is needed besides the introduction of the resonant angle (see Petit et al. 2020). However, this result is not immediately generalized to higher order MMRs. The question of the existence of an integrable model for first-order three-planet MMRs is thus of theoretical interest. Besides that, such a model could help better understand the capture into these resonances as well as their role into the instability of tightly packed systems with eccentric orbits. To this day, there is no direct evidence of exoplanets trapped into first-order three-planet MMR. Yet it seems that they could help stabilize some systems as suggested for Kepler-11 (Migaszewski et al. 2012).
In section 2, I introduce the different networks of first-order three-planet resonances and discuss their relative importance with respect to other two-planet and three-planet MMRs. I show in section 3 how to obtain a one degree of freedom Hamiltonian for the dynamics of an isolated first-order three-planet MMR, adapting the results from the two-planet case. I then analyse the properties of the resonance thanks to this integrable model (section 4). Finally, I present in section 5 some examples of application of the model to the problem of stability of eccentric tightly packed systems and the capture into these resonances, while leaving the detailed studies for future work.

Resonance networks
First-order three-planet resonances occur when the mean motions verify the relationship k 1 n 1 + k 2 n 2 + k 3 n 3 = 0, where k 1 +k 2 +k 3 = 1 and n j are the Keplerian mean motions. The resonance equation defines a plane in the frequency space (n 1 , n 2 , n 3 ). Because the gravitational interactions are scale invariant, we can restrict ourselves to a two-dimensional plane corresponding to the period ratios ν 12 and ν 23 where The resonance relation can be rewritten to remove k 2 and use period ratios in place of the mean motions. In the general case, one obtains which is the equation of the locus of the unperturbed resonance defined by k 1 and k 3 . One can remark that except for −k −1 3 on the right-hand side, it is similar to the zeroth-order resonance relationship used in (Petit et al. 2020). However, unlike the zeroth-order case, there are three different families of resonances depending on the signs of k 1 and k 3 . Noting p = |k 1 | and q = |k 3 |, we have Note that the configuration k 1 < 0, k 3 > 0 is not a solution since the period ratios are smaller than 1. We plot the loci of the first-order three-planet resonances on Figure 1 for p + q ≤ 12. The colors are associated with the three different families, blue curves correspond to (4), red curves to (5) and green curves to (6). The blue and red networks are parallel to the zeroth order network while the green network is restricted to the region where the planet pairs are tight. Indeed, all the green resonances with a given k = p + q, are outside the square (0, k/(k + 1)) × (0, (k − 1)/k)), which means that for large indexes they become irrelevant. Nevertheless, one should also note that for p = q + 1, the corresponding green resonance is actually a two planet MMR between planet 1 and 3. An important difference with the case of zeroth order resonances treated in (Petit et al. 2020), is that the resonance locus equation depends specifically on the two integers k 1 and k 3 rather than on their ratio. This property forbids to create a continuous coordinate constant on the resonance loci. As a result, some simplifications possible in the zeroth-order case are no longer possible.
Despite their small width, we can observe the resonant network onto dynamical maps. Similarly to Gallardo et al. (2016) or Charalambous et al. (2018), I perform a series of simulations of three equal masses planet systems. The initial period ratios ν 12 and ν 23 span (0.55, 0.76) forming a grid of 800 × 800 systems. Planet  Fig. 1 Loci of the first-order three planet resonances with p + q ≤ 12. The colors correspond to the three different families of resonances defined in Eqs. (4, blue), (5, red) and (6, green). The dashed black lines correspond to first-order two-planet MMRs (obliques lines correspond to resonances between the the first and third planets). One can note that for non tight systems, the majority of the resonances belong to the "red" family (eq. 5).
masses are 10 −4 M , the orbits are coplanar and initially circular. The innermost planet semi-major axis is 1 au, the star mass 1M and the starting mean longitudes are random. The masses and the period ratios correspond to typical spacings for exoplanet multiple systems. The simulations are integrated for 10 4 yrs using the integrator WHFast (Rein and Tamayo 2015) from the REBOUND library (Rein and Liu 2012). I measure the chaotic nature of the dynamics thanks to the Mean Exponential Growth Factor of Nearby Orbits (MEGNO, Cincotta et al. 2003). The MEGNO oscillates around 2 for quasi-periodic orbits and grow to infinity for chaotic orbits.
The dynamic map is displayed on Figure 2. One can spot the two-planet MMRs 3:2 and 4:3, responsible for the largest chaotic regions. The largest zeroth-order three-planet resonances, are highlighted as dashed black lines. Beyond these features, most of the obliques lines showing weak chaos correspond to first-order three-planet resonances.  First-order three-planet MMRs appears as thin moderately chaotic lines as described by the legend. Zeroth-order threeplanet MMRs are larger but harder to spot, hence dashed black lines were added to highlight them.
3 Integrable first-order three-planet MMRs Hamiltonian 3.1 Second-order averaging We consider a system of three planets of masses m 1 , m 2 , and m 3 orbiting a star of mass m 0 . The canonical positions r j and momentar j are expressed in canonical heliocentric coordinates (Poincaré 1905;Laskar 1991). The orbits are assumed to be close to circular and coplanar. Let the semi-major axes a j , the eccentricities e j , the mean longitudes λ j , and the periapses longitude j be the orbital elements defining the orbits. A set of canonical coordinates for the system is given by the modified Delaunay coordinates (e.g. Laskar 1991): where µ = Gm 0 and G is the gravitational constant. We note that the gravitational parameter µ is the same for all three planets as in Laskar and Petit (2017). This is possible by considering the so-called democratic-heliocentric formulation of the planetary Hamiltonian (e.g. Morbidelli 2002). The couples of variables (C j , − j ) can also be replaced by their associated complex variables 1 For small eccentricities, we have x j Λ j /2e j e ι j . The system total angular momentum G and Angular Momentum Deficit (AMD Laskar 1997) C are given by The Hamiltonian H describing the dynamics can be split into an integrable part describing the motion on unperturbed Keplerian orbits, and a perturbation, describing the planet interactions. Here, ε is a dimensionless parameter of the order of the planet-to star-mass-ratio to reflect the scale difference between the two parts of the Hamiltonian. I limit my study to term at first-order in eccentricity and inclination. This assumption justifies only treating the planar case since inclination terms only appear at order two in the inclinations. The relevant terms of the perturbing Hamiltonian εH 1 take the form (Laskar and Robutel 1995;Murray and Dermott 1999)  where, in the last three sums, c.c. designates the complex conjugate of the term, and s (α) are the Laplace coefficients (see e.g. Laskar and Robutel 1995). S l ij and U j ij are the terms at second order in eccentricities and the coefficients f k are taken from Murray and Dermott (1999) and defined explicitly in Appendix A. I emphasize that while the dependency in x j is written explictly for first-order terms, I choose to keep it implicit in the second order terms as this helps keep the computations clearer. The inclusion of the second order terms is necessary since they create linear terms in eccentricity during the averaging process described below in Section 3.3. The last sum in (12) corresponds to the indirect terms due to the star's reflex motion. However, as in the two planet case, at first-order in eccentricity, the indirect terms are only relevant for resonances where p = ±1 or q = ±1. As a result, I will not consider the indirect part in the computation of the resonance model 2 .
In the unperturbed case, the system is said to be in a three-body MMR if the mean motions verify an equation of the form k 1 n 1 + k 2 n 2 + k 3 n 3 = 0. The sum k = k 1 + k 2 + k 3 is the 'order' of the resonance. The sum K = |k 1 | + |k 2 | + |k 3 | is the 'index' of the resonance. Because of the d'Alembert rules (e.g. Morbidelli 2002), the leading order term in the perturbation is of order k in eccentricity.
Three-planet resonances emerge in the perturbative Hamiltonian after canonical transformations that aim at removing the dependency at first-order in ε of εH 1 into the fast, non-resonant angles λ j . Such averaging is possible if the system is situated 'far enough' from any two-planet MMR using the classical approach from perturbation theory, the Lie series method (Deprit 1969). In order to perform the averaging and compute the resonant coefficients, I follow the same method and notations as in (Petit et al. 2020) and I refer the reader to this work for a more extensive description. Nevertheless is it convenient to recall the notations used. I note εχ 1 the generating Hamiltonian, solution of the homological equation where εH 1 is the average of εH 1 over the mean longitudes λ and {·, ·} is the Poisson bracket 3 . Noting h , the complex Fourier coefficients of H 1 with respect to the mean longitudes, εχ 1 has for expression Due to the expression of εH 1 given in eq. (12), the denominators k · n are of the form k j n j + k j n j (in particular they only include terms relative to single pairs of planets) and are not 'too small' because I assume the system to be far from two-planet MMRs. Thus the formal series (20) is well defined; one can stop the summation at indices k of sufficiently high order so that the remaining Fourier terms in εH 1 have sizes smaller than ε 2 , which is ensured by the exponential decay of the Fourier coefficients. The averaged coodinates (noted with a superscript 1) are obtained by applying the transformation defined by the flow at time -1 of the Hamiltonian function εχ 1 to the osculating coordinates The new Hamiltonian is where the second-order term in ε can be expressed as ε 2 H 2 still contains terms depending on fast angles. The study of a particular threeplanet MMR can be done by a second averaging over the non-resonant angles. In theory, this results in another change of coordinates, which are ε 2 close to the first-order averaged coordinates. In practice I drop these ε 2 -order corrections to the averaged variables. I also drop the terms of order ε 3 and greater. To keep the notations light, from now on, unless specified differently I will only consider the averaged coordinates and drop the superscripts.

Change of variables to resonant angles
The three families of resonances Eqs. (4-6) can be treated within the same framework. As in the two-planet case, there are three resonant angles (one per planet) associated with a first-order three-planet resonance. Indeed, due to d'Alembert rules, the terms in the Hamiltonian must depend on an angle of the form ϕ j = k 1 λ 1 +k 2 λ 2 +k 3 λ 3 − j , where j = 1, 2, 3. Following the same strategy as in (Delisle Petit et al. 2017), I make a linear, canonical change of coordinates that keeps the highest degree of symmetry. The new set of angles is which give for actions where C is the total AMD. G is the system total angular momentum and by analogy with the two-planet case (Michtchenko et al. 2008), Γ 1 and Γ 3 are two spacing parameters. I also define the set of complex coordinates (−ιx j , x j ) associated with The angles θ 1 , θ 3 and θ G are fast angles for the resonance motion. The secondorder formal averaging over the fast angles removes them form the Hamiltonian that thus only depend the angles ϕ j . As a result, Γ 1 , Γ 3 and G are constant of the motion for a particular resonance defined by k 1 , k 2 and k 3 . Inverting the action coordinates transformation, we have One can note that the variations of the semi-major axes only depend on the variation of the total AMD and not on the individual planet AMDs. This property is similar to the two-planet case (Delisle et al. 2014;Petit et al. 2017). In order to make progress, one needs to develop the actions close to the resonant manifold. In the standard derivation of the integrable model for two planet in MMR (e.g. Deck et al. 2013;Delisle et al. 2014), the action variables are renormalized by the scaling parameter. This reflects the fact that the dynamics can be rescaled by the orbital timescale of the system. One can then write the actions Λ j = Λ j,0 + dΛ j such that the constant values (Λ 1,0 , Λ 2,0 , Λ 3,0 ) verifies the Keplerian resonance relationship. In our case, there are two scaling factors describing the spacing of the two pairs of planets. As a result, I skip the renormalization phase but keep the idea to develop around the Keplerian resonance. Indeed, for a given value of the constants of motion Γ 1,0 and Γ 3,0 , there is a unique angular momentum G 0 = Λ 1,0 + Λ 2,0 + Λ 3,0 such that the circular orbit configuration with such semi-major axes verify the resonance relationship 4 . Noting ∆G = G 0 − G, we can write

Computation of the perturbative coefficients
We need to select from Eq. (23), the terms that remain after the averaging of the non resonant angles. Since εH 1 is independent of λ j , the term {εχ, εH 1 } only contains terms depending on two planet mean longitudes and does not contribute to the resonant dynamics. Due to the form of the terms of εH 1 , εH 1 and εχ 1 , there are twelve terms (and their complex conjugates) contributing to the resonant Hamiltonian of a first-order three-planet MMR. For the resonance defined by k 1 and k 3 , four terms are created by the combination of a zeroth order term and a first order term For these terms, the Poisson bracket reduces to ∂ ∂Λ2 ∂ ∂λ2 − ∂ ∂λ2 ∂ ∂Λ2 . Furthermore, a combination of a second order term with a first order term can create a term depending on the resonant angle. There are eight different terms to consider −ιV −k1 12,>x 2 e ι(k1λ1−(k1−1)λ2) ιV k3−1 23,<x 2 e ι(−(k3−1)λ2+k3λ3) −ιS k1 12 e ιk1(λ1−λ2) Since I limit myself to an expansion at first-order in eccentricity, the Poisson bracket reduces to either ι ∂ ∂x2 ∂ ∂x2 or −ι ∂ ∂x2 ∂ ∂x2 . The terms in Eqs. (30) are absolutely necessary as they have a comparable size to the terms defined in Eqs. (29).
As in the case of zeroth-order three-planet MMRs, ones need to compute the terms of Eq. (23) that depend only on the angles ϕ j . The expression of ε 2 H 2,res as a function of V k ij,≶ , W k ij , S k ij and U k ij is given in Appendix C. After regrouping the terms, the perturbation part of the Hamiltonian can be expressed as where 3≶ (α 12 , α 23 ), and ≶ is a shorthand notation for the two symbols < and >. g k1,k3 j≶ are functions given in Appendix C that are linear combinations of Laplace coefficients and their derivatives and depend only on the semi-major axis ratio. Like the Laplace coefficients, these functions show an exponential decay with |k 1 | and |k 3 | as well as with the orbital spacing.

Generalized Sessin-Henrard transformation
In the two-planet case, a final canonical transformation is needed to obtain a one degree of freedom, integrable Hamiltonian. This transformation, introduced by Sessin and Ferraz-Mello (1984) and Henrard et al. (1986), can be interpreted as a two-dimensional rotation of the canonical complex variables x j . Because the perturbation is linear in the x j variables, I can generalize the Sessin-Henrard transformation by introducing a three-dimensional rotation, creating two additional integrals of motion. I define the canonical change of variables 5 where and This transformation is canonical since the rotation matrix is orthogonal. I define the associated action-angle coordinates to the complex coordinates y j as (I j , ψ j ), such that y j = I j e ιψj . The AMD still verifies C = y 1ȳ1 + y 2ȳ2 + y 3ȳ3 = I 1 +I 2 +I 3 . As in the two-planet case, at the first order in eccentricities, the resonant variable is proportional to linear combination of the complex eccentricities with coefficients depending only on the semi-major axis ratios Finally, I expand at the second order in (C − ∆G) the Keplerian part around the resonance for unperturbed orbits and drop the first order secular terms. The Hamiltonian becomes Since H does not depend on ψ 2 and ψ 3 , the actions coordinates I 2 and I 3 are constant of the motion, up to perturbations of order ε 3 .

Structure and width of the resonances
One can recognize in Eq. (38) the second fundamental model of resonance proposed by Henrard and Lemaitre (1983). This Hamiltonian, also called an Andoyer Hamiltonian of degree one appears in the description of first order two planet MMR, Cassini states and various problems of celestial mechanics. The width of the resonance as well as the dynamics of this integrable Hamiltonian are well-known and I refer to Ferraz-Mello (2007) for a detailed description. For the determination of the resonance widths and structure, I will adapt the formalism from Petit et al. (2017).
Outer circulation center X 1 Libration center X 1 Inner circulation center X 2 Unstable point X 3 Separatrices The main difference with the two-planet first-order MMR is that the resonant term is of second order in planet-to-star mass ratio. I define as well as I 0 = (ε 2 κ) −2/3 (∆G − I 2 − I 3 ), X = (ε 2 κ) −1/3 √ 2I 1 cos ψ 1 and Y = (ε 2 κ) −1/3 √ 2I 1 sin ψ 1 . In these renormalized variables, the Hamiltonian can be written The fixed points verify X 3 −2I 0 X −2 = 0 and Y = 0. The cubic equation has three solution if and only if I 0 > 3/2. In this case the resonance is open and separatrices exist. Noting X 1 > X 2 ≥ X 3 the abscissas of the fixed points, X 3 defines the unstable hyperbolic point where the separatrices meet. I also note X * 1 and X * 2 the intersections of the separatrices with the X-axis. Since the separatrices are on the same energy level as the unstable fixed point, we have In Figure 3, I plot the structure of the resonance for I 0 = 3 in the (X, Y ) plane as well as the position of the fixed points and separatrix intersections for varying I 0 .  Figure 2 around the resonance studied in b). The dashed red curve corresponds to the locus of the resonance −3n 1 + 6n 2 − 2n 3 for Keplerian motion and the yellow segment to the extent of the zoomed slice. b) Dynamical map in the plane (ν 12 , e 2 ). The red curves corresponds to the inner and outer separatrices, the orange curve to the center of the resonance.
The width of the resonance can be expressed implicitly as a function of the minimum value of I 1 for a resonant orbit given a value of I 0 . Indeed, the resonance width in terms of I 1 is The value of X 3 is linked the minimum value of I 1 to enter the resonance island thanks to the relationship There are two limit cases of interest, the case of circular orbits and the large eccentricity regime ). In the case of circular orbits, the resonance width becomes (∆I 1 ) cir = 4 × 2 1/6 (ε 2 κ) 2/3 . For eccentric orbits, we have X 3 (ε 2 κ) −1/3 √ 2C min , which leads to (∆I 1 )ecc = 4ε κ √ 2C min . One can then obtain the width of the resonance in term of period ratio by using Eqs (2) and (28) ∆I 1 , and ∆ν 23 = 3ν 23 I also compare this analytical result to a dynamical map at the vicinity of the resonance −3n 1 + 6n 2 − 2n 3 . I choose to keep ν 23 constant while changing ν 12 and the eccentricities of the planets. The one-dimensional cut in the period ratio plane is shown on Fig. 4a, that is a zoom of Fig. 2 around the resonance of interest. I initialize a regular grid of 400 × 400 initial conditions with the averaged variables 0.635 ≤ ν 12 ≤ 0.640 and 0 ≤ e 2 ≤ 0.04. The averaged outer planet period ratio is kept constant ν 23,0 = 0.65, the orbits are coplanar, the averaged eccentricities e 1 , e 3 and the mean longitudes are chosen such that I 2 = I 3 = 0 and arg(y 1 ) = 0. As in Figure 2, the planet masses are 10 −4 M , the systems are integrated 10 4 orbits and I measure the MEGNO along the trajectory. The resonance is small with respect to the rapid variations of the osculating coordinates. Therefore, these averaged initial conditions are converted into osculating orbital elements before the start of the integration.
The dynamical map is shown on Figure 4b. Qualitatively, the shape of the resonance is similar to the structure observed for two-planet MMR, confirming that the second fundamental model of resonance describes well first-order threeplanet MMR. One can see that the outer separatrix does not go to zero eccentricity on the dynamical map, as predicted by the analytical model. Similarly the position of the inner separatrix matches well the dynamical map at low eccentricities.
However, the analytical model slightly underestimates the size of the resonance (by roughly a factor 2 for a mean eccentricity e 2 0.01). A possible explanation could be that higher order terms in eccentricity becomes relevant for non circular orbits and could increase the width of the resonance as discussed by Hadden (2019) for the two-planet MMR case. At higher eccentricities,the width of the resonance is not proportional to √ e 2 , as it would be in the case of the second fundamental model of resonance. Despite this slight discrepancy, the qualitative behaviour predicted by the model is correct.

Role of first-order three-planet resonances on the planet dynamics
This analytical model for first-order three-planet resonances gives an opportunity to study how they affect exoplanet dynamics. We saw that while these resonances are weak with respect to low-order two-planet MMR and zeroth-order three-planets MMRs, they can affect the dynamics locally. In particular the similarities of the dynamics with the two-planet case suggest that their overlap can lead to destabilization of tightly packed systems and that resonance capture is possible in the case of convergent migration, thanks to the open separatrix at low eccentricity. In this section I show two examples of applications of my analytical model to exoplanet dynamics problems. While both problems deserve a detailed study, the full development is left for future works.

Contribution to the instability of tightly packed systems
The motivation for this study emerged while trying to extend the results on the stability of tightly packed systems from Petit et al. (2020) to eccentric systems. Indeed, as the eccentricity increases, the first-order three-planet resonance width increases to the point where they should play a role in driving the instability. The main concept recently developed to treat the question of complex resonance overlap involving varied resonances sizes is the notion of resonance network density or optical depth (Quillen 2011; Hadden and Lithwick 2018;Petit et al. 2020). It consists in comparing the sum of the width of the resonances to the space they occupy. Under the approximation that the resonances are well distributed, when the total volume of the resonances becomes larger than the available space, large scale chaos can develop. As explained in section 2, the main difference between the zeroth-order and first-order resonance networks is that first-order resonances intersect and there is no transformation of coordinates from the period ratios to a coordinate constant onto the resonance loci.
We can nevertheless estimate the density of the first-order three planet resonance network. For a given resonance described by k 1 and k 3 , I define f k1,k3 (n) = k 1 n 1 − (k 1 + k 3 − 1)n 2 + k 3 n 3 .
The resonant coefficient ε 2 κ mostly depends on the sum k 1 + k 3 , it thus makes sense to compute the density of the subnetwork of resonances verifying k 1 + k 3 . We have f k1+1,k3−1 − f k1,k3 = n 1 − n 3 , it results that the subnetwork density is It is in principle possible to obtain an overlap criterion from this expression by summing over the subnetworks. However, considering the first order resonances alone is incorrect and one need to take into account the zeroth-order resonances. Moreover, as the eccentricity increases, more and more two-planet resonances are activated and become the dominant source of chaos. Due to the complexity of this study, this problem will be the topic of further works.

Capture into resonant chains
Disk-driven planet migration leads to capture into two-planet first-order MMRs (e.g. Terquem and Papaloizou 2007;Cresswell and Nelson 2008;Izidoro et al. 2017). Contrarily to the general case (like zeroth-order or second and higher order resonances) where the capture into resonance is essentially probabilistic (Henrard 1982), low-eccentricity convergent migration leads to automatic capture into resonance because no separatrix crossing is needed to enter the resonance (Batygin 2015). The resonance fixed point can however be overstable or the migration can be too fast to allow for capture, depending on the planet masses and migration rate (e.g. Ogihara and Kobayashi 2013;Deck and Batygin 2015).
Since the three-planet first-order resonances dynamics are similar to the twoplanet case, it is natural to consider whether capture into these resonances is possible. Besides being much weaker than their two-planet counterpart, the main difference regarding resonance capture is the fact that the Keplerian resonance locus is a curve in the period ratio plane rather than a fixed value on a onedimensional axis. As a result, after capture, the orbital separation can still tighten, leading to the crossing of other three-planet resonances or even two-planet ones. While this problem has not been studied to the same level of details than the capture into two-planet MMRs, Charalambous et al. (2018) have shown on the example of systems with configuration similar to Trappist-1 that three-planet firstorder resonances can act as guides during migration (their Figures 7 and 8).
In order to illustrate the capture mechanism, I integrate for 1 Myr a three 10 −4 M equal-mass planet system with semi-major axis and eccentricity damping  Fig. 6 Capture into the resonance −3n 1 + 6n 2 − 2n 3 . Left: Time series of the three classical resonant angles ϕ j (24d) and of the argument of the resonant variable ψ 1 = arg(y 1 ). The capture into the resonance occurs around 0.4 Myr as seen by the libration of ψ 1 . After 0.6 Myr, the system is also trapped into the resonance 2n 2 − 7n 2 + 6n 3 (see Figure 7) Right: Evolution of the resonant complex eccentricities as well as the resonant variable y 1 (34) normalized by 2/Λ 2 in order to be comparable with the eccentricities. The moving average is also applied to these variables.
implemented thanks to REBOUNDx (Tamayo et al. 2019). Since the resonances are very weak, it is necessary to adopt very slow damping timescales so that capture is possible. We follow the constant timescale damping prescription from Goldreich and Schlichting (2014). The eccentricity damping timescale is τe = 10 3 yr for all planets. The semi-major axis damping timescales are τ a,1 = +∞, τ a,2 = 20 Myr and τ a,3 = 10 Myr, so that the differential migration rate between the planets is equal to 20 Myr. 3,000 snapshots are recorded and the semi-major axis and complex eccentricities are converted from osculating to averaged variables.
The evolution of the period ratios and eccentricities are plotted in Figure 5. Due to the rapid evolution of the eccentricities within the resonances, a moving average over 20 snapshots was applied. In the period ratio space, we see that the system evolves towards the resonance −3n 1 +6n 2 − 2n 3 and reaches it around t = 0.4 Myr.   Fig. 7 Capture into the resonance 2n 1 − 7n 2 + 6n 3 . Left: Time series of the three classical resonant angles ϕ j (24d) and of the argument of the resonant variable ψ 1 = arg(y 1 ). The capture into this resonance occurs around 0.6 Myr as seen by the libration of ψ 1 . Right: Evolution of the resonant complex eccentricities as well as the resonant variable y 1 (34) normalized by 2/Λ 2 in order to be comparable with the eccentricities. The moving average is also applied to these variables.
Then the system evolves along the resonance locus and stops at 0.6Myr when the resonance −3n 1 + 6n 2 − 2n 3 crosses the resonance 2n 1 − 7n 2 + 6n 3 . The breaks into the evolution can also be seen onto the eccentricities time series. We see that e 1 and e 2 first rise after the first capture before an increase to a constant value after the secondary capture. We can also see that the resonance angles start to librate after the captures in Figures 6 and 7. It should be noted that in Figure 6, ψ 1 librates after the first capture while ϕ 3 is not librating. The capture can also be observed onto the complex eccentricity plane as shown on the right panels of Figures 6 and 7. Since the angles −3λ 1 + 6λ 2 − 2λ 3 − j and 2λ 1 − 7λ 2 + 6λ 3 − j both librate, the zeroth-order angle 5λ 1 − 13λ 2 + 8λ 3 also librates. Indeed, the zeroth-order resonance passes through the intersection point of the two first-order resonances. It is probable that it helps stabilizing the configuration.
This example shows that capture is possible into these resonances. However, the resonance itself does not stop the migration because it is not a hard barrier in the period ratio space. An intersection of resonances can however lead to capture at fixed period ratios. Importantly, these results were obtained in a very controlled environment with ad-hoc migrations timescales, similarly to the experiments of Charalambous et al. (2018). For a faster migration or smaller planets, the resonance width will be too small to allow for capture as the system will swipe through in less than a libration period. The necessary conditions leading to such capture may thus not be present in the context of a realistic formation scenario.

Conclusion
I have shown that the dynamics of all first-order three-planet resonances can be approximated by a novel integrable model using the same strategy than the firstorder two planet MMRs. After a second-order averaging in the planet masses over the non-resonant angles, the dynamics can be reduced to a second fundamental model of resonance (Henrard and Lemaitre 1983) thanks to a generalized Sessin-Henrard transformation of the eccentricity variables (Sessin and Ferraz-Mello 1984;Henrard et al. 1986). My model is the first analytical solution in the presence of three massive planets instead of two-planet and a test particle. From there, I derived the width and the shape of the resonances as a function of the period ratios, planet masses and eccentricities of the planet. In the context of an isolated resonance, the predicted analytical width is within a factor 2 of the observed width on numerical simulations. The analytical model matches really well the width of the resonance for close to circular orbits. I showed that the analytical width could be used in future work in order to predict more accurately the lifetime of systems with moderate eccentricities. The resonance structure of first-order three-planet MMR makes capture into these resonances possible, although it is not clear which regime can most likely lead to it as my results were obtained in a controlled environment.
As shown on Figure 2, isolated first-order three-planet resonances are mostly important just outside from the two-planet and three-planet MMR overlap limits (Hadden and Lithwick 2018;Petit et al. 2020). In the context of exoplanets, this typically corresponds to compact systems with a planet to star mass ratio of the order 10 −4 with period ratios between 1.33 and 1.7. Observationally, it is not possible to confirm directly if a system is trapped into such resonance because the resonant angles depend on the pericentre longitudes that are poorly constrained in general, particularly for close to circular systems. However, in the context of very tight systems, a dynamical study can restrict the stable configurations to the resonant island, suggesting that such resonances exist. In particular, such a scenario was proposed in the context of Kepler 11 (Migaszewski et al. 2012). The possibility that exoplanet systems are trapped into first-order three-planet MMRs deserves more investigation and will be the topic of future works.
(l − 1)λ 1 − lλ 2 + 1 whereas I define such a term as lλ 1 − (l + 1)λ 2 + 1 . As a result I adapt the expressions of f k to reflect this change of convention. I use

B Derivatives of the perturbation terms
As described in the main text, the zeroth-and first-order factors in eccentricity of the Hamiltonian can be expressed as The derivatives of W l ij as a function of Λ 2 take the form 1/2 (α 23 ) (54d)