Issue 
A&A
Volume 625, May 2019



Article Number  A7  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201935259  
Published online  29 April 2019 
The role of dissipative evolution for threeplanet, nearresonant extrasolar systems
^{1}
Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Université Côte d’Azur, Nice, France
email: gabriele.pichierri@oca.eu
^{2}
Division of Geological and Planetary Sciences, California Institute of Technology,
1200 E. California Blvd., Pasadena, CA 91125, USA
Received:
12
February
2019
Accepted:
13
March
2019
Early dynamical evolution of closein planetary systems is shaped by an intricate combination of planetary gravitational interactions, orbital migration, and dissipative effects. While the process of convergent orbital migration is expected to routinely yield resonant planetary systems, previous analyses have shown that the semimajor axes of initially resonant pairs of planets will gradually diverge under the influence of longterm energy damping, producing an overabundance of planetary period ratios in slight excess of exact commensurability. While this feature is clearly evident in the orbital distribution of closein extrasolar planets, the existing theoretical picture is limited to the specific case of the planetary threebody problem. In this study, we generalise the framework of dissipative divergence of resonant orbits to multiresonant chains, and apply our results to the current observational census of wellcharacterised threeplanet systems. Focusing on the 2:1 and 3:2 commensurabilities, we identify three threeplanet systems, whose current orbital architecture is consistent with an evolutionary history wherein convergent migration first locks the planets into a multiresonant configuration and subsequent dissipation repels the orbits away from exact commensurability. Nevertheless, we find that the architecture of the overall sample of multiplanetary systems is incompatible with this simple scenario, suggesting that additional physical mechanisms must play a dominant role during the early stages of planetary systems’ dynamical evolution.
Key words: planets and satellites: dynamical evolution and stability / planets and satellites: formation
© G. Pichierri et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The search for exoplanets in recent years has uncovered a multitude of planetary systems, the study of which is the key to an understanding of planetary formation and evolution. Currently, the exoplanet population is dominated by Kepler’s transit detections, making the planetary physical radii and orbital periods the better constrained parameters of the sample. Concerning the first aspect, much work has been done recently to understand how photoevaporation sculpts the physical radii of planets (Fulton et al. 2017, and references therein). In this work we address the second, complementary problem of the orbital period distribution. One of the most notable aspects of the Kepler data is that the distribution of the period ratios of neighbouring planets in multiplanets systems shows two seemingly conflicting features: on the one hand, it appears relatively broad and smooth, without any single, unmistakably emerging feature; on the other hand, a slight preference for nearresonant configurations is evident upon close examination. In fact, it is often pointed out that there is a lack of planet pairs in correspondence with period ratios very close to lowinteger ratios, and a definite excess just wide of these values, especially the 2:1 and 3:2, see Fig. 1.
Numerical simulations show that compact chains of mean motion resonances are a common outcome of slow, convergent orbital transport of planets within protoplanetary discs. Although details of discdriven migration remain an active topic of research, it is clear that such a process should play some role in the dynamical history of planetary systems. For example, it is not easy to envision a formation narrative which does not require convergent migration for systems such as Trappist1, a star famously hosting seven planets with period ratios very close to small integer ratios (Gillon et al. 2016, 2017, Luger et al. 2017). Indeed, Gillon et al. (2017) performed Nbody integrations with the orbital fits as initial conditions and these went unstable over timescales 10 000 times shorter than the estimated age of the system; in contrast, Tamayo et al. (2017) remarked that if an initial condition which results from capture into resonance through migration is chosen, then the system is stable over timescales two orders of magnitude longer then the ones found in Gillon et al. (2017). They also note that the addition of tidal eccentricity damping should help maintain the evolution stable over the system’s age. Other good examples of systems necessarily sculpted by migration are the four subNeptune planets of Kepler223 (Mills et al. 2016) and the nowclassic example of Laplacelike resonance in GJ876 (Rivera et al. 2010; Batygin et al. 2015).
In light of the fact that convergent migration should lock planets into mean motion commensurability, it is pertinent to ask how we can explain the lack of planets with exactly resonant period ratios and the excess just wide of them. Analytical models of resonance do predict that a pair of planets in a first order mean motion resonance need not satisfy the exact resonance condition (where a_{1} and a_{2} are the semimajor axes of the inner and outer planet, respectively, and k is a positive integer), but they can reside wide of resonance while the resonant angles are still librating. This divergence of the resonant equilibrium configurations happens at vanishingly low eccentricities and is linked to a fast precession of the perihelia, which is well understood analytically. However some Kepler systems are so wide of resonance that, after the resonant configuration is attained and the disc of gas is dissipated, an auxiliary mechanism might need to be invoked which actively drives these planets farther away from the exact resonance. As we see in Sect. 3, observations show that a significant fraction of nearly resonant systems lie up to 50 times wider from the resonance than the typical resonant width, and at lower eccentricities than are expected forsuch planets captured in resonance via migration in protoplanetary discs. These observations can potentially be interpreted as evidence for dissipative processes acting on the planetary systems after the disc phase.
Papaloizou & Terquem (2010) considered the specific case of the Kdwarf HD 40307, which hosts^{1} 3 hot superEarths/miniNeptunes with both pairs wide of the 2:1 mean motion resonance, and planetary masses obtained with Radial Velocity. They showed that as tidal interaction between the planets and the star reduces the eccentricities, the system maintains the libration of the resonant angles even when the period ratios are considerably far away from exact commensurability. Subsequently, Batygin & Morbidelli (2013a), Lithwick & Wu (2012) and Delisle & Laskar (2014) showed that two planets in mean motion resonance repel each other as energy is lost during tidal evolution. They thus proposed this as a viable mechanism to explain the observed distribution of period ratios in exoplanetary systems. We note that, for two planets, this repulsion can be easily understood if one considers that any process that dissipates the energy, E ∝−m_{1}∕a_{1} − m_{2}∕a_{2}, and at the same time conserves angular momentum, , should give rise to such an evolution. Indeed, this study applies to any dissipative evolution that maintains constant the angular momentum, not just tidal dissipation, and not just resonant coupling (Delisle et al. 2012).
Thanks to these studies, the case of two planet system is well understood. However the data also contains numerous systems of more than two planets (Fig. 1). Accordingly, in this paper we aim to expand the study to detected extrasolar systems of three planets. More specifically, we envision the following scenario for the formation and evolution of these planetary systems. First, three planets are embedded in the protoplanetary disc in which they formed; they interact with the disc, which ultimately results in a resonant capture. Then, after the disc is slowly depleted, the dissipative effects mentioned above are introduced, leading to orbital divergence.
Naturally this is a simplified and idealised scenario. In reality, we still do not know with enough accuracy the final configuration obtained by multisystems migrating in a disc of gas. One approach towards a better approximation would be to perform full hydrodynamic simulations of planets immersed in their protoplanetary disc accounting for various disc parameters (such as disc surface density, turbulence, opacity, etc.). This approach would however be very expensive computationally. Moreover, to date we have virtually no direct observations of the specific physical processes acting during planet formation and evolution in the early epochs of the disc phase, so these simulations, no matter how exhaustive in terms of the implementation of the plausible physics, cannot yet be directly constrained by the available data. In any case, the fact that slow convergent orbital transport strongly favours resonant captured states is well supported both analytically and numerically, as well as by the specific observations of multiplanets systems mentioned above.
In this paper we focus on slow convergent TypeI migration in a disc of gas, and adopt simple synthetic analytical formulæ for the work and the torque generated by the disc on the planets (Cresswell & Nelson 2006, 2008); the requirement that exact prescriptions for the interaction between the planets and the disc be implemented will be relaxed, invoking the aforementioned arguments favouring the plausibility of mean motion resonant capture. A similar reasoning can be applied for the postdisc phase. In order to simulate the dissipative forces that can act on the planetary system, we will implement tidal dissipation. Of course, the tidal parameters for these planets are not known (as we do not yet have a precise understanding of the interior structure of these bodies or the specific physical mechanisms that dominate the dissipation), which would pose additional questions concerning for example the timescales over which this type of dissipation takes place. However, the specific choice of tidal dissipation is only one possible example of a process such that Ė < 0 and . We conclude that our specific implementation of TypeI migration and tidal dissipation after the disc removal is therefore not restrictive, which makes our results generalisable to any other equivalent processes. In the light of these considerations, we ask if it is possible to reproduce the observed orbital configuration of real exoplanetary systems which reside close to resonance, assuming that planets are captured into resonance and undergo dissipative evolution after the disc phase. In other words we ask if the aforementioned physical processes are compatible with the distribution of nearresonant period ratios that emerges from available data.
In order to answer this question, we examined the NASA Exoplanet Archive^{2} and selected confirmed threeplanet systems for which both planet pairs lie close to a first order mean motion resonance, in particular the 2:1 and 3:2 resonances, as these seem to be the most common in the Kepler data. Our aim is to analyse these systems’ orbital parameters and to evaluate quantitatively how close they are to a multiresonant chain, which would be indicative of a dynamical history characterised by the physical processes described above.
Evidence suggesting that planets around Trappist1 and Kepler223 truly reside in a resonant configuration has recently been marshalled from the observed libration of the threebody Laplace angles. To this end, recall that if two neighbouring pairs of planets in a multiplanet systems are in the k^{(in)}:(k^{(in)} − 1) and k^{(out)}:(k^{(out)} − 1) resonances respectively (so that the resonant angles k^{(in)}λ_{2} − (k^{(in)} − 1)λ_{1} − ϖ_{2} and k^{(out)}λ_{3} − (k^{(out)} − 1)λ_{2} − ϖ_{2} are librating), then the Laplance angle φ_{L} = (k^{(in)} − 1)λ_{1} − (k^{(in)} + k^{(out)} − 1)λ_{2} + k^{(out)}λ_{3} will be automatically librating as well. The advantage of examining this threebody angle over the twobody resonant angles is that the latter contain the longitudes of the pericenters ϖ, whose precession rates are poorly constrained by the data, while the former only includes the mean longitudes λ whose derivatives in time are directly deduced by the transit observations. However, solutions for which the resonant angles were originally in libration around a resonant equilibrium point can become circulating when the eccentricity of the equilibrium point becomes small enough under the effect of tidal damping (Delisle et al. 2015), and, similarly, even a small distance from the equilibrium point could be responsible for breaking the libration of the threebody Laplace angle when the equilibrium eccentricity becomes small enough. Therefore, even if such circulations of the angles were observed, this would not be in disagreement with the envisioned scenario of resonant capture and subsequent dissipative evolution. In other words, the libration of the Laplace angle is a sufficient, but not necessary condition for past resonant capture in a chain of firstorder resonances.
We therefore perform here a different analysis of the observed data, where we do not attempt to verify that a given system resides formally in resonance at the present day, but instead we evaluate the distance of a system from the considered resonance chain and the probability that this proximity is due to mere chance. In order to do this, we look for resonant solutions that provide the closest match to the observed planetary orbital configurations, that is the semimajor axis ratios. It is worth anticipating here the following important point. As it will be clear later (see Sect. 3.2), in the case of only two resonant planets residing wide of resonance it is always possible to find a resonant configuration which matches the observed data. This is because the eccentricities of these planets are at the present day not well constrained observationally, making the total orbital momentum of the system a free parameter: it is therefore always possible to find a value of that reproduces the observed a_{2}∕a_{1} with resonancelocked orbits. However, this is not the case for systems of three planets, since we still have only one free parameter (whose value is linked to the initial captured state, not constrained observationally) but two observables, that is the two pairs’ semimajor axis ratios.
As detailed below, we carried out our study of finding orbital configurations that match the observed data using both an analytical and a numerical approach. The semimajor axes of the planets may be inferred from the orbital periods once the stellar mass is known, however this quantity is not yet well constrained in all cases. Nonetheless, all we will be interested in will be the semimajor axes ratios a_{2}∕a_{1} and a_{3}∕a_{2}, which can be obtained without any knowledge of the mass of the star directly from the period ratios and using Kepler’s third law. This is tantamount to renormalising all separations by some arbitrary length, which does not affect the underlying physics since the dynamics only depends on the ratios of the semimajor axes and not on their individual values (only the timescale of the evolution does).
For the purposes of this study, we can limit ourselves to an analysis to first order in the planetary eccentricities. Indeed, the eccentricities that are expected for planets that have been captured into mean motion resonance by slow convergent migration in a disc are of order , where τ_{a} and τ_{e} are the timescales of migration and eccentricity damping respectively, and h = H∕r ~ 0.05 is the aspect ratio of the disc (Goldreich & Schlichting 2014, Pichierri et al. 2018). Since discs with high aspect ratios are not expected, the limit of small eccentricity is justified, and even more so in the phase of dissipative tidal evolution, which further damps the eccentricities. Moreover, given that these are transiting planets, and that during the disc phase any mutual inclination of the planets would be damped out, we assume coplanar orbits for simplicity. Another useful piece of information which is available to us is the radii of the planets. This could in principle be used to infer the planetary masses (e.g. Wu & Lithwick 2013). However the radius–mass relationship in Kepler planets is marked by extreme scatter (Weiss et al. 2013), and we therefore choose to keep the planetary masses as a free parameter. More specifically, we are only interested in the mass ratios m_{1} ∕m_{2} and m_{2} ∕m_{3}, since, as we will show, they are the only dynamically significant quantities that can affect the values of the semimajor axis ratios (see also Appendix A).
The remainder of this paper is organised as follows. In Sect. 2 we obtain an analytical model for three planets in a chain of first order mean motion resonances, valid in the limit of small eccentricities. With this analytical model, we find the stable resonant configurations and map them in terms of the orbital elements. Finally we obtain an analytical confirmation of resonant repulsion for threeplanets systems undergoing dissipation. In Sect. 3 we detail our study, employing both analytical and numerical methods. We select systems of three planets near mean motion resonances, focusing on the 2:1 and 3:2 resonances, and we analyse their orbital configuration using the available data in order to evaluate if they are consistent with the process of resonant capture and subsequent dissipative evolution. We present our results in Sect. 4 and we finally conclude by discussing their significance in Sect. 5.
Fig. 1 Observations of planethosting stars reveal that multiplanetary systems are not rare, hosting over 1600 confirmed planets (panel a). The period ratio distribution of neighbouring planets is shown in panel b. One can observe an overall broad distribution as well as a number of peaks slightly wide of resonant ratios, especially the 2:1 and 3:2 commensurabilities. Data was obtained from the Nasa Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu/). Panel a: histogram of multiplanetary system by number of planets in each system. Panel b: distribution of period ratios of neighbouring planets in exoplanetary systems. 
2 Planetary Hamiltonian
The Hamiltonian of two resonant planets in the limit of low eccentricities has been studied extensively in the literature (e.g. Batygin & Morbidelli 2013b, and references therein). Collectively these studies have pointed out that even if both planets are massive and to first order in eccentricity it is possible to reduce the problem to an Hamiltonian that is analogous to the wellknown Hamiltonian of the restricted, circular threebody problem of a massless particle in resonance with a massive unperturbed body. In particular, such a Hamiltonian is integrable and is equivalent to the socalled second fundamental model for resonance (Henrard & Lamaitre 1983). This is, however, not the case for three planets. Nonetheless, it is useful to extend an analytical description of the resonant dynamics at low amplitude of libration of the resonant angles in the case of three planets orbiting a star. In this section, we introduce the Hamiltonian of the system, derive curves representing the loci of its stable equilibrium points, and show how these can provide a description of a system along the dissipative evolution. We will apply this model to real Kepler system in Sect. 3.2.
Consider three planets of masses m_{1}, m_{2} and m_{3}, orbiting around a star of mass M_{*} in a canonical heliocentric reference frame (Poincare 1892). Indices 1, 2 and 3 will refer to the inner, middle and outer planet, respectively.As usual, we consider the planetary Hamiltonian, which we write as (1)
where the keplerian part is given by (2)
and describes the (integrable) motion of the three planets due to their interaction with the star, to which the small perturbation is added, which includes all the mutual interactions between the planets. We now assume that the inner pair of planets is close to a k^{in}:(k^{in} − 1) mean motion resonance, and that the outer pair of planets is close to a k^{out}:(k^{out} − 1) mean motion resonance, where k^{in}, k^{out} > 1 are two positive integers. In other words, we assume the resonance conditions n_{1} ∕n_{2} ≃ k^{in}∕(k^{in} − 1), n_{2} ∕n_{3} ≃ k^{out}∕(k^{out} − 1), where for each planet is the mean motion. Since we are interested in the resonant interaction between the planets only, we will average the Hamiltonian over the fast evolving angles so that only combinations of the resonant angles k^{in} λ_{2} − (k^{in} − 1)λ_{1} − ϖ_{1}, k^{in} λ_{2} − (k^{in} − 1)λ_{1} − ϖ_{2}, k^{out} λ_{3} − (k^{out} − 1)λ_{2} − ϖ_{2}, and k^{out} λ_{3} − (k^{out} − 1)λ_{2} − ϖ_{3} remain in the Hamiltonian, where λ is the mean longitude of a planet, and ϖ is its longitude of the periastron.
The resonant perturbing Hamiltonian expanded to first order in the eccentricities reads (3)
where the orbital elements are constructed from heliocentric positions and barycentric velocities (Poincare 1892). The coefficients f_{res} are typically of order unity, and it is straightforward to determine the strength of each resonant harmonic, and incorporate direct and indirect terms. They depend (weakly) on the semimajor axis ratios, and their expressions may be found in Murray & Dermott (2000). We therefore write the Hamiltonian after the averaging procedure as (4)
and then drop the higher order terms. We note that terms that describe the mutual influence of the innermost and outermost planet are not included in as this is a higher order effect. Note also that by dropping the higher order terms the problem is reduced to a planar one. In order to maintain the canonical nature of the equations of motion, we introduce for each planet the modified Delaunay actionangle variables (Λ, Γ, λ, γ) (omitting the subscripts 1,2,3 for simplicity), which are given in terms of the orbital elements by (5)
where μ = M_{*}m∕(M_{*} + m) is the reduced mass, and ℓ = E − e sin E is the mean anomaly (E being the eccentric anomaly). In these variables, the Keplerian part of the Hamiltonian (1) takes the form (6)
while the resonant Hamiltonian writes (7)
we note that in going from Eqs. (3) to (7) we have made use of the approximation , which holds at first order in e.
This Hamiltonian is clearly not integrable. However, one can perform a series of changes of variables that allow us to reduce by two the number of degrees of freedom. The first canonical transformation is (8)
it is straightforward to check using the Poisson bracket criterion (Morbidelli 2002) that it is indeed canonical. Now, the new angle κ does not appear explicitly in the Hamiltonian, which makes a constant of motion. The significance of has already been discussed for two planets (e.g. Michtchenko et al. 2008, Batygin & Morbidelli 2013b), and it has to do with the location of exact resonance. As we have already mentioned, neighbouring planets can still be in resonance while their semimajor axis ratios do not satisfy exactly the resonant condition , therefore the observed a_{i,obs} do not alone reveal how close the planets are to resonance, nor do they represent the nominal ā_{i} that do satisfy it. However by calculating from a_{i,obs} the value of the constant of motion , and imposing in the formula (9)
the condition of exact resonance, , for all pairs i = 1, 2, we derive the nominal value of . From this, one easily obtains ā_{3} from , and then recursively , and finally .
It is worth briefly recalling here why even in resonance the planets’ semimajor axes do not coincide exactly with their nominal values. As an example, for the inner planet, the condition for resonance is that the resonant angle k^{in} λ_{2} − (k^{in} − 1)λ_{1} − ϖ_{1} is librating, , which together with the condition of exact nominal resonance k^{in}n_{2} − (k^{in} − 1)n_{1} = 0 would imply ; however from the Hamiltonian (Eq. (7)) we have which grows as e_{1} ↘ 0, meaning that at low eccentricities , which in turn forces in order to maintain the libration of the resonant angle. The resonant equilibrium points will therefore correspond to semimajor axes a_{i} which may well deviate farther and farther from ā_{i} as e_{i} approaches 0. We can however already greatly simplify the calculations given that we will only consider deviations of the semimajor axis ratios from the nominal ones of no more than 5% (moreover, very small values of the eccentricities are observationally disfavoured for Kepler systems, Hadden & Lithwick 2017). In this limit, we can expand the Keplerian part to second order in , where is the nominal resonant value of Λ_{i}, and write (10)
which, inserting the definition of δΛ_{i} and dropping the unimportant constant term and the higher order terms, reduces to: (11)
where is the nominal mean motion and can be interpreted as the inverse of the moment of inertia of a circular orbit around the star. As we will see below, for the purposes of our study, considering the expanded Keplerian Hamiltonian up to order does not introduce any significant inaccuracy in our calculations. Concerning the resonant Hamiltonian (Eq. (7)), we can evaluate it at the nominal values as it is already of order .
Finally, one last canonical change of variable is made: (12)
Again, we see that the new angle θ′ does not appear in the Hamiltonian, making Ω another constant of motion of the system (we note that here Ω does not denote the longitude of the node which does not appear in our model, since the problem is planar). We are therefore left with a fourdegreeoffreedom Hamiltonian which depends parametrically on the constants of motion , Ω. We already mentioned the meaningof ; for Ω, one can easily show that , the total angular momentum of the system, which is to be expected knowing that it is a conserved quantity.
2.1 Resonant equilibrium points
Let us briefly summarise our work so far. We have obtained a 4degreesoffreedom Hamiltonian which is a function of the actions and the angles , and depends parametrically on the values of and Ω (which are linked to the orbital elements as expressed in Eqs. (5) and (12)); the Hamiltonian in these variables reads (13)
where the explicit dependence of each term can be obtained by direct substitution. We now consider the stable equilibria of this system. We look for equilibrium points of this Hamiltonian by simultaneously solving the set of equations (14)
We note that by the functional form of the Hamiltonian, any combination of values in {0, π} for the angles immediately satisfies the last line. These are known as the symmetric equilibria. Asymmetric equilibria are possible (e.g. Beaugé et al. 2006), but they do not play a role at the low eccentricities at which we are limiting ourselves here.
Plugging in specific values for the angles in {0, π} reduces the problem of solving the four equations that appear in the first line to find the stable equilibria of the system. We note that although the Hamiltonian depends on both Ω and , the latter assumes a natural value for any specific problem at hand (that is, any values of m_{1}, m_{2}, m_{3} and of k^{in}, k^{out}) by rescaling the units so that for example ā_{1} = 1. To trace out the loci of the resonant equilibria, we then simply change the value of Ω (which corresponds to changing the angular momentum , at constant ) and solve Eq. (14) to find , which are then translated into orbital elements working backwards through the canonical transformations.
We show in Fig. 2 one example of equilibrium curves for three equalmass planets down to eccentricities of order 10^{−3}, where we also show that the expanded Keplerian Hamiltonian provides an accurate description of the system. These curves are then matched against the result of full Nbody numerical simulations of a system with the same physical parameters which starts deep in resonance and evolves dissipatively so to slowly follow the resonant equilibrium points (transparent lines)^{3}. These Nbody integrations with the addition of dissipative effects will be detailed below in Sect. 3.3.
Fig. 2 Equilibrium curves showing the loci of the stable resonant equilibrium points calculated as explained in the text, in the case of a 3:2 – 3:2 mean motion resonance chain, with m_{1} = m_{2} = m_{3} = 10^{−4}M_{*}. The full curves are calculated using the expanded Keplerian Hamiltonian (Eq. (11)), while the dashed curves are calculated using the unexpanded Keplerian Hamiltonian (Eq. (6)), showing very little difference down to very small eccentricities and for reasonable values of the nearly exactly resonant semimajor axis ratio. The location of the nominal resonant semimajor axis ratio (3/2)^{2/3} is shown bya vertical orange line. We also superimpose the numerically computed evolution of a threeplanet system deep in the 3:2 mean motion resonance (for both pairs) and undergoing dissipative evolution depicted with transparent lines: the system follows the locations of the equilibrium points, which are close to the curves calculated analytically. 
2.2 Resonant repulsion for threeplanets systems
The equilibrium curves in the a_{i+1}∕a_{i} vs. e_{i} plane show that the resonant repulsion during energy dissipation is expected also for threeplanets systems. For systems of two planets, it is well known that for first order resonances the resonant equilibria always reside wide of the exact resonant ratio of the semimajor axes. That is, because the resonant condition requires and since the perihelion precession is always retrograde, , one necessarily has , that is . More concretely, at low enough eccentricities and at semimajor axis ratios close to the nominal ones, one finds directly using the resonant Hamiltonian expanded tofirst order in e that , with , : this means that the lower are the eccentricities, the wider are the equilibria from the exact commensurability. At higher eccentricities the secular terms, of , become more important, and they contribute a positive contribution (that is constant in e) in ; however, one still finds at higher eccentricities as well (e.g. Pichierri et al. 2018).
For systems with three planets, since we used a first order expansion in e in the analytical model and therefore we are not including the mutual interaction between the inner planet and the outer planet, the perihelion precession will still be retrograde and it will remain true that for each pair of planets the resonant equilibria lie wide of the exact nominal resonance, and that, in the limit of small enough eccentricities, the separations grow with diminishing eccentricities. This is indeed what we see in Fig. 2, where the analytically computed resonant equilibria agree very well with our numerical simulations^{4}.
Consider now a resonant threeplanet system that is close to some resonant equilibrium point and is subjected to (tidal) dissipation. Assuming that the dissipative evolution is slow compared to that of the resonant variables, which has a characteristic timescale given by the libration period at vanishing amplitude of libration, the system will remain bound to the equilibrium curves. Since the eccentricities are damped by the dissipation, we conclude that the semimajor axes are expected to diverge. We also conclude that systems of three planets that are close to a given firstorder mean motion resonance but for which one or both pairs is narrow of the resonance can only be explained by a resonant configuration if the amplitude of libration around the resonant equilibrium point is large, and temporarily takes the planets to period ratios that are lower than the exact resonant period ratios.
We finally note here a property of these curves that will be used later. The Hamiltonian (Eq. (13)) can be rescaled by a parameter which encapsulates all of the information regarding how the dynamics scales with mass ratios and physical sizes of the orbits. This is analogous to the rescaling found in Batygin & Morbidelli (2013b) for the 2planets case, and only works when using an expanded Hamiltonian and for semimajor axes close to the nominal resonant ones, which are our working assumptions anyway. Therefore, after rescaling all planetary masses by a certain factor , the corresponding loci of the resonant equilibria are also simply rescaled, and can be immediately calculated. More specifically, one can easily see that for given semimajor axis ratios, the values of the eccentricities that correspond to the resonant equilibrium point are just rescaled by , since the eccentricities and the planetary masses appear as a product in the perturbing Hamiltonian. This can be easily understood using the previous formula , and noticing that fixing the semimajor axis ratio ultimately fixes by the resonance condition; therefore, rescaling the planetary masses, at fixed (i.e. at constant semimajor axis ratios), the eccentricities are simply rescaled by the same factor. This implies that for a given equilibrium configuration of the semimajor axis ratio a_{2,eq}∕a_{1,eq}, the corresponding equilibrium of the ratio a_{3,eq}∕a_{2,eq} will be independent of , that is independent of the absolute value of the planetary masses. Only the ratios m_{1} ∕m_{2} and m_{2} ∕m_{3} are significant, meaning that if one changes one of these ratios, the equilibrium in a_{3,eq}∕a_{2,eq} corresponding to the same a_{2,eq}∕a_{1,eq} will have changed (see Appendix A for an explicit presentation of this rescaling procedure).
3 A scenario for dissipative evolution of threeplanet systems
In this section, we select nearresonant systems of 3 planets from the NASA Exoplanet Archive catalogue, and discuss whether or not their observed orbital configuration is compatible with the dynamical evolution driven by the following physical processes. As we already mentioned, planets are expected to dynamically interact with the protoplanetary disc in which they formed. This has two effects: a damping of the eccentricities (and inclinations), and an exchange in angular momentum between the planets and the disc which causes the planets’ semimajor axes to change (planetary migration). Planets captured in mean motion resonance are usually attributed to inward migration (that is, the planet looses orbital angular momentum to the disc), and we will consider this case in this paper, but these results are general to any form of convergent evolution. In our implementation, when the inner planet reaches the inner edge of the disc, migration is stopped, and the second incoming planet will cross a commensurability with the first; finally, the third planet will cross a commensurability with the second. We note that the timescale over which the eccentricity is damped is usually of order ~ 100 times shorter than that over which the semi major axis changes. Planets are therefore expected to approach these commensurabilities with vanishingly low eccentricity. As we can see from Fig. 2, however, this configuration with semimajor axis ratio wide of resonance and low eccentricity is very close to the resonant equilibrium points. The planets keep approachingeach other and their eccentricities keep increasing due to the curvature of the locus of the equilibrium points in the (a_{i+1}∕a_{i}, e_{i}) diagram, until an equilibrium is reached when the damping of e balances such resonant eccentricity pumping. This is why planets are expected to be (close to) a resonant equilibrium point in the first place, and a chain of resonances can be formed. The disc is then slowly depleted and the planets maintain their configuration.
Following the depletion of the gas, dissipation is introduced which removes orbital energy at constant angular momentum; this is done here implementing tidal dissipation but again the method is general. During this phase of dissipative evolution, the planets will follow again the equilibrium curves of the resonant Hamiltonian for changing Ω, this time decreasing their eccentricities and hence increasing the semi major axis ratios a_{i+1} ∕a_{i} for each planet–planet pair. We note that Ω changes because changes (since do the semimajor axes as a consequence of the dissipation of energy) and has to stay constant for this kind of dissipation.
3.1 Choice of systems
Recall that the only orbital parameters that are well constrained by transit data are the orbital periods, which allow us to obtain the semimajor axis ratios even without knowing the mass of the star. The orbital periods listed in the NASA Exoplanet Archive catalogue are, for the cases considered below, obtained by fase folding the observed signal. Since we will be considering shortperiod planets this is equivalent to obtaining the proper value of the periods, so that we can directly compare the corresponding observed semimajor axis ratios with the ones coming from our averaged model of resonance^{5}. The masses of the planets could be obtained starting from the estimates for the planetary radii, and making use of a scaling relation such as the one found in Wu & Lithwick (2013), m_{pl} = 3m_{⊕}(R_{pl}∕R_{⊕}), where m_{pl}, R_{pl} are the mass and radius of the exoplanet, and R_{⊕}, m_{⊕} those of the Earth. However this is only a statistical law and the uncertainties on the mean densities usually preclude accurate estimates for the masses, which are indeed not yet known. We will therefore use the masses as free parameters of our study. In fact, as we already saw, the only significant quantities for our study are the mass ratios between the planet; this follows from the discussion at the end of Sect. 2.2. In practice, we will choose a total planetary mass m_{tot} = m_{1} + m_{2} + m_{3} by using for each system the mean planetary radius (R_{1} + R_{2} + R_{3})∕3, the aforementioned scaling relationship from Wu & Lithwick (2013) to obtain an average planetary mass m_{pl,avg}, and setting m_{tot} = 3m_{pl,avg}. Again, this is simply a choice that we are forced to make in order to run Nbody simulations, but it does not in any way change our result, which is therefore not sensitive to the uncertainties on the radii (or to the lack of their knowledge, as will be the case for YZ Ceti). Note however that since individual Kepler systems seem to show a homogeneity in planetary radii and masses (Weiss et al. 2018; Millholland et al. 2017), this choice likely constitutes a good approximation to the real architecture of these systems.
Given a pair of neighbouring planets with periods P_{1} and P_{2} respectively, which are close to a given first order resonance, P_{1}∕P_{2} ≃ (k − 1)∕k, one can define (see for example Lithwick et al. 2012; Hadden & Lithwick 2014) (15)
called the normalised distance from (exact) resonance. When Δ > 0 the planets reside wide of the k : k − 1 resonance, while when Δ < 0 the planets are narrow of the resonance.
We will be selecting planetary systems of three planets with both pairs close to some first order mean motion resonance such that Δ ≲ 0.05 holds for both pairs, with k = 2, 3 as they appear to be the most common resonances. We recall that the normalised width of a resonance is of order (Deck et al. 2013; Batygin 2015), where the average planetary mass for Kepler systems is of order m ~ 3 × 10^{−5}M_{*}, and moreover most planets in each system appear to be quite homogeneous in mass (Weiss et al. 2018; Millholland et al. 2017). This gives a typical resonance width of order Δ ~ 10^{−3} in normalised units, meaning that in selecting systems with 0 < Δ ≲ 0.05 we are generously including planetary pairs with separation 50 times larger than the typical resonant width. Moreover, the available data shows that for systems close to mean motion resonance, the distribution of Δ favours values between 0 and ≲0.05 (Hadden & Lithwick (2017), see also Fig. 3). Additionally, we will require that Δ > 0, which is justified by our results in the previous section.
Of all threeplanet systems, only 8 satisfy Δ < 0.05 for both pairs, that is, appear to be close to a multiresonant chain (they are Kepler31, Kepler53, Kepler114, Kepler207, Kepler289, Kepler305, Kepler326, YZ Ceti). The architecture of these systems is shown in Fig. 4; of these, only 3 satisfy Δ > 0 for both pairs. These are Kepler31, Kepler305 and YZ Ceti. For these systems, we consider whether or not their observed orbital configuration can be consistent with the scenario envisioned above.
Fig. 3 Distribution of the (signed) normalised distance from first order mean motion resonances with k =2, 3 in all exoplanetary systems (yellow) and for threeplanet systems (purple). We note a clear peak to the right of the value Δ = 0 (corresponding to the exact nominal resonance, indicated by a vertical orange line), which is most prominent for 0 < Δ ≲0.05. Out of the 358 pairs plotted here, there are 123 total pairs in this configuration. For the threeplanets systems only, 48 pairs have 0 <Δ ≲ 0.05, out of the 121 shown in the purple histogram. 
Fig. 4 Orrery of the threeplanets systems sufficiently close to first order mean motion resonances k : k − 1 with k =2, 3, with a normalised distance to resonance Δ < 0.05 for both pairs. For each system, we place a circle in correspondence of the period of each planet, and indicate between pairs of planets the first order mean motion resonance in which they are envisioned to reside (below) and the normalised distance to that resonance Δ (above); the sign of Δ indicates if the pair of planets are narrow (Δ < 0) or wide (Δ > 0) of the resonance. For our analysis, we will only consider systems for which Δ > 0 for both pairs (the systems enclosed by a box). The size of the circle is an indication of the estimated radius of the planet (the small dot in the top right corner demonstrates the size of Earth). For YZ Ceti this information is not available, but this does not pose a problem for our study. 
3.2 Analytical maps
With our analytical model of resonance at hand, we can construct analytical maps of resonant equilibrium points for different resonant chains and different planetary mass ratios. For the purposes of the current study, we proceed as follows. For an arbitrary system, we assume to have observations for the orbital period ratios and obtain the values of k^{(in)} and k^{(out)}. We then pick both mass ratios m_{2}∕m_{1} and m_{3}∕m_{2} and construct the equilibrium curves as explained in Sect. 2.1 (in practice, we work with the aforementioned Hamiltonian rescaled by the common planetary mass factor , see Appendix A). Then, we find the resonant equilibrium point (i.e. the value of ) that corresponds to a value of a_{3} ∕a_{1} equal to the observed semimajor axis ratio, and therefore put . The determined value of fixes the eccentricities of all three planets, since they are all linked by the equilibrium curves.
Recall that we only have one free parameter to select the chosen equilibrium configuration: the angular momentum; however, we have two observables that we want to match, which are both semimajor axis ratios a_{2} ∕a_{1} and a_{3} ∕a_{2}. This is unlike the case of only two resonant planets, where one has one free parameter (again the angular momentum ) and only one observable (the single a_{2}∕a_{1} ratio): in this case, it would always be possible to find a suiting value of which gives a resonant equilibrium configuration such that the semimajor axis ratio is equal to the observed one (provided that the latter is wide of the nominal value ). In the case of three planets, choosing such that automatically fixes the equilibrium values of both semimajor axis ratios a_{2} ∕a_{1} and a_{3} ∕a_{2}. Considering for example the corresponding equilibrium , we obtain the weighted difference (16)
between and the observed value . The same can be done for a_{2}∕a_{1}, which gives a similar (absolute) result, . Maintaining this procedure, we loop over different planetary mass ratios.
We applied this procedure to the three systems selected above, starting with Kepler305, which resides close to a 3:2–2:1 mean motion resonance chain. First of all, to better represent what these analytical maps intend to show, we draw in Fig. 5 equilibrium curves (equivalent to those shown in Fig. 2) which describe the locations of the resonant equilibria for this resonant chain and for one specific choice of mass ratios m_{1} ∕m_{2} = m_{2}∕m_{3} = 1. The observed values of the semimajor axis ratios are indicated by dashed vertical lines; then, we indicate with a red dot the location of the specific equilibrium point that is selected when we impose ; finally, we obtain as defined in Eq. (16) (and similarly for ). As explained above, imposing automatically fixes all equilibrium eccentricities e_{1,eq}, e_{2,eq}, e_{3,eq}, and we can store their maximum max{e_{1,eq}, e_{2,eq}, e_{3,eq}} in to better describe the orbital configuration at the selected equilibrium point. We choose to consider the quantity maxe_{eq} rescaled by a common planetary mass factor in order to obtain results that are independent of the planettostar mass ratio, since again the latter quantity does not effect equilibrium semimajor axis ratio configurations, and therefore does not affect . Panels a and b in Fig. 6 are the result of this procedure spanning different planetary mass ratios, showing maps of and using the observed semimajor axis ratios of Kepler305 (the bottom panels c and d show the results of numerical simulations which will be detailed in Sect. 3.3, and are only intended to validate the analytical results). We show analogous results for the system YZ Cet, residing close to a 3:2 – 3:2 chain, in Fig. 7, and for For Kepler31, a chain (close to the) 2:1 – 2:1 mean motion resonances, in Fig. 8.
3.3 Numerical simulations
In order to check the validity of our analytical calculations, we turned to numerical simulations by performing the following study. We simulated the process of capture into a chain of firstorder mean motion resonances by placing the planets relatively wide of the desired resonances, according to the specific values of k^{in} and k^{out} of each case, and simulating the effects of the protoplanetary discs by adding fictitious forces which mimic the interaction with the disc (Cresswell & Nelson 2006, 2008) to the Nbody integrator swift_symba. To ensure convergent migration for all planetary mass ratios, we stopped the migration of the inner planet by adding at the desired location a socalled planetary trap, which reproduces the effect of the inner edge of the disc and describes a disc cavity around a star (Masset et al. 2006, Pichierri et al. 2018). As we mentioned above, the mass ratios were kept as free parameters. Since again we are interested in mass ratios of order unity, in our simulations we limited ourselves to m_{1} ∕m_{2} and m_{2} ∕m_{3} between ~ 0.2 and ~ 5, and repeated the same set of simulations.
Recall that the timescale for planetary migration depends on the mass of the planet, meaning that changing the mass ratios will change the relative speeds at which each planet’s semimajor axis decreases, which is practically inconvenient. Therefore, we used a fictitious τ_{a} which is kept equal for all planets and constant along the different simulations of resonant capture. This has the sole effect of making it easier to automate the simulations, and does not affect our results. We need only to make sure that at the end of the discmigration phase the semimajor axis ratios are smaller than the observed ones, since the subsequent evolution dominated by tidal dissipation will only cause the semimajor axis ratios to expand. We note that this is always possible, since one can obtain different final eccentricities at the captured state by changing the ratio of the eccentricity damping τ_{e} and the migration rate τ_{a}, and thus obtain different corresponding equilibrium values of the semimajor axis ratios; the latter approach the exact commensurabilities as the eccentricities grow (Fig. 2), and since for the systems that we are studying both pairs reside wide enough of the nominal resonance, the final eccentricities need not be too high, of order a few 10^{−2} for a typical planetary mass of order 10^{−5}M_{*}.
After the desired resonant state is obtained, we slowly depleted the gas. Finally, we added the effects of tidal dissipation (following Mardling & Lin 2002), using arbitrary quality factors for the planets but large enough to ensure that the dissipative evolution be slow compared to the resonant evolution of the two planets’ pairs. This allows us to perform efficient integrations without breaking the adiabatic approximation which keeps the system close to the resonant equilibrium points. We note that we have little to no information on the internal structures of exoplanets, so we would not be able to confidently assign realistic eccentricity damping timescales anyway. Moreover, as we have already mentioned, tides are only one example of dissipation (that is, loss of orbital energy E) at constant angular momentum , so that these results are in fact generalisable to any dissipative process such that Ė < 0 and . Therefore, a resonant system undergoing any such process is expected to follow the loci of the resonant equilibria, and the divergence of the semimajor axes is obtained as a general result. We now explain how we obtain maps similar to those drawn in the previous section from these numerical simulations.
Consider a choice of the mass ratios, and a simulation of the dissipative evolution of two pairs of resonant planets. The semimajor axis ratios a_{3}∕a_{1} and a_{3}∕_{2} will increase in time. When, for two consecutive outputs of the simulation, the ratio a_{3} ∕a_{1} crosses the observed one , we store the corresponding value of a_{3}∕a_{2} from the simulation (an average of it at the two consecutive outputs). Since the system might be librating around the equilibrium points with some amplitude, and there are additional short period terms, this will happen many times for a single simulation, and we obtain a list of a_{3}∕a_{2} values. Then, we report the average of this list, and again compare this quantity with the observed (since they are obtained from the mean period extracted from the data) by computing the relative difference as in Eq. (16). We then loop over different choices of planetary mass ratios and obtain a map that can be compared with the analytical maps of the previous section. A similar procedure can be applied to a_{2} ∕a_{1} (which gives again similar values to that of a_{3}∕a_{2} as we mentioned in Sect. 3.2), as well as the quantity .
This analysis has been performed for the three selected systems. For Kepler305 and YZ Cet, we show the resulting plots on thebottom panels c and d of Figs. 6 and 7 respectively, and notice very good agreement with the analytical results. The noise that is observed in the panels c relative to the quantity is due to the fact that the numerically simulated systems are undergoing fast oscillation while they cross the observed value , but the typical value of is similar to the one found analytically.
The case of Kepler31, which resides close to a 2:1–2:1 mean motion resonant chain, is a bit different, since the 2:1 resonance capture might be only temporary if the librations around equilibrium are overstable (Goldreich & Schlichting 2014). This behaviour has been already investigated thoroughly in the case of two planets (Delisle et al. 2015, Deck & Batygin 2015), however it has been shown to be dependent on the specific implementation of the discplanet forces, and to disappear in some cases (Xu et al. 2018). In this work we do not intend to expand on these matters, since the formulas that mimic the planetdisc interactions represent only approximate implementations of the real forces that are felt by the planets from the disc, which themselves remain observationally unconstrained. We therefore take a practical approach, and note that in the numerical simulations where resonant capture was successful (typically for m_{1} ∕m_{2}, m_{2}∕m_{3} ≳ 1) the numerical results agree very well with the analytical ones; moreover we still observe that the theoretical value of varies extremely little with m_{1}∕m_{2} and m_{2}∕m_{3} (Fig. 8), so the latter simulations can be considered as enough support for the analytical calculations.
Fig. 5 Locations of the resonant equilibrium points in the a_{i+1}∕a_{i} vs. planes, i = 1, 2, for three equalmass planets in a 3:2 – 2:1 mean motion resonance chain, close to which Kepler305 resides. Orange vertical lines show the exact nominal commensurability, while dashed vertical lines show the observed a_{2} ∕a_{1} and a_{3} ∕a_{2} in the case of Kepler305. As explained in the text, we select one equilibrium configuration (indicated by the red dot in both panels) by requiring that (a_{3} ∕a_{1})_{eq} ≡ (a_{3}∕a_{1})_{obs}, which automatically fixes all orbital elements a_{2}∕a_{1}, a_{3} ∕a_{2}, e_{1}, e_{2} and e_{3} along the equilibrium curves. 
Fig. 6 Top row: analytical maps constructed for Kepler305 as explained in Sect. 3.2. Panel a: we plot the quantity , which represents how close the system is now to some resonant equilibrium point, for different mass ratios m_{1} ∕m_{2} and m_{2} ∕m_{3} (each point in this plot is constructed by repeating the procedure described in Fig. 5). We notice that changes verylittle with the mass ratios, and is of the order of ~0.002. Comparing with Fig. 9, we see that this can be the case by pure chance only in ~ 15% of randomly generated systems close to the 3:2–2:1 mean motion resonance chain. Panel b: we show a map of the quantity representing the equilibrium orbital configuration that is selected at each fixed value of m_{1} ∕m_{2} and m_{2} ∕m_{3} by imposing (a_{3}∕a_{1})_{eq} ≡ (a_{3}∕a_{1})_{obs}. Bottom row: numerical maps constructed for Kepler305 as explained in Sect. 3.3. We show numerical maps of in panel c and in panel d, analogous to the analytical plots above (over a slightly smaller range of mass ratios for simplicity). These are intended to validate the analytical maps, and show very good agreement between corresponding panels. 
Fig. 7 Same as in Fig. 6, but for the system YZ Cet, residing close to a 3:2–3:2 chain. The value of across all planetary mass ratios can be matched against the corresponding curve in Fig. 9, where we find that ~ 80% of randomlygenerated systems lie this close to the 3:2–3:2 chain. 
Fig. 8 Same as in Fig. 6, but for the system Kepler31. Only analytical maps are shown in this case since in some simulations capture into resonance was unsuccessful due to overstability of the captured state for different planetary mass ratios. As explained in the text, this issue is modeldependent and is not within the scope of our analysis. Moreover, inthe simulations where capture was successful, the results agree well with the analytical calculations, showing . Comparing with Fig. 9 we see that there is a ~20% probability that Kepler31 lies this close to the 2:1 – 2:1 chain by pure chance. 
4 Results
4.1 Probabilistic measure of a resonant configuration in Kepler305, YZ Cet and Kepler31
Using the maps of shown in Figs. 6–8 for the three selected systems Kepler305, YZ Cet and Kepler31, we draw the following conclusions. First of all, one might expect that the quantity should change with the different choices of mass ratios, thus allowing one to make a prediction on their (so far unknown) values of m_{1} ∕m_{2} and m_{2} ∕m_{3} under the assumption that these systems are indeed in resonance and evolving dissipatively. Followup monitoring of these systems could then produce new observations from which to obtain the real masses of the planets, and so validate or disprove the hypothesis. However, in practice we find that these analytical maps show very little dependence on m_{1} ∕m_{2} and m_{2} ∕m_{3} spanning reasonable values. Note also that for all three systems is small, but never vanishing, which would represent an analytically computed equilibrium configuration such that , that is, a resonant equilibrium point which satisfies and . But even if this happened to be the case, the level curve would still span a broad range of mass ratios: given moreover the uncertainty in the observed period ratios of exoplanetary systems, this would make any determination of m_{2}∕m_{1} or m_{3}∕m_{2} using observed data, in general, inconclusive.
Secondly, we note that we do obtain in all three cases small values for , meaning that these systems are indeed close to some equilibrium point of the Hamiltonian (Eq. (13)) and therefore could potentially reside in a multiresonant chain. However, these values by themselves do not contain any meaningful information. The quantity should indeed be calibrated if we intend to use it as a measure of the probability that the actual system (with its real unknown eccentricities) is in resonance, which in turn would yield a measure of how consistent the orbital architecture of such a system is with the envisioned scenario described above. To this end, for various resonant chains we randomly generate period ratios of fictional systems such that 0 < Δ < 0.05 for each pair, and extract the corresponding (calculated for the choice of mass ratios m_{1}∕m_{2} = m_{2}∕m_{3} = 1 for simplicity, since as we saw above the value depends extremely weakly on the mass ratios). From the cumulative distribution of that arises from this procedure we can obtain the probability that any given system has a given (small) purely by chance.
Since we are interested mainly in the 2:1 and 3:2 mean motion resonances, we show in Fig. 9 these cumulative distributions for systems close to any possible combinations of these resonances. The results showthat the proximity of YZ Cet to the 3:2–3:2 resonance is not statistically significant, since in ~ 80% of randomly generate systems close to the 3:2–3:2 chain we obtain an equivalent or smaller value of . Instead, Kepler305 and Kepler31 are likely to be in resonance at the 1σ level (i.e. the probability that their value of is smaller than the determined value by chance is less than 32%): for Kepler305 there is a ~15% chance that this particular system lies this close to resonance by chance, while for Kepler31 the probability is ~ 20%.
We should remark that these specific values for the probabilities that each system is this close to exact resonance by chance are calibrated by the choice 0 < Δ < 0.05 for both pairs of planets, which is used in generating the fictional systems. This value is however not arbitrary. For, it must be consistent with the choice made in Sect. 3.1, which produced only these three systems with both the inner and outer planet pair this close to first order mean motion resonance: there, the choice Δ < max Δ = 0.05 was dictated by the observation of the location of the peak wide of nominal resonance (Hadden & Lithwick 2017 and our Fig. 3), so restricting the interval in Δ values with smaller max Δ might have resulted in excluding potential systems. On the other hand, increasing max Δ in the generation of the fictional systems would have the only effect to decrease the calculated probabilities. Therefore, we conclude that our choice of max Δ = 0.05 is not arbitrary, and gives a reasonable upper bound to the probabilities that each system finds itself so close to exact resonance by pure chance.
For completeness, we report the observed variation of the threebody Laplace angle in these systems, since its libration can be in principle a sufficient condition to conclude that they are indeed resonant. For Kepler305 we checked that the Laplace angle φ_{L} = 2λ_{1} − 4λ_{2} + 2λ_{3} satisfies given the observed transits periods; for Kepler31, φ_{L} = λ_{1} − 3λ_{2} + 2λ_{3} satisfies ; finally for YZ Cet, φ_{L} = 2λ_{1} − 5λ_{2} + 3λ_{3} satisfies . As we argued in the Introduction, in case of libration of the resonant angles around the equilibrium point (and hence libration of the Laplace angle) the average of the (a, e) oscillation falls on the equilibrium point, while in case of circulation, the average falls off the equilibrium point curve. Consequently, the circulation of the Laplace angle implies that the libration amplitude is larger than the distance of the equilibrium point from the axis e = 0, and that cannot be zero. However, this does not mean that the system did not reach that point via divergent migration: being the (a, e) equilibrium so close to to e = 0, even a minute perturbation can induce circulation of the Laplace angle. Hence the libration of the Laplace angle is a sufficient but not necessary condition to conclude that a system’s dynamical history has been shaped by resonant capture and subsequent resonant repulsion driven by dissipation.
Fig. 9 Cumulative distribution functions for for randomly generated systems close to chains of any possible combinations of the 2:1 and 3:2 mean motion resonances. We indicate the chains that represent selected systems from Fig. 4; for each of them, a point indicates the observed . From this, we obtain on the vertical axis the probability that these systems could have this value of by pure chance. 
4.2 The 5:4 – 4:3 resonant chain on Kepler60 and other nearresonant systems with k > 3
While in this work we have concentrated on the 2:1 and 3:2 mean motion commensurabilities, it is worthwhile to point out that more compact resonant chains are possible, and Kepler60 represents a notable example. This system hosts three planets with mean observed periods of ≃5.49, ≃ 8.29 and ≃ 16.74 days respectively. Their masses have been constrained via TTV to be ≃4 M_{⊕} for all planets (JontofHutter et al. 2016). The mean motions of the planets satisfy , hinting at a resonant configuration. Indeed, Goździewski et al. (2016) found that the TTV signal for these planets is consistent with a true threebody Laplacelike resonance as well as a chain of 5:4–4:3 twobody mean motion resonances. Using the system’s parameters we can find a resonant equilibrium configuration as in Sect. 3.2, by imposing a_{3} ∕a_{1} to be equal to the observed . This gives of order 4 × 10^{−5}. Using an analogous argument to that of Fig. 9, we find that there is only a 0.25% probability that Kepler60 lies this close to a 5:4–4:3 resonant chain by pure chance. The eccentricities that we find at the selected resonant equilibrium point are of order e_{1} ≃ 0.02, e_{2} ≃ 0.03, e_{3} ≃ 0.01 for the observed planetary masses. These numbers are quite close to the ones consistent for the twobody mean motion resonance chain solution found in Goździewski et al. (2016). Note in passing that their solution is for nonvanishing libration amplitude of the four resonant angles (however their mean values are the same found here for a stable resonant equilibrium).
For completeness, we cite other nearresonant systems of three planets with k > 3 that are found in the catalogue. The only ones which satisfy our criterion Δ < 0.05 for both pairs are K2239 (close to a 3:2–4:3 chain), Kepler289 (close to a 2:1–2:1 chain), Kepler226 (close to a 4:3–3:2 chain) and Kepler431 (close to a 5:4–4:3). Of these, only the latter two satisfy Δ > 0 for both pairs.
5 Conclusions
In this work, we have generalised the formalism of dissipative divergence of resonant orbits to multiresonant chains. The analyticalstudy performed in Sect. 2 allows us to predict the orbital configurations of systems of planets deep in a chain of first order mean motion resonances, and therefore, even though at a lesser degree of precision, of systems that are in resonance with a finite amplitude of libration. Then, we showed in Sect. 2.2 that under the effect of slow dissipation a nearlyresonant system is expected to follow the loci of the equilibrium points of the resonant Hamiltonian (Eq. (13)) maintaining the amplitude of libration in an adiabatic manner. Therefore, if the orbital architecture of a system is found near one of these equilibrium points, it is strongly suggestive that the envisioned scenario of slow convergent orbital migration leading to capture into resonance and subsequent orbital divergence due to dissipative evolution really occurred for the system. In the light of the results presented above, we can draw the following conclusions.
On the one hand, we must face the fact that the orbital architecture of a significant fraction of the systems of three planets is actually not consistent with these physical mechanisms. More precisely, the majority of the systems are not close to resonance, implying that either they never captured in resonance in the first place, or they escaped from resonance due to a violent instability (Izidoro et al. 2017) losing any memory of their resonant dynamical past. To ponder these two possibilities, consider first of all that some form of orbital transport is expected to take place: for example, the majority of planets with R_{pl} ≳ 1.6 R_{⊕} have H/He gaseous atmospheres that cannot be explained by production of volatiles after the formation of the planet (Rogers 2015), implying that these planets formed while the protoplanetary disc was still present. The associated planet–disc interaction would then force the planets’ orbital elements to change, in other words, force the planets to migrate. However, orbital migration may not be convergent (Migaszewski 2015, 2016), that is, not leading to resonant capture. Moreover, some mechanisms have been proposed to inhibit the capture even in the case of convergent migration, such as turbulence in the disc or edependent migration rates. Nevertheless, these processes alone do not adequately explain the lack of resonance in the exoplanet sample (e.g. Batygin & Adams 2017; Deck & Batygin 2015; Xu et al. 2018). For these reasons, it is more likely that capture into mean motion resonance is a common outcome of the early epochs of discplanet interaction, but the subsequent evolution after the disc removal is subject to instabilities which break the compact configuration. This approach seems to be able to reproduce the observed distribution of period ratios if these instabilities are extremely common (Izidoro et al. 2017). The primary mechanism through which planets are ejected from resonance, however, remains elusive, and a topic of active research.
On the other hand, systems with orbital properties that are compatible with a (near) resonant state do exist in the exoplanet census. These include the already known examples mentioned above of Trappist1, Kepler223, GJ867 and Kepler60, and, potentially, some of the systems analysed in this paper. That is, while it is difficult to prove definitively that a given system is now in resonance in a formal sense (the resonant angles are in libration), in this work we have developed a method to quantitatively test the consistency of a given orbital architecture with a dynamical history characterised by resonant capture and subsequent dissipative evolution. This is achieved through the calculation of the quantity defined in Eq. (16), which is obtained directly from the observed semimajor axis ratios and, as we have shown, depends very weakly on the mass ratios between the planets, making the observational uncertainties on the latter quantities irrelevant. Then, using the approach illustrated in Fig. 9, this “indicator” can be turned into a quantitative probability that the system is related to the considered resonant chain. In this sense, we have found that there is a ~ 15, ~ 20 and ~ 80% probability that Kepler305, Kepler31 and YZ Cet respectively find themselves this close to resonance merely by chance. Multiplying these probabilities we find that there is only a ~2% probability that all three of these systems lie close to resonance just by chance. This suggests strongly that at least some of them should have had a resonant dynamical history. Although the sample is clearly too small to make any meaningful inference, the probabilities of resonant association that we have found indicate that between 1/3 and 2/3 of the systems with 0 < Δ < 0.05 show memory of the processes of resonant capture; this is consistent with the histogram of Fig. 3, where the peak wide of the resonance is about 2 times higher than the underlying randomlike flat distribution.
The architecture of many planetary systems observed by transit is not well constrained by observations. Opportunities for more extensive characterisation will come from missions such as the Transiting Exoplanet Survey Satellite (TESS) or the PLAnetary Transits and Oscillations of stars (PLATO), which are designed to target bright stars to allow for followup via further groundbased and spacebased observations (with methods such as radial velocity). This will allow for a better quantification of planetary masses, radii, ages of the systems and eccentricities. In the light of this augmented perception that we can expect to acquire, our study outlines the groundwork for further dynamical characterisation of the physical processes that shaped the presentday architectures of extrasolar planetary systems.
Appendix A Reduced Hamiltonian to a common planetary mass factor
In the course of the paper we make implicit use of a reduced Hamiltonian which incorporates the planetary masses through a common planettostar mass factor . In this appendix, we detail the construction of this Hamiltonian and its use in the paper.
Consider three planets, whose physical parameters are labelled 1, 2, and 3 for the inner, middle and outer planet respectively, orbiting around a star of mass M_{*} on the same plane. Suppose that the planets are (close to) a chain of meanmotion resonance, with nominal semimajor axes ā and that the deviations of the semimajor axes from the nominal values are small, and assume that the eccentricities are small enough, so that an analysis to first order in e is valid. These are the working assumptions throughout Sect. 2. Having fixed the planet–planet mass ratios m_{1} ∕m_{2} = β_{1} and m_{2} ∕m_{3} = β_{2}, we introduce the average planet–star mass ratio (A.1)
Inverting this expression we easily get all planetary masses in terms of , (A.2)
with coefficients c depending on M_{*}, β_{1} and β_{2} only.
We introduce the modified Delaunay actionangle variables (Λ, Γ, λ, γ) as in Eq. (5), but we rescale the actions by the common mass factor : this gives the following definition for the Λ’s (we maintain the same notation as the nonrescaled actions for simplicity) (A.3)
and the same formal definition of Γ = Λe^{2}∕2 at lowest order in e. We now introduce the Hamiltonian for the problem, that is the sum of the Keplerian Hamiltonian (Eq. (6)) and the resonant interaction Hamiltonian (Eq. (7)) to first order in e. However, since we have rescaled the actions by , in order for the Hamilton equations to be conserved we must also rescale the Hamiltonian itself by . As in Sect. 2, since we are not considering large deviations in the semimajor axes from their nominal values, and since the resonant Hamiltonian is already of order , we evaluate the resonant Hamiltonianon the nominal values defined fromā using Eq. (A.3). It is then easy to see that the rescaled Keplerian Hamiltonian takes the form (A.4)
and is therefore independent of , while the rescaled resonant part will have a multiplicative coefficient : (A.5)
From here, the sequence of changes of variables detailed in Sect. 2 can be performed using the same formal definitions for the new rescaled variables.
Already from the Hamiltonian written in terms of the rescaled variables Λ and Γ ∝ e^{2} one can see the following. Assuming a fixed equilibrium value of the semimajor axes (that is, of the Λ’s), a change in the planetstar mass factor will have the only effect to rescale the equilibrium values of all by the same quantity. In this configuration, the equilibria of the semimajor axis ratios and remain independent of .
References
 Batygin, K. 2015, MNRAS, 451, 2589 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Adams, F. C. 2017, AJ, 153, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Morbidelli, A. 2013a, AJ, 145, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Morbidelli, A. 2013b, A&A, 556, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Batygin, K., Deck, K. M., & Holman, M. J. 2015, AJ, 149, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Beaugé, C., Michtchenko, T. A., & FerrazMello, S. 2006, MNRAS, 365, 1160 [NASA ADS] [CrossRef] [Google Scholar]
 Cresswell, P., & Nelson, R. P. 2006, A&A, 450, 833 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cresswell, P., & Nelson, R. P. 2008, A&A, 482, 677 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deck, K. M., & Batygin, K. 2015, ApJ, 810, 119 [Google Scholar]
 Deck, K. M., Payne, M., & Holman, M. J. 2013, ApJ, 774, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Delisle, J.B., & Laskar, J. 2014, A&A, 570, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J.B., Laskar, J., Correia, A. C. M., & Boué, G. 2012, A&A, 546, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J.B., Correia, A. C. M., & Laskar, J. 2015, A&A, 579, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Goldreich, P., & Schlichting, H. E. 2014, AJ, 147, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Goździewski,K., Migaszewski, C., Panichi, F., & Szuszkiewicz, E. 2016, MNRAS, 455, L104 [Google Scholar]
 Hadden, S., & Lithwick, Y. 2014, ApJ, 787, 80 [Google Scholar]
 Hadden, S., & Lithwick, Y. 2017, AJ, 154, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Henrard, J., & Lamaitre, A. 1983, Celest. Mech., 30, 197 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750 [NASA ADS] [CrossRef] [Google Scholar]
 JontofHutter,D., Ford, E. B., Rowe, J. F., et al. 2016, ApJ, 820, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Lithwick, Y., & Wu, Y. 2012, ApJ, 756, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nat. Astron., 1, 0129 [Google Scholar]
 Mardling, R. A., & Lin, D. N. C. 2002, ApJ, 573, 829 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Michtchenko, T. A., Beaugé, C., & FerrazMello, S. 2008, MNRAS, 387, 747 [NASA ADS] [CrossRef] [Google Scholar]
 Millholland, S., Wang, S., & Laughlin, G. 2017, ApJ, 849, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Migaszewski, C. 2015, MNRAS, 453, 1632 [NASA ADS] [CrossRef] [Google Scholar]
 Migaszewski, C. 2016, MNRAS, 458, 2051 [Google Scholar]
 Morbidelli, A. 2002, Modern Celestial Mechanics: Aspects of Solar System Dynamics (London: Taylor & Francis) [Google Scholar]
 Murray, C. D., & Dermott, S. F. 2000, Solar System Dynamics (Cambridge, UK: Cambridge University Press) [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Terquem, C. 2010, MNRAS, 405, 573 [NASA ADS] [Google Scholar]
 Pichierri, G., Morbidelli, A., & Crida, A. 2018, Celest. Mech. Dyn. Astron., 130, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Poincare, H. 1892, Les Méthodes nouvelles de la mécanique céleste (Paris: GauthierVillars et fils) [Google Scholar]
 Rivera, E. J., Laughlin, G., Butler, R. P., et al. 2010, ApJ, 719, 890 [Google Scholar]
 Rogers, L. A. 2015, ApJ, 801, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Tamayo, D., Rein, H., Petrovich, C., & Murray, N. 2017, ApJ, 840, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Wisdom, J. 1980, AJ, 85, 1122 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y., & Lithwick, Y. 2013, ApJ, 772, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, W., Lai, D., & Morbidelli, A. 2018, MNRAS, 481, 1538 [NASA ADS] [CrossRef] [Google Scholar]
We note that since the publication of the aforementioned paper, more planets have been observed in the same system, including two confirmed planets HD 40307 f and HD 40307 g. For this reason, we will not consider this system in the current work, although we draw inspiration from the analysis of Papaloizou & Terquem (2010).
We note that even away from nominal resonance all four resonant angles can continue to librate when the system is sufficiently close to the resonant equilibrium point, unlike what has been erroneously stated in Sect. 4 of Batygin & Morbidelli (2013a).
At higher eccentricities the main term which might shift the equilibria in a_{2}∕a_{1}, a_{3}∕a_{2} to the left of exact resonance is the second order (secular) term which describes the interaction between the inner planet and outer planet; however, we checked that adding this term to the Hamiltonian, even at high eccentricities and for a very massive outer planet the picture does not change.
We should remark however that even the periods are not known with arbitrary precision, meaning that there might be small discrepancies in the values that are used in different works. In this paper, we will use the ones listed in the NASA Exoplanet Archive catalogue without considering error bars. This is enough for the scope of our analysis.
All Figures
Fig. 1 Observations of planethosting stars reveal that multiplanetary systems are not rare, hosting over 1600 confirmed planets (panel a). The period ratio distribution of neighbouring planets is shown in panel b. One can observe an overall broad distribution as well as a number of peaks slightly wide of resonant ratios, especially the 2:1 and 3:2 commensurabilities. Data was obtained from the Nasa Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu/). Panel a: histogram of multiplanetary system by number of planets in each system. Panel b: distribution of period ratios of neighbouring planets in exoplanetary systems. 

In the text 
Fig. 2 Equilibrium curves showing the loci of the stable resonant equilibrium points calculated as explained in the text, in the case of a 3:2 – 3:2 mean motion resonance chain, with m_{1} = m_{2} = m_{3} = 10^{−4}M_{*}. The full curves are calculated using the expanded Keplerian Hamiltonian (Eq. (11)), while the dashed curves are calculated using the unexpanded Keplerian Hamiltonian (Eq. (6)), showing very little difference down to very small eccentricities and for reasonable values of the nearly exactly resonant semimajor axis ratio. The location of the nominal resonant semimajor axis ratio (3/2)^{2/3} is shown bya vertical orange line. We also superimpose the numerically computed evolution of a threeplanet system deep in the 3:2 mean motion resonance (for both pairs) and undergoing dissipative evolution depicted with transparent lines: the system follows the locations of the equilibrium points, which are close to the curves calculated analytically. 

In the text 
Fig. 3 Distribution of the (signed) normalised distance from first order mean motion resonances with k =2, 3 in all exoplanetary systems (yellow) and for threeplanet systems (purple). We note a clear peak to the right of the value Δ = 0 (corresponding to the exact nominal resonance, indicated by a vertical orange line), which is most prominent for 0 < Δ ≲0.05. Out of the 358 pairs plotted here, there are 123 total pairs in this configuration. For the threeplanets systems only, 48 pairs have 0 <Δ ≲ 0.05, out of the 121 shown in the purple histogram. 

In the text 
Fig. 4 Orrery of the threeplanets systems sufficiently close to first order mean motion resonances k : k − 1 with k =2, 3, with a normalised distance to resonance Δ < 0.05 for both pairs. For each system, we place a circle in correspondence of the period of each planet, and indicate between pairs of planets the first order mean motion resonance in which they are envisioned to reside (below) and the normalised distance to that resonance Δ (above); the sign of Δ indicates if the pair of planets are narrow (Δ < 0) or wide (Δ > 0) of the resonance. For our analysis, we will only consider systems for which Δ > 0 for both pairs (the systems enclosed by a box). The size of the circle is an indication of the estimated radius of the planet (the small dot in the top right corner demonstrates the size of Earth). For YZ Ceti this information is not available, but this does not pose a problem for our study. 

In the text 
Fig. 5 Locations of the resonant equilibrium points in the a_{i+1}∕a_{i} vs. planes, i = 1, 2, for three equalmass planets in a 3:2 – 2:1 mean motion resonance chain, close to which Kepler305 resides. Orange vertical lines show the exact nominal commensurability, while dashed vertical lines show the observed a_{2} ∕a_{1} and a_{3} ∕a_{2} in the case of Kepler305. As explained in the text, we select one equilibrium configuration (indicated by the red dot in both panels) by requiring that (a_{3} ∕a_{1})_{eq} ≡ (a_{3}∕a_{1})_{obs}, which automatically fixes all orbital elements a_{2}∕a_{1}, a_{3} ∕a_{2}, e_{1}, e_{2} and e_{3} along the equilibrium curves. 

In the text 
Fig. 6 Top row: analytical maps constructed for Kepler305 as explained in Sect. 3.2. Panel a: we plot the quantity , which represents how close the system is now to some resonant equilibrium point, for different mass ratios m_{1} ∕m_{2} and m_{2} ∕m_{3} (each point in this plot is constructed by repeating the procedure described in Fig. 5). We notice that changes verylittle with the mass ratios, and is of the order of ~0.002. Comparing with Fig. 9, we see that this can be the case by pure chance only in ~ 15% of randomly generated systems close to the 3:2–2:1 mean motion resonance chain. Panel b: we show a map of the quantity representing the equilibrium orbital configuration that is selected at each fixed value of m_{1} ∕m_{2} and m_{2} ∕m_{3} by imposing (a_{3}∕a_{1})_{eq} ≡ (a_{3}∕a_{1})_{obs}. Bottom row: numerical maps constructed for Kepler305 as explained in Sect. 3.3. We show numerical maps of in panel c and in panel d, analogous to the analytical plots above (over a slightly smaller range of mass ratios for simplicity). These are intended to validate the analytical maps, and show very good agreement between corresponding panels. 

In the text 
Fig. 7 Same as in Fig. 6, but for the system YZ Cet, residing close to a 3:2–3:2 chain. The value of across all planetary mass ratios can be matched against the corresponding curve in Fig. 9, where we find that ~ 80% of randomlygenerated systems lie this close to the 3:2–3:2 chain. 

In the text 
Fig. 8 Same as in Fig. 6, but for the system Kepler31. Only analytical maps are shown in this case since in some simulations capture into resonance was unsuccessful due to overstability of the captured state for different planetary mass ratios. As explained in the text, this issue is modeldependent and is not within the scope of our analysis. Moreover, inthe simulations where capture was successful, the results agree well with the analytical calculations, showing . Comparing with Fig. 9 we see that there is a ~20% probability that Kepler31 lies this close to the 2:1 – 2:1 chain by pure chance. 

In the text 
Fig. 9 Cumulative distribution functions for for randomly generated systems close to chains of any possible combinations of the 2:1 and 3:2 mean motion resonances. We indicate the chains that represent selected systems from Fig. 4; for each of them, a point indicates the observed . From this, we obtain on the vertical axis the probability that these systems could have this value of by pure chance. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.