The role of dissipative evolution for three-planet, near-resonant extrasolar systems

TL;DR

This paper extends the understanding of how dissipative forces influence the orbital evolution of three-planet systems near mean-motion resonances, revealing that additional physical processes are likely involved in shaping observed architectures.

Contribution

It generalizes the dissipative divergence framework from three-body to multi-resonant chains and applies it to observed three-planet systems, highlighting limitations of current models.

Findings

Some systems' architectures are consistent with resonant locking followed by dissipation.

Most systems' configurations cannot be explained by simple migration and dissipation alone.

Additional physical mechanisms are likely important in early planetary system evolution.

Abstract

Early dynamical evolution of close-in 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 semi-major axes of initially resonant pairs of planets will gradually diverge under the influence of long-term 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 close-in extrasolar planets, the existing theoretical picture is limited to the specific case of the planetary three-body problem. In this study, we generalise the framework of dissipative divergence of resonant orbits to multi-resonant chains, and apply our results to the current observational census of well-characterised three-planet systems. Focusing on the 2:1 and 3:2 commensurabilities, we identify three 3-planet systems, whose current orbital architecture is consistent with an evolutionary history wherein convergent migration first locks the planets into a multi-resonant configuration and subsequent dissipation repels the orbits away from exact commensurability. Nevertheless, we find that the architecture of the overall sample of multi planetary 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.

Equations

$$\mathcal{H}={\mathcal{H}}_{kepl}+{\mathcal{H}}_{pert},$$

$${\mathcal{H}}_{kepl}=-\frac{\mathcal{G}{M}_{\ast }{m}_1}{2a_1}-\frac{\mathcal{G}{M}_{\ast }{m}_2}{2a_2}-\frac{\mathcal{G}{M}_{\ast }{m}_3}{2a_3},$$

$$\begin{split}\mathcal{H}_{res}&=-\frac{\mathcal{G}m_{1}m_{2}}{a_{2}}\left(f_{res}^{(1,in)}e_{1}\cos\big{(}k^{in}\lambda_{2}-(k^{in}-1)\lambda_{1}-\varpi_{1}\big{)}+\right.\\ &\quad\quad\quad\left.+f_{res}^{(2,in)}e_{2}\cos\big{(}k^{in}\lambda_{2}-(k^{in}-1)\lambda_{1}-\varpi_{2}\big{)}\right)+\\ &-\frac{\mathcal{G}m_{2}m_{3}}{a_{3}}\left(f_{res}^{(1,out)}e_{2}\cos\big{(}k^{out}\lambda_{3}-(k^{out}-1)\lambda_{2}-\varpi_{2}\big{)}+\right.\\ &\quad\quad\quad\left.+f_{res}^{(2,out)}e_{3}\cos\big{(}k^{out}\lambda_{3}-(k^{out}-1)\lambda_{2}-\varpi_{3}\big{)}\right),\end{split}$$

$$\bar{\mathcal{H}}={\mathcal{H}}_{kepl}+{\mathcal{H}}_{res}+\mathcal{O}(e^2,I^2),$$

$$\Lambda$$

$$\Gamma$$

$${\mathcal{H}}_{kepl}=-\sum _{i=1}^3\frac{{\mathcal{G}}^2(M_{\ast }+m_i{)}^2{\mu }_i^3}{2{\Lambda }_i^2}\simeq -\sum _{i=1}^3\frac{m_i^3}{2}{\left(\frac{\mathcal{G}{M}_{\ast }}{{\Lambda }_i}\right)}^2,$$

$$\begin{split}\mathcal{H}_{res}&=-\frac{\mathcal{G}^{2}M_{*}m_{1}m_{2}^{3}}{\Lambda_{2}^{2}}\\ &\quad\quad\times\left(f_{res}^{(1,in)}\sqrt{\frac{2\Gamma_{1}}{\Lambda_{1}}}\cos\big{(}k^{in}\lambda_{2}-(k^{in}-1)\lambda_{1}+\gamma_{1}\big{)}\right.\\ &\quad\quad\quad\left.+f_{res}^{(2,in)}\sqrt{\frac{2\Gamma_{2}}{\Lambda_{2}}}\cos\big{(}k^{in}\lambda_{2}-(k^{in}-1)\lambda_{1}+\gamma_{2}\big{)}\right)+\\ &\quad-\frac{\mathcal{G}^{2}M_{*}m_{2}m_{3}^{3}}{\Lambda_{3}^{2}}\\ &\quad\quad\times\left(f_{res}^{(1,out)}\sqrt{\frac{2\Gamma_{2}}{\Lambda_{2}}}\cos\big{(}k^{out}\lambda_{3}-(k^{out}-1)\lambda_{2}+\gamma_{2}\big{)}\right.\\ &\quad\quad\quad\left.+f_{res}^{(2,out)}\sqrt{\frac{2\Gamma_{3}}{\Lambda_{3}}}\cos\big{(}k^{out}\lambda_{3}-(k^{out}-1)\lambda_{2}+\gamma_{3}\big{)}\right);\end{split}$$

$$\begin{array}{rl}\mathcal{K}& ={\Lambda }_1+\frac{k^{in}-1}{k^{in}}{\Lambda }_2+\frac{(k^{in}-1)(k^{out}-1)}{k^{in}{k}^{out}}{\Lambda }_3,\kappa ={\lambda }_1,\\ {\Theta }^{(1)}& =\frac{1}{k^{in}}{\Lambda }_2+\frac{k^{out}-1}{k^{in}{k}^{out}}{\Lambda }_3,{\theta }^{(1)}=k^{in}{\lambda }_2-(k^{in}-1){\lambda }_1,\\ {\Theta }^{(2)}& =\frac{1}{k^{out}}{\Lambda }_3,{\theta }^{(2)}=k^{out}{\lambda }_3-(k^{out}-1){\lambda }_2;\end{array}$$

$$\begin{array}{rl}\frac{\mathcal{K}}{{\Lambda }_3}& =\frac{{\mu }_1}{{\mu }_3}\sqrt{\frac{M_{\ast }+m_1}{M_{\ast }+m_3}\frac{a_1}{a_3}}+\frac{k^{in}-1}{k^{in}}\frac{{\mu }_2}{{\mu }_3}\sqrt{\frac{M_{\ast }+m_2}{M_{\ast }+m_3}\frac{a_2}{a_3}}+\\ & +\frac{(k^{in}-1)(k^{out}-1)}{k^{in}{k}^{out}},\end{array}$$

$$\begin{array}{rl}{\mathcal{H}}_{kepl}& =-\sum _{i=1}^3\frac{{\mathcal{G}}^2(M_{\ast }+m_i{)}^2{\mu }_i^3}{2}\\ & \times \left(\frac{1}{{\bar{\Lambda }}_i^2}-2\frac{1}{{\bar{\Lambda }}_i^3}\delta {\Lambda }_i+3\frac{1}{{\bar{\Lambda }}_i^4}\delta {\Lambda }_i^2+\mathcal{O}(\delta {\Lambda }_i^3)\right),\end{array}$$

$${\mathcal{H}}_{kepl}=\sum _{i=1}^3\left(4{\bar{n}}_i{\Lambda }_i-\frac{3}{2}{\bar{h}}_1{\Lambda }_i^2\right),$$

$$\Omega$$

$${\Psi }_1^{(1)}$$

$${\Psi }_1^{(2)}$$

$${\Psi }_2^{(1)}$$

$${\Psi }_2^{(2)}$$

$$\begin{array}{rl}\bar{\mathcal{H}}& ={\mathcal{H}}_{kepl}+{\mathcal{H}}_{res},\\ {\mathcal{H}}_{kepl}& ={\mathcal{H}}_{kepl}\left({\Psi }_1^{(1)},{\Psi }_1^{(2)};\mathcal{K},\Omega \right),\\ {\mathcal{H}}_{res}& ={\mathcal{H}}_{res}\left({\Psi }_1^{(1)},{\Psi }_1^{(2)},{\Psi }_2^{(1)},{\Psi }_2^{(2)},{\psi }_1^{(1)},{\psi }_1^{(2)},\delta {\gamma }^{(1)},\delta {\gamma }^{(2)};\mathcal{K},\Omega \right),\end{array}$$

$$\begin{array}{rl}\frac{\partial \bar{\mathcal{H}}}{\partial {\Psi }_1^{(1)}}& =0,\text{leavevmode}\frac{\partial \bar{\mathcal{H}}}{\partial {\Psi }_1^{(2)}}=0,\text{leavevmode}\frac{\partial \bar{\mathcal{H}}}{\partial {\Psi }_2^{(1)}}=0,\text{leavevmode}\frac{\partial \bar{\mathcal{H}}}{\partial {\Psi }_2^{(2)}}=0,\\ \frac{\partial \bar{\mathcal{H}}}{\partial {\psi }_1^{(1)}}& =0,\text{leavevmode}\frac{\partial \bar{\mathcal{H}}}{\partial {\psi }_1^{(2)}}=0,\text{leavevmode}\frac{\partial \bar{\mathcal{H}}}{\delta {\gamma }^{(1)}}=0,\text{leavevmode}\frac{\partial \bar{\mathcal{H}}}{\delta {\gamma }^{(2)}}=0.\end{array}$$

$${\Delta }_{k/(k-1)}:=\frac{k-1}{k}\frac{P_2}{P_1}-1,$$

$$\bar{\delta }(a_3/a_2):=\frac{({(a_3/a_2)|}_{eq}-{(a_3/a_2)|}_{obs})}{{(a_3/a_2)|}_{obs}}$$

$$\tilde{m}=\frac{m_1+m_2+m_3}{3M_{\ast }}=\frac{m_1(1+{\beta }_1^{-1}+{\beta }_1^{-1}{\beta }_2^{-1})}{3M_{\ast }}.$$

$$\begin{array}{r}{m}_1=c_1\tilde{m}:=\frac{3{\beta }_1{\beta }_2{M}_{\ast }}{1+{\beta }_2+{\beta }_1{\beta }_2}\tilde{m},\\ m_2=c_2\tilde{m}:=\frac{3{\beta }_1{M}_{\ast }}{1+{\beta }_2+{\beta }_1{\beta }_2}\tilde{m},\\ m_3=c_3\tilde{m}:=\frac{3M_{\ast }}{1+{\beta }_2+{\beta }_1{\beta }_2}\tilde{m},\end{array}$$

$$\begin{array}{r}{\Lambda }_1=\frac{3{\beta }_1{\beta }_2{M}_{\ast }}{1+{\beta }_2+{\beta }_1{\beta }_2}\sqrt{\mathcal{G}{M}_{\ast }{a}_1}=c_1\sqrt{\mathcal{G}{M}_{\ast }{a}_1},\\ {\Lambda }_2=\frac{3{\beta }_1{M}_{\ast }}{1+{\beta }_2+{\beta }_1{\beta }_2}\sqrt{\mathcal{G}{M}_{\ast }{a}_2}=c_2\sqrt{\mathcal{G}{M}_{\ast }{a}_2},\\ {\Lambda }_3=\frac{3M_{\ast }}{1+{\beta }_2+{\beta }_1{\beta }_2}\sqrt{\mathcal{G}{M}_{\ast }{a}_3}=c_3\sqrt{\mathcal{G}{M}_{\ast }{a}_3},\end{array}$$

$${\mathcal{H}}_{kepl}=-\sum _{i=1}^3\frac{c_i^3}{2}{\left(\frac{\mathcal{G}{M}_{\ast }}{{\Lambda }_i}\right)}^2,$$

$$\begin{split}\mathcal{H}_{res}&=\tilde{m}\left[-\frac{\mathcal{G}^{2}M_{*}c_{1}c_{2}^{3}}{\bar{\Lambda}_{2}^{2}}\right.\\ &\quad\quad\times\left(f_{res}^{(1,in)}\sqrt{\frac{2\Gamma_{1}}{\bar{\Lambda}_{1}}}\cos\big{(}k^{in}\lambda_{2}-(k^{in}-1)\lambda_{1}+\gamma_{1}\big{)}\right.\\ &\quad\quad\quad\left.+f_{res}^{(2,in)}\sqrt{\frac{2\Gamma_{2}}{\bar{\Lambda}_{2}}}\cos\big{(}k^{in}\lambda_{2}-(k^{in}-1)\lambda_{1}+\gamma_{2}\big{)}\right)+\\ &\quad-\frac{\mathcal{G}^{2}M_{*}c_{2}c_{3}^{3}}{\bar{\Lambda}_{3}^{2}}\\ &\quad\quad\times\left(f_{res}^{(1,out)}\sqrt{\frac{2\Gamma_{2}}{\bar{\Lambda}_{2}}}\cos\big{(}k^{out}\lambda_{3}-(k^{out}-1)\lambda_{2}+\gamma_{2}\big{)}\right.\\ &\quad\quad\quad\left.\left.+f_{res}^{(2,out)}\sqrt{\frac{2\Gamma_{3}}{\bar{\Lambda}_{3}}}\cos\big{(}k^{out}\lambda_{3}-(k^{out}-1)\lambda_{2}+\gamma_{3}\big{)}\right)\right];\end{split}$$

Full text

11institutetext: Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, France

11email: gabriele.pichierri@oca.eu 22institutetext: Division of Geological and Planetary Sciences, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

The role of dissipative evolution for three-planet, near-resonant extrasolar systems

Gabriele Pichierri 11   

Konstantin Batygin 22   

Alessandro Morbidelli 11

Early dynamical evolution of close-in 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 semi-major axes of initially resonant pairs of planets will gradually diverge under the influence of long-term 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 close-in extrasolar planets, the existing theoretical picture is limited to the specific case of the planetary three-body problem. In this study, we generalise the framework of dissipative divergence of resonant orbits to multi-resonant chains, and apply our results to the current observational census of well-characterised three-planet systems. Focusing on the 2:1 and 3:2 commensurabilities, we identify three 3-planet systems, whose current orbital architecture is consistent with an evolutionary history wherein convergent migration first locks the planets into a multi-resonant configuration and subsequent dissipation repels the orbits away from exact commensurability. Nevertheless, we find that the architecture of the overall sample of multi planetary 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.

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 multi-planets 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 near-resonant 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 low-integer ratios, and a definite excess just wide of these values, especially the 2:1 and 3:2, see Figure 1.

Numerical simulations show that compact chains of mean motion resonances are a common outcome of slow, convergent orbital transport of planets within protoplanetary disks. Although details of disk-driven 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 Trappist-1, a star famously hosting 7 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 $N$-body 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 sub-Neptune planets of Kepler-223 (Mills et al., 2016) and the now-classic example of Laplace-like resonance in GJ-876 (Rivera et al., 2010; Batygin et al., 2015).

In light of the fact that convergent migration should lock planets into mean motion commensurability, how can we 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 $a_{1}/a_{2}=((k-1)/k)^{2/3}$ (where $a_{1}$ and $a_{2}$ are the semi-major 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 disk 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 will see in Section 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 for such planets captured in resonance via migration in protoplanetary disks. These observations can potentially be interpreted as evidence for dissipative processes acting on the planetary systems after the disk phase.

Papaloizou & Terquem (2010) considered the specific case of the K-dwarf HD 40307, which hosts111 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). 3 hot super-Earths/mini-Neptunes 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 (2013), 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. Note that, for two planets, this repulsion can be easily understood if one considers that any process that dissipates the energy, $E\propto-m_{1}/a_{1}-m_{2}/a_{2}$, and at the same time conserves angular momentum, $L\propto m_{1}\sqrt{a_{1}}+m_{2}\sqrt{a_{2}}$, 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 (Figure 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 disk in which they formed; they interact with the disk, which ultimately results in a resonant capture. Then, after the disk 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 multi-systems migrating in a disk of gas. One approach towards a better approximation would be to perform full hydrodynamic simulations of planets immersed in their protoplanetary disk accounting for various disk parameters (such as disk 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 disk 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 multi-planets systems mentioned above.

In this paper we focus on slow convergent Type-I migration in a disk of gas, and adopt simple synthetic analytical formulæ for the work and the torque generated by the disk on the planets (Cresswell & Nelson 2006, 2008); the requirement that exact prescriptions for the interaction between the planets and the disk 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 post-disk 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 $\dot{E}<0$ and $\dot{\mathcal{L}}=0$. We conclude that our specific implementation of Type-I migration and tidal dissipation after the disk removal are therefore not restrictive, which makes our results generalisable to any other equivalent processes. In the light of these considerations, we ask the following question: assuming that planets are captured into resonance and undergo dissipative evolution after the disk phase, is it possible to reproduce the observed orbital configuration of real exoplanetary systems which reside close to resonance? In other words: are the aforementioned physical processes compatible with the distribution of near-resonant period ratios that emerges from available data?

1.1 Methods and physical setup

In order to answer this question, we examined the NASA Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu) and selected confirmed 3-planet 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 multi-resonant chain, which would be indicative of a dynamical history characterised by the physical processes described above.

Evidence suggesting that planets around Trappist-1 and Kepler-223 truly reside in a resonant configuration has recently been marshalled from the observed libration of the three-body Laplace angles. To this end, recall that if two neighbouring pairs of planets in a multi-planet systems are in the $k^{(in)}$:$(k^{(in)}-1)$ and $k^{(out)}$:$(k^{(out)}-1)$ resonances respectively (so that the resonant angles $k^{(in)}\lambda_{2}-(k^{(in)}-1)\lambda_{1}-\varpi_{2}$ and $k^{(out)}\lambda_{3}-(k^{(out)}-1)\lambda_{2}-\varpi_{2}$ are librating), then the Laplance angle $\varphi_{L}=(k^{(in)}-1)\lambda_{1}-(k^{(in)}+k^{(out)}-1)\lambda_{2}+k^{(out)}\lambda_{3}$ will be automatically librating as well. The advantage of examining this three-body angle over the two-body resonant angles is that the latter contain the longitudes of the pericenters $\varpi$, whose precession rates are poorly constrained by the data, while the former only includes the mean longitudes $\lambda$ 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 three-body 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 first-order 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 semi-major axis ratios. It is worth anticipating here the following important point. As it will be clear later (see Section 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 $\mathcal{L}$ a free parameter: it is therefore always possible to find a value of $\mathcal{L}$ that reproduces the observed $a_{2}/a_{1}$ with resonance-locked orbits. However, this is not the case for systems of three planets, since we still have only one free parameter $\mathcal{L}$ (whose value is linked to the initial captured state, not constrained observationally) but two observables, that is the two pairs’ semi-major 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 semi-major 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 semi-major 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 semi-major 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 disk are of order $\sqrt{\tau_{e}/\tau_{a}}\sim h$, where $\tau_{a}$ and $\tau_{e}$ are the timescales of migration and eccentricity damping respectively, and $h=H/r\sim 0.05$ is the aspect ratio of the disk (Goldreich & Schlichting 2014, Pichierri et al. 2018). Since disks 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 disk 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 semi-major axis ratios (see also the Appendix A).

The remainder of this paper is organised as follows. In Section 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 three-planets systems undergoing dissipation. In Section 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 Section 4 and we finally conclude by discussing their significance in Section 5.

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 (2013), 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 well-known Hamiltonian of the restricted, circular three-body problem of a massless particle in resonance with a massive unperturbed body. In particular, such a Hamiltonian is integrable and is equivalent to the so-called 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 Section 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

[TABLE]

where the keplerian part is given by

[TABLE]

and describes the (integrable) motion of the three planets due to their interaction with the star, to which the small perturbation $\mathcal{H}_{pert}$ 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},\leavevmode\nobreak\ k^{out}>1$ are two positive integers. In other words, we assume the resonance conditions $n_{1}/n_{2}\simeq k^{in}/(k^{in}-1)$, $n_{2}/n_{3}\simeq k^{out}/(k^{out}-1)$, where for each planet $n=\sqrt{\mathcal{G}(M_{*}+m)a^{-3}}$ 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}\lambda_{2}-(k^{in}-1)\lambda_{1}-\varpi_{1}$, $k^{in}\lambda_{2}-(k^{in}-1)\lambda_{1}-\varpi_{2}$, $k^{out}\lambda_{3}-(k^{out}-1)\lambda_{2}-\varpi_{2}$, and $k^{out}\lambda_{3}-(k^{out}-1)\lambda_{2}-\varpi_{3}$ remain in the Hamiltonian, where $\lambda$ is the mean longitude of a planet, and $\varpi$ is its longitude of the periastron.

The resonant perturbing Hamiltonian expanded to first order in the eccentricities reads

[TABLE]

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 semi-major axis ratios, and their expressions may be found in Murray & Dermott (2000). We therefore write the Hamiltonian after the averaging procedure as

[TABLE]

and then drop the higher order terms. Note that terms that describe the mutual influence of the innermost and outermost planet are not included in $\mathcal{H}_{res}$ 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 action-angle variables $(\Lambda,\Gamma,\lambda,\gamma)$ (omitting the subscripts 1,2,3 for simplicity), which are given in terms of the orbital elements by

[TABLE]

where $\mu=M_{*}m/(M_{*}+m)$ is the reduced mass, and $\ell=E-e\sin E$ is the mean anomaly ($E$ being the eccentric anomaly). In these variables, the Keplerian part $\mathcal{H}_{kepl}$ of the Hamiltonian (2.1) takes the form

[TABLE]

while the resonant Hamiltonian writes

[TABLE]

note that in going from (2.3) to (2.7) we have made use of the approximation $e\simeq\sqrt{2\Gamma/\Lambda}$, 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

[TABLE]

it is straightforward to check using the Poisson bracket criterion (Morbidelli, 2002) that it is indeed canonical. Now, the new angle $\kappa$ does not appear explicitly in the Hamiltonian, which makes $\mathcal{K}$ a constant of motion. The significance of $\mathcal{K}$ has already been discussed for two planets (e.g. Michtchenko et al. 2008, Batygin & Morbidelli 2013), 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 semi-major axis ratios do not satisfy exactly the resonant condition $a_{i}/a_{i+1}=((k-1)/k)^{2/3}$, therefore the observed $a_{i,obs}$ do not alone reveal how close the planets are to resonance, nor do they represent the nominal $\bar{a}_{i}$ that do satisfy it. However by calculating from $a_{i,obs}$ the value of the constant of motion $\mathcal{K}$, and imposing in the formula

[TABLE]

the condition of exact resonance, $a_{i}/a_{i+1}=((k-1)/k)^{2/3}$, for all pairs $i=1,2$, we derive the nominal value of $\bar{\Lambda}_{3}$. From this, one easily obtains $\bar{a}_{3}$ from $\bar{a}_{3}=(\bar{\Lambda}_{3}/m_{3})^{2}/(\mathcal{G}M_{*})$, and then recursively $\bar{a}_{2}=((k^{out}-1)/k^{out})^{2/3}\bar{a}_{3}$, and finally $\bar{a}_{1}=((k^{in}-1)/k^{in})^{2/3}\bar{a}_{2}$.

It is worth briefly recalling here why even in resonance the planets’ semi-major 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}\lambda_{2}-(k^{in}-1)\lambda_{1}-\varpi_{1}$ is librating, $0\sim k^{in}\dot{\lambda}_{2}-(k^{in}-1)\dot{\lambda}_{1}-\dot{\varpi}_{1}=k^{in}n_{2}-(k^{in}-1)n_{1}-\dot{\varpi}_{1}$, which together with the condition of exact nominal resonance $k^{in}n_{2}-(k^{in}-1)n_{1}=0$ would imply $\dot{\varpi}_{1}\sim 0$; however from the Hamiltonian (2.7) we have $\dot{\varpi}_{1}=-\dot{\gamma}_{1}\propto\Gamma_{1}^{-1/2}\sim 1/e_{1}$ which grows as $e_{1}\searrow 0$, meaning that at low eccentricities $\dot{\varpi}_{1}\nsim 0$, which in turn forces $k^{in}\dot{\lambda}_{2}-(k^{in}-1)\dot{\lambda}_{1}=k^{in}n_{2}-(k^{in}-1)n_{1}\nsim 0$ in order to maintain the libration of the resonant angle. The resonant equilibrium points will therefore correspond to semi-major axes $a_{i}$ which may well deviate farther and farther from $\bar{a}_{i}$ as $e_{i}$ approaches 0. We can however already greatly simplify the calculations given that we will only consider deviations of the semi-major 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 $\delta\Lambda_{i}=\Lambda_{i}-\bar{\Lambda}_{i}$, where $\bar{\Lambda}_{i}=\mu_{1}\sqrt{\mathcal{G}(M_{*}+m_{i})\bar{a}_{1}}$ is the nominal resonant value of $\Lambda_{i}$, and write

[TABLE]

which, inserting the definition of $\delta\Lambda_{i}$ and dropping the unimportant constant term and the higher order terms, reduces to:

[TABLE]

where $\bar{n}_{i}=\sqrt{\mathcal{G}(M_{*}+m_{i})/\bar{a}^{3}}$ is the nominal mean motion and $\bar{h}_{i}=\bar{n}_{i}/\bar{\Lambda}_{i}=1/(\mu_{i}\bar{a}_{i}^{2})$ 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 $\mathcal{O}(\delta\Lambda^{2})$ does not introduce any significant inaccuracy in our calculations. Concerning the resonant Hamiltonian (2.7), we can evaluate it at the nominal values $\Lambda=\bar{\Lambda}$ as it is already of order $\mathcal{O}(e)$.

Finally, one last canonical change of variable is made:

[TABLE]

Again, we see that the new angle $\theta^{\prime}$ does not appear in the Hamiltonian, making $\Omega$ another constant of motion of the system (note that here $\Omega$ 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 four-degree-of-freedom Hamiltonian which depends parametrically on the constants of motion $\mathcal{K}$, $\Omega$. We already mentioned the meaning of $\mathcal{K}$; for $\Omega$, one can easily show that $\mathcal{K}+\Omega=(\Lambda_{1}-\Gamma_{1})+(\Lambda_{2}-\Gamma_{2})+(\Lambda_{3}-\Gamma_{3})\equiv\mathcal{L}$, 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 4-degrees-of-freedom Hamiltonian $\bar{\mathcal{H}}$ which is a function of the actions \big{(}\Psi_{1}^{(1)},\Psi_{1}^{(2)},\Psi_{2}^{(1)},\Psi_{2}^{(2)}\big{)} and the angles \big{(}\psi_{1}^{(1)},\psi_{1}^{(2)},\delta\gamma^{(1)},\delta\gamma^{(2)}\big{)}, and depends parametrically on the values of $\mathcal{K}$ and $\Omega$ (which are linked to the orbital elements as expressed in (2) and (2)); the Hamiltonian in these variables reads

[TABLE]

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

[TABLE]

Note that by the functional form of the Hamiltonian, any combination of values in $\{0,\pi\}$ 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,\pi\}$ reduces the problem of solving the four equations that appear in the first line to find the stable equilibria of the system. Note that although the Hamiltonian depends on both $\Omega$ and $\mathcal{K}$, 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 e.g. $\bar{a}_{1}=1$. To trace out the loci of the resonant equilibria, we then simply change the value of $\Omega$ (which corresponds to changing the angular momentum $\mathcal{L}$, at constant $\mathcal{K}$) and solve equations (2.14) to find \big{(}\Psi_{1,eq}^{(1)},\Psi_{1,eq}^{(2)},\Psi_{2,eq}^{(1)},\Psi_{2,eq}^{(2)}\big{)}, which are then translated into orbital elements working backwards through the canonical transformations.

We show in Figure 2 one example of equilibrium curves for three equal-mass 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 $N$-body 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).222 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 Section 4 of Batygin & Morbidelli (2013).

These $N$-body integrations with the addition of dissipative effects will be detailed below in Section 3.3.

2.2 Resonant repulsion for three-planets 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 three-planets 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 semi-major axes. That is, because the resonant condition requires $k\dot{\lambda}_{2}-(k-1)\dot{\lambda}_{1}-\dot{\varpi}_{1,2}\simeq 0$ and since the perihelion precession is always retrograde, $\dot{\varpi}_{1,2}<0$, one necessarily has $k\dot{\lambda}_{2}-(k-1)\dot{\lambda}_{1}<0$ i.e. $a_{2}/a_{1}>(k/(k-1))^{2/3}$. More concretely, at low enough eccentricities and at semi-major axis ratios close to the nominal ones, one finds directly using the resonant Hamiltonian expanded to first order in $e$ that $\dot{\varpi}_{1}=\alpha f_{res}^{(1)}n_{1}(m_{2}/M_{*})e_{1}^{-1}$, $\dot{\varpi}_{2}=-f_{res}^{(2)}n_{2}(m_{1}/M_{*})e_{2}^{-1}$ with $f_{res}^{(1)}<0$, $f_{res}^{(2)}>0$: this means that the lower are the eccentricities, the wider are the equilibria from the exact commensurability. At higher eccentricities the secular terms, of $\mathcal{O}(e^{2})$, become more important, and they contribute a positive contribution (that is constant in $e$) in $\dot{\varpi}$; however, one still finds $\dot{\varpi}_{1,2}<0$ 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 Figure 2, where the analytically computed resonant equilibria agree very well with our numerical simulations.333 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.

Consider now a resonant 3-planet 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 semi-major axes are expected to diverge. We also conclude that systems of three planets that are close to a given first-order 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 (2.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 e.g. in Batygin & Morbidelli (2013) for the 2-planets case, and only works when using an expanded Hamiltonian and for semi-major axes close to the nominal resonant ones, which are our working assumptions anyway. Therefore, after rescaling all planetary masses by a certain factor $\tilde{m}$, 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 semi-major axis ratios, the values of the eccentricities that correspond to the resonant equilibrium point are just rescaled by $\tilde{m}$, since the eccentricities and the planetary masses appear as a product in the perturbing Hamiltonian. This can be easily understood using the previous formula $\dot{\varpi}\propto(m_{pl}/M_{*})e^{-1}$, and noticing that fixing the semi-major axis ratio ultimately fixes $\dot{\varpi}$ by the resonance condition; therefore, rescaling the planetary masses, at fixed $\dot{\varpi_{1}}=\dot{\varpi_{2}}=\dot{\varpi_{3}}$ (i.e. at constant semi-major axis ratios), the eccentricities are simply rescaled by the same factor. This implies that for a given equilibrium configuration of the semi-major axis ratio $a_{2,eq}/a_{1,eq}$, the corresponding equilibrium of the ratio $a_{3,eq}/a_{2,eq}$ will be independent of $\tilde{m}$, 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 the Appendix A for an explicit presentation of this rescaling procedure).

3 A scenario for dissipative evolution of 3-planet systems

In this section, we select near-resonant 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 disk 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 disk which causes the planets’ semi-major 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 disk), 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 disk, 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 $\sim 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 Figure 2, however, this configuration with semi-major axis ratio wide of resonance and low eccentricity is very close to the resonant equilibrium points. The planets keep approaching each 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 disk 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 $\Omega$, this time decreasing their eccentricities and hence increasing the semi major axis ratios $a_{i+1}/a_{i}$ for each planet-planet pair. Note that $\Omega$ changes because $\mathcal{K}$ changes (since do the semi-major axes as a consequence of the dissipation of energy) and $\mathcal{L}$ 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 semi-major 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 short-period planets this is equivalent to obtaining the proper value of the periods, so that we can directly compare the corresponding observed semi-major axis ratios with the ones coming from our averaged model of resonance.444 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. 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_{\oplus}(R_{pl}/R_{\oplus})$, where $m_{pl}$, $R_{pl}$ are the mass and radius of the exoplanet, and $R_{\oplus}$, $m_{\oplus}$ 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 Section 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 $N$-body 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}\simeq(k-1)/k$, one can define (see for example Lithwick et al. 2012; Hadden & Lithwick 2014)

[TABLE]

called the normalised distance from (exact) resonance. Note that when $\Delta>0$ the planets reside wide of the $k:k-1$ resonance, while when $\Delta<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 $|\Delta|\lesssim 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 $\sim(m/M_{*})^{2/3}$ (Deck et al. 2013; Batygin 2015), where the average planetary mass for Kepler systems is of order $m\sim 3\times 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 $\Delta\sim 10^{-3}$ in normalised units, meaning that in selecting systems with $0<\Delta\lesssim 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 $\Delta$ favours values between [math] and $\lesssim 0.05$ (Hadden & Lithwick (2017), see also Figure 3). Additionally, we will require that $\Delta>0$, which is justified by our results in the previous section.

Of all 3-planet systems, only 8 satisfy $|\Delta|<0.05$ for both pairs, that is, appear to be close to a multi-resonant chain (they are Kepler-31, Kepler-53, Kepler-114, Kepler-207, Kepler-289, Kepler-305, Kepler-326, YZ Ceti). The architecture of these systems is shown in Figure 4; of these, only 3 satisfy $\Delta>0$ for both pairs. These are Kepler-31, Kepler-305 and YZ Ceti. For these systems, we consider whether or not their observed orbital configuration can be consistent with the scenario envisioned above.

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 Section 2.1 (in practice, we work with the aforementioned Hamiltonian rescaled by the common planetary mass factor $\tilde{m}$, see the Appendix A). Then, we find the resonant equilibrium point (i.e. the value of $\mathcal{L}$) that corresponds to a value of $a_{3}/a_{1}$ equal to the observed semi-major axis ratio, and therefore put $\left.(a_{3}/a_{1})\right|_{eq}\equiv\left.(a_{3}/a_{1})\right|_{obs}$. The determined value of $\mathcal{L}$ 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 semi-major 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 $\mathcal{L}$) 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 $\mathcal{L}$ which gives a resonant equilibrium configuration such that the semi-major axis ratio is equal to the observed one (provided that the latter is wide of the nominal value $(k/(k-1))^{2/3}$). In the case of three planets, choosing $\mathcal{L}$ such that $\left.(a_{3}/a_{1})\right|_{eq}\equiv\left.(a_{3}/a_{1})\right|_{obs}$ automatically fixes the equilibrium values of both semi-major axis ratios $a_{2}/a_{1}$ and $a_{3}/a_{2}$. Considering for example the corresponding equilibrium $\left.(a_{3}/a_{2})\right|_{eq}$, we obtain the weighted difference

[TABLE]

between $\left.(a_{3}/a_{2})\right|_{eq}$ and the observed value $\left.(a_{3}/a_{2})\right|_{obs}$. The same can be done for $a_{2}/a_{1}$, which gives a similar (absolute) result, $|\bar{\delta}(a_{2}/a_{1})|\simeq|\bar{\delta}(a_{3}/a_{2})|$. Maintaining this procedure, we loop over different planetary mass ratios.

We applied this procedure to the three systems selected above, starting with Kepler-305, 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 Figure 5 equilibrium curves (equivalent to those shown in Figure 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 semi-major 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 $\left.(a_{3}/a_{1})\right|_{eq}\equiv\left.(a_{3}/a_{1})\right|_{obs}$; finally, we obtain $\bar{\delta}(a_{3}/a_{2})$ as defined in (3.2) (and similarly for $\bar{\delta}(a_{2}/a_{1})$). As explained above, imposing $\left.(a_{3}/a_{1})\right|_{eq}\equiv\left.(a_{3}/a_{1})\right|_{obs}$ 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 $\max e_{eq}$ rescaled by a common planetary mass factor $\tilde{m}$ in order to obtain results that are independent of the planet-to-star mass ratio, since again the latter quantity does not effect equilibrium semi-major axis ratio configurations, and therefore does not affect $\bar{\delta}(a_{3}/a_{2})$. Panels (a) and (b) in Figure 6 are the result of this procedure spanning different planetary mass ratios, showing maps of $\bar{\delta}(a_{2}/a_{1})$ and $\max e/\tilde{m}$ using the observed semi-major axis ratios of Kepler-305 (the bottom panels (c) and (d) show the results of numerical simulations which will be detailed in Section 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 Figure 7, and for For Kepler-31, a chain (close to the) 2:1 – 2:1 mean motion resonances, in Figure 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 first-order 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 disks by adding fictitious forces which mimic the interaction with the disk (Cresswell & Nelson 2006, 2008) to the $N$-body 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 so-called planetary trap, which reproduces the effect of the inner edge of the disk and describes a disk 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 $\sim 0.2$ and $\sim 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 semi-major axis decreases, which is practically inconvenient. Therefore, we used a fictitious $\tau_{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 disk-migration phase the semi-major axis ratios are smaller than the observed ones, since the subsequent evolution dominated by tidal dissipation will only cause the semi-major axis ratios to expand. 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 $\tau_{e}$ and the migration rate $\tau_{a}$, and thus obtain different corresponding equilibrium values of the semi-major axis ratios; the latter approach the exact commensurabilities as the eccentricities grow (Figure 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. 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 $\mathcal{L}$, so that these results are in fact generalisable to any dissipative process such that $\dot{E}<0$ and $\dot{\mathcal{L}}=0$. Therefore, a resonant system undergoing any such process is expected to follow the loci of the resonant equilibria, and the divergence of the semi-major 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 semi-major 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 $\left.(a_{3}/a_{1})\right|_{obs}$, 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 $\left.(a_{3}/a_{2})\right|_{obs}$ (since they are obtained from the mean period extracted from the data) by computing the relative difference as in (3.2). 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 Section 3.2), as well as the quantity $\max e/\tilde{m}$.

This analysis has been performed for the three selected systems. For Kepler-305 and YZ Cet, we show the resulting plots on the bottom panels (c) and (d) of Figures 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 $\bar{\delta}(a_{3}/a_{2})$ is due to the fact that the numerically simulated systems are undergoing fast oscillation while they cross the observed value $\left.(a_{3}/a_{1})\right|_{obs}$, but the typical value of $\bar{\delta}(a_{3}/a_{2})$ is similar to the one found analytically.

The case of Kepler-31, 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 disk-planet 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 planet-disk interactions represent only approximate implementations of the real forces that are felt by the planets from the disk, 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}\gtrsim 1$) the numerical results agree very well with the analytical ones; moreover we still observe that the theoretical value of $\bar{\delta}(a_{3}/a_{2})$ varies extremely little with $m_{1}/m_{2}$ and $m_{2}/m_{3}$ (Figure 8), so the latter simulations can be considered as enough support for the analytical calculations.

4 Results

4.1 Probabilistic measure of a resonant configuration in Kepler-305, YZ Cet and Kepler-31

Using the maps of $\bar{\delta}(a_{3}/a_{2})$ shown in Figures 6, 7 and 8 for the three selected systems Kepler-305, YZ Cet and Kepler-31, we draw the following conclusions. First of all, one might expect that the quantity $\bar{\delta}(a_{3}/a_{2})$ 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. Follow-up 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 $\bar{\delta}(a_{3}/a_{2})$ is small, but never vanishing, which would represent an analytically computed equilibrium configuration such that $\bar{\delta}(a_{3}/a_{2})=0$, that is, a resonant equilibrium point which satisfies $\left.(a_{3}/a_{2})\right|_{eq}=\left.(a_{3}/a_{2})\right|_{obs}$ and $\left.(a_{2}/a_{1})\right|_{eq}=\left.(a_{2}/a_{1})\right|_{obs}$. But even if this happened to be the case, the level curve $\bar{\delta}(a_{3}/a_{2})=0$ 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 $\bar{\delta}(a_{3}/a_{2})$, meaning that these systems are indeed close to some equilibrium point of the Hamiltonian (2.13) and therefore could potentially reside in a multi-resonant chain. However, these values by themselves do not contain any meaningful information. The quantity $\bar{\delta}(a_{3}/a_{2})$ 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<\Delta<0.05$ for each pair, and extract the corresponding $\bar{\delta}(a_{3}/a_{2})$ (calculated for the choice of mass ratios $m_{1}/m_{2}=m_{2}/m_{3}=1$ for simplicity, since as we saw above the $\bar{\delta}(a_{3}/a_{2})$ value depends extremely weakly on the mass ratios). From the cumulative distribution of $|\bar{\delta}(a_{3}/a_{2})|$ that arises from this procedure we can obtain the probability that any given system has a given (small) $\bar{\delta}(a_{3}/a_{2})$ purely by chance.

Since we are interested mainly in the 2:1 and 3:2 mean motion resonances, we show in Figure 9 these cumulative distributions for systems close to any possible combinations of these resonances. The results show that the proximity of YZ Cet to the 3:2 – 3:2 resonance is not statistically significant, since in $\sim 80\%$ of randomly generate systems close to the 3:2 – 3:2 chain we obtain an equivalent or smaller value of $\bar{\delta}(a_{3}/a_{2})$. Instead, Kepler-305 and Kepler-31 are likely to be in resonance at the $1\sigma$ level (i.e. the probability that their value of $\bar{\delta}(a_{3}/a_{2})$ is smaller than the determined value by chance is less than $32\%$): for Kepler-305 there is a $\sim 15\%$ chance that this particular system lies this close to resonance by chance, while for Kepler-31 the probability is $\sim 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<\Delta<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 Section 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 $|\Delta|<\max\Delta=0.05$ was dictated by the observation of the location of the peak wide of nominal resonance (Hadden & Lithwick (2017) and our Figure 3), so restricting the interval in $\Delta$ values with smaller $\max\Delta$ might have resulted in excluding potential systems. On the other hand, increasing $\max\Delta$ 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\Delta=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 three-body Laplace angle $\dot{\varphi_{L}}$ in these systems, since its libration can be in principle a sufficient condition to conclude that they are indeed resonant. For Kepler-305 we checked that the Laplace angle $\varphi_{L}=2\lambda_{1}-4\lambda_{2}+2\lambda_{3}$ satisfies $\dot{\varphi}_{L}\simeq 0.5\degree{\leavevmode\nobreak\ \mathrm{days}}^{-1}$ given the observed transits periods; for Kepler-31, $\varphi_{L}=\lambda_{1}-3\lambda_{2}+2\lambda_{3}$ satisfies $\dot{\varphi}_{L}\simeq 0.1\degree{\leavevmode\nobreak\ \mathrm{days}}^{-1}$; finally for YZ Cet, $\varphi_{L}=2\lambda_{1}-5\lambda_{2}+3\lambda_{3}$ satisfies $\dot{\varphi}_{L}\simeq 9.4\degree{\leavevmode\nobreak\ \mathrm{days}}^{-1}$. 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 $\bar{\delta}(a_{3}/a_{2})$ 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.

4.2 The 5:4 – 4:3 resonant chain on Kepler-60 and other near-resonant 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 Kepler-60 represents a notable example. This system hosts three planets with mean observed periods of $\simeq 5.49{\leavevmode\nobreak\ \mathrm{days}}$, $\simeq 8.29{\leavevmode\nobreak\ \mathrm{days}}$ and $\simeq 16.74{\leavevmode\nobreak\ \mathrm{days}}$ respectively. Their masses have been constrained via TTV to be $\simeq 4M_{\oplus}$ for all planets (Jontof-Hutter et al., 2016). The mean motions of the planets satisfy $4\dot{\lambda}_{1}-8\dot{\lambda}_{2}+4\dot{\lambda}_{3}\simeq-0.02\degree{\leavevmode\nobreak\ \mathrm{days}}^{-1}$, hinting at a resonant configuration. Indeed, Goździewski et al. (2016) found that the TTV signal for these planets is consistent with a true three-body Laplace-like resonance as well as a chain of 5:4 – 4:3 two-body mean motion resonances. Using the system’s parameters we can find a resonant equilibrium configuration as in Section 3.2, by imposing $a_{3}/a_{1}$ to be equal to the observed $\left.(a_{3}/a_{1})\right|_{obs}$. This gives $\bar{\delta}(a_{3}/a_{2})$ of order $4\times 10^{-5}$. Using an analogous argument to that of Figure 9, we find that there is only a 0.25% probability that Kepler-60 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}\simeq 0.02$, $e_{2}\simeq 0.03$, $e_{3}\simeq 0.01$ for the observed planetary masses. These numbers are quite close to the ones consistent for the two-body mean motion resonance chain solution found in Goździewski et al. (2016). Note in passing that their solution is for non-vanishing 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 near-resonant systems of three planets with $k>3$ that are found in the catalogue. The only ones which satisfy our criterion $|\Delta|<0.05$ for both pairs are K2-239 (close to a 3:2 – 4:3 chain), Kepler-289 (close to a 2:1 – 2:1 chain), Kepler-226 (close to a 4:3 – 3:2 chain) and Kepler-431 (close to a 5:4 – 4:3). Of these, only the latter two satisfy $\Delta>0$ for both pairs.

5 Conclusions

In this work, we have generalised the formalism of dissipative divergence of resonant orbits to multi-resonant chains. The analytical study performed in Section 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 Section 2.2 that under the effect of slow dissipation a nearly-resonant system is expected to follow the loci of the equilibrium points of the resonant Hamiltonian (2.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_{\mathrm{pl}}\gtrsim 1.6R_{\oplus}$ 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 disk was still present. The associated planet-disk interaction would then force the planets’ orbital elements to change, i.e. 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 disk or $e$-dependent 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 disk-planet interaction, but the subsequent evolution after the disk 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 Trappist-1, Kepler-223, GJ867 and Kepler-60, 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 $\bar{\delta}(a_{3}/a_{2})$ defined in (3.2), which is obtained directly from the observed semi-major 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 Figure 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 $\sim 15\%$, $\sim 20\%$ and $\sim 80\%$ probability that Kepler-305, Kepler-31 and YZ Cet respectively find themselves this close to resonance merely by chance. Multiplying these probabilities we find that there is only a $\sim 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<\Delta<0.05$ show memory of the processes of resonant capture; this is consistent with the histogram of Figure 3, where the peak wide of the resonance is about 2 times higher than the underlying random-like 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 follow-up via further ground-based and space-based 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 present-day 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 planet-to-star mass factor $\tilde{m}$. In this appendix, we detail the construction of this Hamiltonian and its use in the paper.

Consider three planets, whose physical parameters are labeled 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 mean-motion resonance, with nominal semi-major axes $\bar{a}$ and that the deviations of the semi-major 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 Section 2. Having fixed the planet-planet mass ratios $m_{1}/m_{2}=\beta_{1}$ and $m_{2}/m_{3}=\beta_{2}$, we introduce the average planet-star mass ratio

[TABLE]

Inverting this expression we easily get all planetary masses in terms of $\tilde{m}$,

[TABLE]

with coefficients $c$ depending on $M_{*}$, $\beta_{1}$ and $\beta_{2}$ only.

We introduce the modified Delaunay action-angle variables $(\Lambda,\Gamma,\lambda,\gamma)$ as in (2), but we rescale the actions by the common mass factor $\tilde{m}$: this gives the following definition for the $\Lambda$’s (we maintain the same notation as the non-rescaled actions for simplicity)

[TABLE]

and the same formal definition of $\Gamma=\Lambda e^{2}/2$ at lowest order in $e$. We now introduce the Hamiltonian for the problem, that is the sum of the Keplerian Hamiltonian (2.6) and the resonant interaction Hamiltonian (2.7) to first order in $e$. However since we have rescaled the actions by $\tilde{m}$, in order for the Hamilton equations to be conserved we must also rescale the Hamiltonian itself by $\tilde{m}$. As in Section 2, since we are not considering large deviations in the semi-major axes from their nominal values, and since the resonant Hamiltonian is already of order $\mathcal{O}(e)$, we evaluate the resonant Hamiltonian on the nominal values $\bar{\Lambda}$ defined from $\bar{a}$ using (A.3). It is then easy to see that the rescaled Keplerian Hamiltonian takes the form

[TABLE]

and is therefore independent of $\tilde{m}$, while the rescaled resonant part will have a multiplicative coefficient $\tilde{m}$:

[TABLE]

From here, the sequence of changes of variables detailed in Section 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 $\Lambda$ and $\Gamma\propto e^{2}$ one can see the following. Assuming a fixed equilibrium value of the semi-major axes (that is, of the $\Lambda$’s), a change in the planet-star mass factor $\tilde{m}$ will have the only effect to rescale the equilibrium values of all $\sqrt{\Gamma}\propto e$ by the same quantity. In this configuration, the equilibria of the semi-major axis ratios $(a_{2}/a_{1})_{eq}$ and $(a_{3}/a_{2})_{eq}$ remain independent of $\tilde{m}$.