© G. Pichierri et al. 2019

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 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 disc-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 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 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, 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 $a_1/a_2={((k-1)/k)}^{2/3}$ (where a₁ and a₂ 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 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 K-dwarf HD 40307, which hosts¹ 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 (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₁∕a₁ − m₂∕a₂, 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 (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 multi-systems 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 multi-planets systems mentioned above.

In this paper we focus on slow convergent Type-I 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 post-disc 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 $\stackrel{.}{\mathcal{L}}=0$. We conclude that our specific implementation of Type-I 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 near-resonant period ratios that emerges from available data.

In order to answer this question, we examined the NASA Exoplanet Archive² and selected confirmed three-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⁽ⁱⁿ⁾:(k⁽ⁱⁿ⁾ − 1) and k^(out):(k^(out) − 1) resonances respectively (so that the resonant angles k⁽ⁱⁿ⁾λ₂ − (k⁽ⁱⁿ⁾ − 1)λ₁ − ϖ₂ and k^(out)λ₃ − (k^(out) − 1)λ₂ − ϖ₂ are librating), then the Laplance angle φ_L = (k⁽ⁱⁿ⁾ − 1)λ₁ − (k⁽ⁱⁿ⁾ + k^(out) − 1)λ₂ + k^(out)λ₃ 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 ϖ, 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 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 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 $\mathcal{L}$ a free parameter: it is therefore always possible to find a value of $\mathcal{L}$ that reproduces the observed a₂∕a₁ 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₂∕a₁ and a₃∕a₂, 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 disc are of order $\sqrt{{\tau }_e/{\tau }_a}~h$, 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₁ ∕m₂ and m₂ ∕m₃, 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 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 three-planets 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 planet-hosting stars reveal that multi-planetary 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 multi-planetary 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 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 Sect. 3.2.

Consider three planets of masses m₁, m₂ and m₃, 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 $$\mathcal{H}={\mathcal{H}}_{\text{kepl}}+{\mathcal{H}}_{\text{pert}},$$(1)

where the keplerian part is given by $${\mathcal{H}}_{\text{kepl}}=-\frac{\mathcal{G}{M}_{*}{m}_1}{2a_1}-\frac{\mathcal{G}{M}_{*}{m}_2}{2a_2}-\frac{\mathcal{G}{M}_{*}{m}_3}{2a_3},$$(2)

and describes the (integrable) motion of the three planets due to their interaction with the star, to which the small perturbation ${\mathcal{H}}_{\text{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ⁱⁿ:(kⁱⁿ − 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ⁱⁿ,  k^out > 1 are two positive integers. In other words, we assume the resonance conditions n₁ ∕n₂ ≃ kⁱⁿ∕(kⁱⁿ − 1), n₂ ∕n₃ ≃ 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ⁱⁿ λ₂ − (kⁱⁿ − 1)λ₁ − ϖ₁, kⁱⁿ λ₂ − (kⁱⁿ − 1)λ₁ − ϖ₂, k^out λ₃ − (k^out − 1)λ₂ − ϖ₂, and k^out λ₃ − (k^out − 1)λ₂ − ϖ₃ 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 $$\begin{array}{ll}{\mathcal{H}}_{\text{res}}=\hfill & -\frac{\mathcal{G}{m}_1{m}_2}{a_2}(f_{\text{res}}^{(1,\text{in})}{e}_1\cos (k^{\text{in}}{\lambda }_2-(k^{\text{in}}-1){\lambda }_1-{\varpi }_1)\hfill \\ \hfill & +f_{\text{res}}^{(2,\text{in})}{e}_2\cos (k^{\text{in}}{\lambda }_2-(k^{\text{in}}-1){\lambda }_1-{\varpi }_2))\hfill \\ \hfill & -\frac{\mathcal{G}{m}_2{m}_3}{a_3}(f_{\text{res}}^{(1,\text{out})}{e}_2\cos (k^{\text{out}}{\lambda }_3-(k^{\text{out}}-1){\lambda }_2-{\varpi }_2)\hfill \\ \hfill & +f_{\text{res}}^{(2,\text{out})}{e}_3\cos (k^{\text{out}}{\lambda }_3-(k^{\text{out}}-1){\lambda }_2-{\varpi }_3)),\hfill \end{array}$$(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 semi-major axis ratios, and their expressions may be found in Murray & Dermott (2000). We therefore write the Hamiltonian after the averaging procedure as $$\overline{\mathcal{H}}={\mathcal{H}}_{\text{kepl}}+{\mathcal{H}}_{\text{res}}+\mathcal{O}(e^2,I^2),$$(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 ${\mathcal{H}}_{\text{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 (Λ, Γ, λ, γ) (omitting the subscripts 1,2,3 for simplicity), which are given in terms of the orbital elements by $$\begin{array}{lll}\Lambda \hfill & =\mu \sqrt{\mathcal{G}(M_{*}+m)a}\simeq m\sqrt{\mathcal{G}{M}_{*}a},\hfill & \lambda =\mathcal{l}+\varpi ,\hfill \\ \Gamma \hfill & =\Lambda (1-\sqrt{1-e^2})\simeq \Lambda e^2/2,\hfill & \gamma =-\varpi ,\hfill \end{array}$$(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 ${\mathcal{H}}_{\text{kepl}}$ of the Hamiltonian (1) takes the form $${\mathcal{H}}_{\text{kepl}}=-{\displaystyle \sum _{i=1}^3}\frac{{\mathcal{G}}^2{(M_{*}+m_i)}^2{\mu }_i^3}{2{\Lambda }_i^2}\simeq -{\displaystyle \sum _{i=1}^3}\frac{m_i^3}{2}{\left(\frac{\mathcal{G}{M}_{*}}{{\Lambda }_i}\right)}^2,$$(6)

while the resonant Hamiltonian writes $$\begin{array}{lll}{\mathcal{H}}_{\text{res}}\hfill & =\hfill & -\frac{{\mathcal{G}}^2{M}_{*}{m}_1{m}_2^3}{{\Lambda }_2^2}\hfill \\ \hfill & \hfill & \times (f_{\text{res}}^{(1,\text{in})}\sqrt{\frac{2{\Gamma }_1}{{\Lambda }_1}}\cos (k^{\text{in}}{\lambda }_2-(k^{\text{in}}-1){\lambda }_1+{\gamma }_1)\hfill \\ \hfill & \hfill & +f_{res}^{(2,in)}\sqrt{\frac{2{\Gamma }_2}{{\Lambda }_2}}\cos (k^{\text{in}}{\lambda }_2-(k^{\text{in}}-1){\lambda }_1+{\gamma }_2))\hfill \\ \hfill & \hfill & -\frac{{\mathcal{G}}^2{M}_{*}{m}_2{m}_3^3}{{\Lambda }_3^2}\hfill \\ \hfill & \hfill & \times (f_{\text{res}}^{(1,\text{out})}\sqrt{\frac{2{\Gamma }_2}{{\Lambda }_2}}\cos (k^{\text{out}}{\lambda }_3-(k^{\text{out}}-1){\lambda }_2+{\gamma }_2)\hfill \\ \hfill & \hfill & +f_{\text{res}}^{(2,\text{out})}\sqrt{\frac{2{\Gamma }_3}{{\Lambda }_3}}\cos (k^{\text{out}}{\lambda }_3-(k^{\text{out}}-1){\lambda }_2+{\gamma }_3)),\hfill \end{array}$$(7)

we note that in going from Eqs. (3) to (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 $$\begin{array}{lll}\hfill & \hfill & \mathcal{K}={\Lambda }_1+\frac{k^{\text{in}}-1}{k^{\text{in}}}{\Lambda }_2+\frac{(k^{\text{in}}-1)(k^{\text{out}}-1)}{k^{\text{in}}{k}^{\text{out}}}{\Lambda }_3,\kappa ={\lambda }_1,\hfill \\ \hfill & \hfill & {\Theta }^{(1)}=\frac{1}{k^{\text{in}}}{\Lambda }_2+\frac{k^{\text{out}}-1}{k^{\text{in}}{k}^{\text{out}}}{\Lambda }_3,{\theta }^{(1)}=k^{\text{in}}{\lambda }_2-(k^{\text{in}}-1){\lambda }_1,\hfill \\ \hfill & \hfill & {\Theta }^{(2)}=\frac{1}{k^{\text{out}}}{\Lambda }_3,{\theta }^{(2)}=k^{\text{out}}{\lambda }_3-(k^{\text{out}}-1){\lambda }_2;\hfill \end{array}$$(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 $\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 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 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 ā_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 $$\begin{array}{lll}\frac{\mathcal{K}}{{\Lambda }_3}\hfill & =\hfill & \frac{{\mu }_1}{{\mu }_3}\sqrt{\frac{M_{*}+m_1}{M_{*}+m_3}\frac{a_1}{a_3}}+\frac{k^{\text{in}}-1}{k^{\text{in}}}\frac{{\mu }_2}{{\mu }_3}\sqrt{\frac{M_{*}+m_2}{M_{*}+m_3}\frac{a_2}{a_3}}\hfill \\ \hfill & \hfill & +\frac{(k^{\text{in}}-1)(k^{\text{out}}-1)}{k^{\text{in}}{k}^{\text{out}}},\hfill \end{array}$$(9)

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 ${\overline{\Lambda }}_3$. From this, one easily obtains ā₃ from ${\overline{a}}_3={({\overline{\Lambda }}_3/m_3)}^2/(\mathcal{G}{M}_{*})$, and then recursively ${\overline{a}}_2={((k^{\text{out}}-1)/k^{\text{out}})}^{2/3}{\overline{a}}_3$, and finally ${\overline{a}}_1={((k^{\text{in}}-1)/k^{\text{in}})}^{2/3}{\overline{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ⁱⁿ λ₂ − (kⁱⁿ − 1)λ₁ − ϖ₁ is librating, $0~k^{\text{in}}{\dot{\lambda }}_2-(k^{\text{in}}-1){\dot{\lambda }}_1-{\dot{\varpi }}_1=k^{\text{in}}{n}_2-(k^{\text{in}}-1)n_1-{\dot{\varpi }}_1$, which together with the condition of exact nominal resonance kⁱⁿn₂ − (kⁱⁿ − 1)n₁ = 0 would imply ${\dot{\varpi }}_1~0$; however from the Hamiltonian (Eq. (7)) we have ${\dot{\varpi }}_1=-{\dot{\gamma }}_1\propto {\Gamma }_1^{-1/2}~1/e_1$ which grows as e₁ ↘ 0, meaning that at low eccentricities ${\dot{\varpi }}_1\nsim 0$, which in turn forces $k^{\text{in}}{\dot{\lambda }}_2-(k^{\text{in}}-1){\dot{\lambda }}_1=k^{\text{in}}{n}_2-(k^{\text{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 ā_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-{\overline{\Lambda }}_i$, where ${\overline{\Lambda }}_i={\mu }_1\sqrt{\mathcal{G}(M_{*}+m_i){\overline{a}}_1}$ is the nominal resonant value of Λ_i, and write $$\begin{array}{lll}{\mathcal{H}}_{\text{kepl}}\hfill & =\hfill & -{\displaystyle \sum _{i=1}^3}\frac{{\mathcal{G}}^2{(M_{*}+m_i)}^2{\mu }_i^3}{2}\hfill \\ \hfill & \hfill & \times \left(\frac{1}{{\overline{\Lambda }}_i^2}-2\frac{1}{{\overline{\Lambda }}_i^3}\delta {\Lambda }_i+3\frac{1}{{\overline{\Lambda }}_i^4}\delta {\Lambda }_i^2+\mathcal{O}(\delta {\Lambda }_i^3)\right),\hfill \end{array}$$(10)

which, inserting the definition of δΛ_i and dropping the unimportant constant term and the higher order terms, reduces to: $${\mathcal{H}}_{\text{kepl}}={\displaystyle \sum _{i=1}^3}\left(4{\overline{n}}_i{\Lambda }_i-\frac{3}{2}{\overline{h}}_1{\Lambda }_i^2\right),$$(11)

where ${\overline{n}}_i=\sqrt{\mathcal{G}(M_{*}+m_i)/{\overline{a}}^3}$ is the nominal mean motion and ${\overline{h}}_i={\overline{n}}_i/{\overline{\Lambda }}_i=1/({\mu }_i{\overline{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 (Eq. (7)), we can evaluate it at the nominal values $\Lambda =\overline{\Lambda }$ as it is already of order $\mathcal{O}(e)$.

Finally, one last canonical change of variable is made: $$\begin{array}{llll}\hfill & \Omega ={\Theta }^{(1)}+{\Theta }^{(2)}-({\Gamma }_1+{\Gamma }_2+{\Gamma }_3),\hfill & \hfill & {\theta }^{\prime }={\theta }^{(1)},\hfill \\ \hfill & {\Psi }_1^{(1)}={\Gamma }_1+{\Gamma }_2+{\Gamma }_3-{\Theta }^{(2)},\hfill & \hfill & {\psi }_1^{(1)}={\theta }^{(1)}+{\gamma }_1,\hfill \\ \hfill & {\Psi }_1^{(2)}={\Theta }^{(2)},\hfill & \hfill & {\psi }_1^{(2)}={\theta }^{(2)}+{\gamma }_2,\hfill \\ \hfill & {\Psi }_2^{(1)}=-{\Gamma }_2-{\Gamma }_3+{\Theta }^{(2)},\hfill & \hfill & \delta {\gamma }^{(1)}={\gamma }_1-{\gamma }_2,\hfill \\ \hfill & {\Psi }_2^{(2)}=-{\Gamma }_3,\hfill & \hfill & \delta {\gamma }^{(2)}={\gamma }_2-{\gamma }_3.\hfill \end{array}$$(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 four-degree-of-freedom Hamiltonian which depends parametrically on the constants of motion $\mathcal{K}$, Ω. We already mentioned the meaningof $\mathcal{K}$; for Ω, 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 $\overline{\mathcal{H}}$ which is a function of the actions $({\Psi }_1^{(1)},{\Psi }_1^{(2)},{\Psi }_2^{(1)},{\Psi }_2^{(2)})$ and the angles $({\psi }_1^{(1)},{\psi }_1^{(2)},\delta {\gamma }^{(1)},\delta {\gamma }^{(2)})$, and depends parametrically on the values of $\mathcal{K}$ and Ω (which are linked to the orbital elements as expressed in Eqs. (5) and (12)); the Hamiltonian in these variables reads $$\begin{array}{ll}\hfill & \overline{\mathcal{H}}={\mathcal{H}}_{\text{kepl}}+{\mathcal{H}}_{\text{res}},\hfill \\ \hfill & {\mathcal{H}}_{\text{kepl}}={\mathcal{H}}_{\text{kepl}}\left({\Psi }_1^{(1)},{\Psi }_1^{(2)};\mathcal{K},\Omega \right),\hfill \\ \hfill & {\mathcal{H}}_{\text{res}}={\mathcal{H}}_{\text{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),\hfill \end{array}$$(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 $$\begin{array}{lll}\hfill & \hfill & \frac{\partial \overline{\mathcal{H}}}{\partial {\Psi }_1^{(1)}}=0, \frac{\partial \overline{\mathcal{H}}}{\partial {\Psi }_1^{(2)}}=0, \frac{\partial \overline{\mathcal{H}}}{\partial {\Psi }_2^{(1)}}=0, \frac{\partial \overline{\mathcal{H}}}{\partial {\Psi }_2^{(2)}}=0,\hfill \\ \hfill & \hfill & \frac{\partial \overline{\mathcal{H}}}{\partial {\psi }_1^{(1)}}=0, \frac{\partial \overline{\mathcal{H}}}{\partial {\psi }_1^{(2)}}=0, \frac{\partial \overline{\mathcal{H}}}{\delta {\gamma }^{(1)}}=0, \frac{\partial \overline{\mathcal{H}}}{\delta {\gamma }^{(2)}}=0.\hfill \end{array}$$(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 $\mathcal{K}$, the latter assumes a natural value for any specific problem at hand (that is, any values of m₁, m₂, m₃ and of kⁱⁿ, k^out) by rescaling the units so that for example ā₁ = 1. To trace out the loci of the resonant equilibria, we then simply change the value of Ω (which corresponds to changing the angular momentum $\mathcal{L}$, at constant $\mathcal{K}$) and solve Eq. (14) to find $({\Psi }_{1,\text{eq}}^{(1)},{\Psi }_{1,\text{eq}}^{(2)},{\Psi }_{2,\text{eq}}^{(1)},{\Psi }_{2,\text{eq}}^{(2)})$, 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 equal-mass planets down to eccentricities of order 10⁻³, 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)³. These N-body 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 m1 = m2 = m3 = 10−4M*. 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 semi-major axis ratio. The location of the nominal resonant semi-major axis ratio (3/2)2/3 is shown bya vertical orange line. We also superimpose the numerically computed evolution of a three-planet 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 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$, that is $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 tofirst order in e that ${\dot{\varpi }}_1=\alpha f_{\text{res}}^{(1)}{n}_1(m_2/M_{*})e_1^{-1}$, ${\dot{\varpi }}_2=-f_{\text{res}}^{(2)}{n}_2(m_1/M_{*})e_2^{-1}$ with $f_{\text{res}}^{(1)}<0$, $f_{\text{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 Fig. 2, where the analytically computed resonant equilibria agree very well with our numerical simulations⁴.

Consider now a resonant three-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 (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 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 $\stackrel{.}{{\varpi }_1}=\stackrel{.}{{\varpi }_2}=\stackrel{.}{{\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₁ ∕m₂ and m₂ ∕m₃ 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 three-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 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’ 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 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 semi-major 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 $\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⁵. 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₁ + m₂ + m₃ by using for each system the mean planetary radius (R₁ + R₂ + R₃)∕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₁ and P₂ respectively, which are close to a given first order resonance, P₁∕P₂ ≃ (k − 1)∕k, one can define (see for example Lithwick et al. 2012; Hadden & Lithwick 2014) $${\Delta }_{k/(k-1)}:=\frac{k-1}{k}\frac{P_2}{P_1}-1,$$(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 $~{(m/M_{*})}^{2/3}$ (Deck et al. 2013; Batygin 2015), where the average planetary mass for Kepler systems is of order m ~ 3 × 10⁻⁵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⁻³ 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 three-planet systems, only 8 satisfy |Δ| < 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 Fig. 4; of these, only 3 satisfy Δ > 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.

Fig. 3 Distribution of the (signed) normalised distance from first order mean motion resonances $\Delta =\frac{k-1}{k}\frac{P_2}{P_1}-1$ with k =2, 3 in all exoplanetary systems (yellow) and for three-planet 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 three-planets systems only, 48 pairs have 0 <Δ ≲ 0.05, out of the 121 shown in the purple histogram.

Fig. 4 Orrery of the three-planets 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⁽ⁱⁿ⁾ and k^(out). We then pick both mass ratios m₂∕m₁ and m₃∕m₂ 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 $\tilde{m}$, see Appendix A). Then, we find the resonant equilibrium point (i.e. the value of $\mathcal{L}$) that corresponds to a value of a₃ ∕a₁ equal to the observed semi-major axis ratio, and therefore put ${(a_3/a_1)|}_{\text{eq}}\equiv {(a_3/a_1)|}_{\text{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₂ ∕a₁ and a₃ ∕a₂. 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₂∕a₁ 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 ${(a_3/a_1)|}_{\text{eq}}\equiv {(a_3/a_1)|}_{\text{obs}}$ automatically fixes the equilibrium values of both semi-major axis ratios a₂ ∕a₁ and a₃ ∕a₂. Considering for example the corresponding equilibrium ${(a_3/a_2)|}_{\text{eq}}$, we obtain the weighted difference $$\overline{\delta }(a_3/a_2):=\frac{({(a_3/a_2)|}_{\text{eq}}-{(a_3/a_2)|}_{\text{obs}})}{{(a_3/a_2)|}_{\text{obs}}}$$(16)

between ${(a_3/a_2)|}_{\text{eq}}$ and the observed value ${(a_3/a_2)|}_{\text{obs}}$. The same can be done for a₂∕a₁, which gives a similar (absolute) result, $|\overline{\delta }(a_2/a_1)|\simeq |\overline{\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 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₁ ∕m₂ = m₂∕m₃ = 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 ${(a_3/a_1)|}_{\text{eq}}\equiv {(a_3/a_1)|}_{\text{obs}}$; finally, we obtain $\overline{\delta }(a_3/a_2)$ as defined in Eq. (16) (and similarly for $\overline{\delta }(a_2/a_1)$). As explained above, imposing ${(a_3/a_1)|}_{\text{eq}}\equiv {(a_3/a_1)|}_{\text{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 maxe_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 $\overline{\delta }(a_3/a_2)$. Panels a and b in Fig. 6 are the result of this procedure spanning different planetary mass ratios, showing maps of $\overline{\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 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 Kepler-31, 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 first-order mean motion resonances by placing the planets relatively wide of the desired resonances, according to the specific values of kⁱⁿ 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 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 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₁ ∕m₂ and m₂ ∕m₃ 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 semi-major 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 disc-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. 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 semi-major 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⁻² for a typical planetary mass of order 10⁻⁵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 $\mathcal{L}$, so that these results are in fact generalisable to any dissipative process such that Ė < 0 and $\stackrel{.}{\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₃∕a₁ and a₃∕₂ will increase in time. When, for two consecutive outputs of the simulation, the ratio a₃ ∕a₁ crosses the observed one ${(a_3/a_1)|}_{\text{obs}}$, we store the corresponding value of a₃∕a₂ 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₃∕a₂ values. Then, we report the average of this list, and again compare this quantity with the observed ${(a_3/a_2)|}_{\text{obs}}$ (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₂ ∕a₁ (which gives again similar values to that of a₃∕a₂ as we mentioned in Sect. 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 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 $\overline{\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 ${(a_3/a_1)|}_{\text{obs}}$, but the typical value of $\overline{\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 disc-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-disc 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₁ ∕m₂, m₂∕m₃ ≳ 1) the numerical results agree very well with the analytical ones; moreover we still observe that the theoretical value of $\overline{\delta }(a_3/a_2)$ varies extremely little with m₁∕m₂ and m₂∕m₃ (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 ai+1∕ai vs. $e_i/\tilde{m}$ planes, i = 1, 2, for three equal-mass planets in a 3:2 – 2:1 mean motion resonance chain, close to which Kepler-305 resides. Orange vertical lines show the exact nominal commensurability, while dashed vertical lines show the observed a2 ∕a1 and a3 ∕a2 in the case of Kepler-305. As explained in the text, we select one equilibrium configuration (indicated by the red dot in both panels) by requiring that (a3 ∕a1)|eq ≡ (a3∕a1)|obs, which automatically fixes all orbital elements a2∕a1, a3 ∕a2, e1, e2 and e3 along the equilibrium curves.

Fig. 6 Top row: analytical maps constructed for Kepler-305 as explained in Sect. 3.2. Panel a: we plot the quantity $\overline{\delta }(a_3/a_2)$, which represents how close the system is now to some resonant equilibrium point, for different mass ratios m1 ∕m2 and m2 ∕m3 (each point in this plot is constructed by repeating the procedure described in Fig. 5). We notice that $\overline{\delta }(a_3/a_2)$ 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 $\text{max\hspace{0.17em}}{e}_{\text{eq}}/\tilde{m}$ representing the equilibrium orbital configuration that is selected at each fixed value of m1 ∕m2 and m2 ∕m3 by imposing (a3∕a1)|eq ≡ (a3∕a1)|obs. Bottom row: numerical maps constructed for Kepler-305 as explained in Sect. 3.3. We show numerical maps of $\overline{\delta }(a_3/a_2)$ in panel c and $\text{max\hspace{0.17em}}{e}_{\text{eq}}/\tilde{m}$ 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 $\overline{\delta }(a_3/a_2)~0.01$ 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 Kepler-31. 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 model-dependent 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 $\overline{\delta }(a_3/a_2)~0.003$. Comparing with Fig. 9 we see that there is a ~20% probability that Kepler-31 lies this close to the 2:1 – 2:1 chain by pure chance.

4 Results

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

Using the maps of $\overline{\delta }(a_3/a_2)$ shown in Figs. 6–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 $\overline{\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₁ ∕m₂ and m₂ ∕m₃ 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₁ ∕m₂ and m₂ ∕m₃ spanning reasonable values. Note also that for all three systems $\overline{\delta }(a_3/a_2)$ is small, but never vanishing, which would represent an analytically computed equilibrium configuration such that $\overline{\delta }(a_3/a_2)=0$, that is, a resonant equilibrium point which satisfies ${(a_3/a_2)|}_{\text{eq}}={(a_3/a_2)|}_{\text{obs}}$ and ${(a_2/a_1)|}_{\text{eq}}={(a_2/a_1)|}_{\text{obs}}$. But even if this happened to be the case, the level curve $\overline{\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₂∕m₁ or m₃∕m₂ using observed data, in general, inconclusive.

Secondly, we note that we do obtain in all three cases small values for $\overline{\delta }(a_3/a_2)$, meaning that these systems are indeed close to some equilibrium point of the Hamiltonian (Eq. (13)) and therefore could potentially reside in a multi-resonant chain. However, these values by themselves do not contain any meaningful information. The quantity $\overline{\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 < Δ < 0.05 for each pair, and extract the corresponding $\overline{\delta }(a_3/a_2)$ (calculated for the choice of mass ratios m₁∕m₂ = m₂∕m₃ = 1 for simplicity, since as we saw above the $\overline{\delta }(a_3/a_2)$ value depends extremely weakly on the mass ratios). From the cumulative distribution of $|\overline{\delta }(a_3/a_2)|$ that arises from this procedure we can obtain the probability that any given system has a given (small) $\overline{\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 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 $\overline{\delta }(a_3/a_2)$. Instead, Kepler-305 and Kepler-31 are likely to be in resonance at the 1σ level (i.e. the probability that their value of $\overline{\delta }(a_3/a_2)$ is smaller than the determined value by chance is less than 32%): for Kepler-305 there is a ~15% chance that this particular system lies this close to resonance by chance, while for Kepler-31 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 three-body Laplace angle $\stackrel{.}{{\phi }_{\text{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 φ_L = 2λ₁ − 4λ₂ + 2λ₃ satisfies ${\dot{\phi }}_{\text{L}}\simeq {0.5}^{°} {\text{days}}^{-1}$ given the observed transits periods; for Kepler-31, φ_L = λ₁ − 3λ₂ + 2λ₃ satisfies ${\dot{\phi }}_{\text{L}}\simeq {0.1}^{°} {\text{days}}^{-1}$; finally for YZ Cet, φ_L = 2λ₁ − 5λ₂ + 3λ₃ satisfies ${\dot{\phi }}_{\text{L}}\simeq {9.4}^{°} {\text{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 $\overline{\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.

Fig. 9 Cumulative distribution functions for $|\overline{\delta }(a_3/a_2)|$ 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 $\overline{\delta }(a_3/a_2)$. From this, we obtain on the vertical axis the probability that these systems could have this value of $\overline{\delta }(a_3/a_2)$ by pure chance.

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 ≃5.49, ≃ 8.29 and ≃ 16.74 days respectively. Their masses have been constrained via TTV to be ≃4 M_⊕ 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}^{°} {\text{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 Sect. 3.2, by imposing a₃ ∕a₁ to be equal to the observed ${(a_3/a_1)|}_{\text{obs}}$. This gives $\overline{\delta }(a_3/a_2)$ of order 4 × 10⁻⁵. Using an analogous argument to that of Fig. 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₁ ≃ 0.02, e₂ ≃ 0.03, e₃ ≃ 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 |Δ| < 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 Δ > 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 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 nearly-resonant 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 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 disc-planet 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 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 $\overline{\delta }(a_3/a_2)$ defined in Eq. (16), 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 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 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 ~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 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 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 mean-motion resonance, with nominal semi-major axes ā 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 Sect. 2. Having fixed the planet–planet mass ratios m₁ ∕m₂ = β₁ and m₂ ∕m₃ = β₂, we introduce the average planet–star mass ratio $$\tilde{m}=\frac{m_1+m_2+m_3}{3M_{*}}=\frac{m_1(1+{\beta }_1^{-1}+{\beta }_1^{-1}{\beta }_2^{-1})}{3M_{*}}.$$(A.1)

Inverting this expression we easily get all planetary masses in terms of $\tilde{m}$, $$\begin{array}{lll}\hfill & \hfill & m_1=c_1\tilde{m}:=\frac{3{\beta }_1{\beta }_2{M}_{*}}{1+{\beta }_2+{\beta }_1{\beta }_2}\tilde{m},\hfill \\ \hfill & \hfill & m_2=c_2\tilde{m}:=\frac{3{\beta }_1{M}_{*}}{1+{\beta }_2+{\beta }_1{\beta }_2}\tilde{m},\hfill \\ \hfill & \hfill & m_3=c_3\tilde{m}:=\frac{3M_{*}}{1+{\beta }_2+{\beta }_1{\beta }_2}\tilde{m},\hfill \end{array}$$(A.2)

with coefficients c depending on M_*, β₁ and β₂ only.

We introduce the modified Delaunay action-angle variables (Λ, Γ, λ, γ) as in Eq. (5), but we rescale the actions by the common mass factor $\tilde{m}$: this gives the following definition for the Λ’s (we maintain the same notation as the non-rescaled actions for simplicity) $$\begin{array}{lll}\hfill & \hfill & {\Lambda }_1=\frac{3{\beta }_1{\beta }_2{M}_{*}}{1+{\beta }_2+{\beta }_1{\beta }_2}\sqrt{\mathcal{G}{M}_{*}{a}_1}=c_1\sqrt{\mathcal{G}{M}_{*}{a}_1},\hfill \\ \hfill & \hfill & {\Lambda }_2=\frac{3{\beta }_1{M}_{*}}{1+{\beta }_2+{\beta }_1{\beta }_2}\sqrt{\mathcal{G}{M}_{*}{a}_2}=c_2\sqrt{\mathcal{G}{M}_{*}{a}_2},\hfill \\ \hfill & \hfill & {\Lambda }_3=\frac{3M_{*}}{1+{\beta }_2+{\beta }_1{\beta }_2}\sqrt{\mathcal{G}{M}_{*}{a}_3}=c_3\sqrt{\mathcal{G}{M}_{*}{a}_3},\hfill \end{array}$$(A.3)

and the same formal definition of Γ = Λe²∕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 $\tilde{m}$, in order for the Hamilton equations to be conserved we must also rescale the Hamiltonian itself by $\tilde{m}$. As in Sect. 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 Hamiltonianon the nominal values $\overline{\Lambda }$ defined fromā using Eq. (A.3). It is then easy to see that the rescaled Keplerian Hamiltonian takes the form $${\mathcal{H}}_{\text{kepl}}=-{\displaystyle \sum _{i=1}^3}\frac{c_i^3}{2}{\left(\frac{\mathcal{G}{M}_{*}}{{\Lambda }_i}\right)}^2,$$(A.4)

and is therefore independent of $\tilde{m}$, while the rescaled resonant part will have a multiplicative coefficient $\tilde{m}$: $$\begin{array}{lll}{\mathcal{H}}_{\text{res}}\hfill & =\hfill & \tilde{m}[-\frac{{\mathcal{G}}^2{M}_{*}{c}_1{c}_2^3}{{\overline{\Lambda }}_2^2}\hfill \\ \hfill & \hfill & \times (f_{\text{res}}^{(1,\text{in})}\sqrt{\frac{2{\Gamma }_1}{{\overline{\Lambda }}_1}}\cos (k^{\text{in}}{\lambda }_2-(k^{\text{in}}-1){\lambda }_1+{\gamma }_1)\hfill \\ \hfill & \hfill & +f_{\text{res}}^{(2,\text{in})}\sqrt{\frac{2{\Gamma }_2}{{\overline{\Lambda }}_2}}\cos (k^{\text{in}}{\lambda }_2-(k^{\text{in}}-1){\lambda }_1+{\gamma }_2))+\hfill \\ \hfill & \hfill & -\frac{{\mathcal{G}}^2{M}_{*}{c}_2{c}_3^3}{{\overline{\Lambda }}_3^2}\hfill \\ \hfill & \hfill & \times (f_{\text{res}}^{(1,\text{out})}\sqrt{\frac{2{\Gamma }_2}{{\overline{\Lambda }}_2}}\cos (k^{\text{out}}{\lambda }_3-(k^{\text{out}}-1){\lambda }_2+{\gamma }_2)\hfill \\ \hfill & \hfill & +f_{\text{res}}^{(2,\text{out})}\sqrt{\frac{2{\Gamma }_3}{{\overline{\Lambda }}_3}}\cos (k^{\text{out}}{\lambda }_3-(k^{\text{out}}-1){\lambda }_2+{\gamma }_3))].\hfill \end{array}$$(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² one can see the following. Assuming a fixed equilibrium value of the semi-major axes (that is, of the Λ’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)}_{\text{eq}}$ and ${(a_3/a_2)}_{\text{eq}}$ remain independent of $\tilde{m}$.

References

All Figures

Fig. 1 Observations of planet-hosting stars reveal that multi-planetary 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 multi-planetary 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 m1 = m2 = m3 = 10−4M*. 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 semi-major axis ratio. The location of the nominal resonant semi-major axis ratio (3/2)2/3 is shown bya vertical orange line. We also superimpose the numerically computed evolution of a three-planet 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 $\Delta =\frac{k-1}{k}\frac{P_2}{P_1}-1$ with k =2, 3 in all exoplanetary systems (yellow) and for three-planet 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 three-planets 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 three-planets 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 ai+1∕ai vs. $e_i/\tilde{m}$ planes, i = 1, 2, for three equal-mass planets in a 3:2 – 2:1 mean motion resonance chain, close to which Kepler-305 resides. Orange vertical lines show the exact nominal commensurability, while dashed vertical lines show the observed a2 ∕a1 and a3 ∕a2 in the case of Kepler-305. As explained in the text, we select one equilibrium configuration (indicated by the red dot in both panels) by requiring that (a3 ∕a1)|eq ≡ (a3∕a1)|obs, which automatically fixes all orbital elements a2∕a1, a3 ∕a2, e1, e2 and e3 along the equilibrium curves.

In the text

Fig. 6 Top row: analytical maps constructed for Kepler-305 as explained in Sect. 3.2. Panel a: we plot the quantity $\overline{\delta }(a_3/a_2)$, which represents how close the system is now to some resonant equilibrium point, for different mass ratios m1 ∕m2 and m2 ∕m3 (each point in this plot is constructed by repeating the procedure described in Fig. 5). We notice that $\overline{\delta }(a_3/a_2)$ 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 $\text{max\hspace{0.17em}}{e}_{\text{eq}}/\tilde{m}$ representing the equilibrium orbital configuration that is selected at each fixed value of m1 ∕m2 and m2 ∕m3 by imposing (a3∕a1)|eq ≡ (a3∕a1)|obs. Bottom row: numerical maps constructed for Kepler-305 as explained in Sect. 3.3. We show numerical maps of $\overline{\delta }(a_3/a_2)$ in panel c and $\text{max\hspace{0.17em}}{e}_{\text{eq}}/\tilde{m}$ 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 $\overline{\delta }(a_3/a_2)~0.01$ 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 Kepler-31. 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 model-dependent 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 $\overline{\delta }(a_3/a_2)~0.003$. Comparing with Fig. 9 we see that there is a ~20% probability that Kepler-31 lies this close to the 2:1 – 2:1 chain by pure chance.

In the text

Fig. 9 Cumulative distribution functions for $|\overline{\delta }(a_3/a_2)|$ 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 $\overline{\delta }(a_3/a_2)$. From this, we obtain on the vertical axis the probability that these systems could have this value of $\overline{\delta }(a_3/a_2)$ by pure chance.

In the text