© The Authors 2024

1 Introduction

Thousands of exo-planetary systems, over one-fourth of which are multi-planetary systems, have been detected. According to their orbital architectures, these multi-planetary systems are categorised into separable, compact, and resonant planetary systems. The separable (secular) systems are protected from a mutual collision and a close encounter by the low angular momentum deficit (AMD), that is, the AMD stability (Laskar & Petit 2017). The compact multi-planetary systems, however, even for those with circular orbits, are generally unstable. Petit et al. (2020) show how the slow chaotic diffusion due to the overlap of three-body resonances dominates the timescale leading to the instability for initially coplanar and circular orbits, and this type of mechanism reproduces the qualitative behaviour found in numerical simulations very well (Hussain & Tamayo 2020). The resonant systems, in which planets are found in mean motion resonance (MMR) with each other, are likely to be stable and of particular interest because they may bear important clues as to the formation and evolution of planetary systems. For the case of two planets in an MMR, the dynamics are well understood through the resonant periodic orbits (e.g. Voyatzis 2008; Voyatzis et al. 2009) and the apsidal corotational resonance (ACR) (e.g. Michtchenko et al. 2008). But the resonant structure could be much richer and complex in multi-planetary systems.

Except for the two-body MMR involving two planets and the hosting star, the three-body resonance (three planets are involved) also plays a role (e.g. Petit 2021; Cerioni et al. 2022). Particularly, a multi-planetary system with orbital periods of neighbouring planets librating at a sequence of near-integer ratios (2/1, 3/2, 4/3, etc.), often referred to as a resonant chain, is of great interest. The well-known examples include three planets of Kepler-51 in a 1:2:3 resonant chain (Masuda 2014), four planets in a 1:2:4:8 resonant chain in HR8799 (Goździewski & Migaszewski 2020) and a 1:2:5:7 resonant chain in K2-32 (Heller et al. 2019), five planets of TOI-178 in a 2:4:6:9:12 (Leleu et al. 2021) and a 4:6:9:12:18 resonant chain in Kepler-80 (MacDonald et al. 2021), and another five planets composing a resonant chain of 3:2 MMRs in K2-138 (Christiansen et al. 2018; Lopez et al. 2019). In these systems, in addition to the transit light curves, a transit-timing variation (TTV), and radial velocity data, the dynamics of the resonant chain have also been taken into account to improve the planetary parameters (mass, orbital period, etc.).

Plenty of numerical studies have been devoted to the formation history and dynamical stability of multi-planetary systems with planets locked in or close to a resonant chain. At the early stage of a planetary system’s formation, planets can be captured into a resonant chain through the convergent orbital migration. The tidal effect and the accretion of matter (increasing planetary mass) may influence the further evolution, driving the planetary system away from the nominal positions of resonances, and destabilising the resonant chain via the secondary resonance. The final surviving systems exhibit complex resonant structures involving in a two-body MMR, three-body MMR, resonant chain, secular resonance, etc. (e.g. MacDonald & Dawson 2018; Pichierri et al. 2019; Morrison et al. 2020; Pichierri & Morbidelli 2020; MacDonald et al. 2022; Wong & Lee 2024).

Due to the high degree of freedom, the analytical method is difficult to be applied to such systems. Delisle (2017) proposes an analytical model for any resonant chain and successfully predicts the resonant configuration of Kepler-223 (Mills et al. 2016). And recently, Antoniadou & Voyatzis (2022) extended the methodology for a pair of resonant exoplanets in coplanar orbits (Antoniadou & Voyatzis 2016; Antoniadou 2016) to the case of triple massive coplanar planets in a resonant chain (general four-body problem) and, particularly, computed the resonant periodic orbits for the Kepler-51 system. Based on the stable periodic orbits of the 1:2:3 resonant chain, they demonstrate three possible scenarios that can safeguard the Kepler-51 planetary system. The comprehensive study on the stable periodic orbits can provide an optimum deduction of orbital elements, in addition to the observational data fitting.

In this paper, we generalise the semi-analytical model for analysing the periodic orbits in the coplanar two-body MMR first developed by Michtchenko et al. (2006) to the case of the resonant chain. For the sake of explaining our methods clearly, we mainly use the 1:2:3 resonance chain in the Kepler-51 system as an example in our analyses. But we emphasise here that the methods and some conclusions presented in this paper are general and can be applied to any other resonant chains of three or even more planets. In fact, we also briefly show our analyses on some other resonant chains in the end of this paper.

The rest part of this paper is organised as follows. In Sec. 2, we revisit the definition of resonant periodic orbits and generalise the definition to arbitrary resonant chains. We introduce our semi-analytical Hamiltonian method and the determination of ACR in the resonant chain in Sec. 3. Then this methodology is applied to the case of the 1:2:3 resonant chain, and the periodic orbits, that are, ACRs are obtained. In Sec. 4, the orbital dynamics around the stable ACRs are investigated. We study how the orbital configuration of planets in a resonant chain can be achieved via convergent planetary migration in Sec. 5. Finally, we draw our conclusion in Sec. 6.

2 The model

2.1 Rotational frame and resonant periodic orbits

The motion of planets in a resonant chain is often described in a rotational frame. Consider three massive planets m₁, m₂, m₃ from the inside out orbiting around the central star M on the same plane, with their orbital periods T₁, T₂, T₃ satisfying the commensurable condition $$T_1:T_2:T_3=k_1:k_2:k_3,\text{ with }{k}_1,k_2,k_3\in \mathcal{Z}.$$(1)

Following Antoniadou & Voyatzis (2022), a non-inertial frame (GXY) is introduced as follows. Let G coincide with the mass centre of M and m₁, the axis GX always points to m₁ from M, and GY be perpendicular to GX. In this rotating frame, m₁ and M are always on the GX axis, while m₂ and m₃ revolve about m₁ with the relative frequencies of n₂ − n₁ and n₃ − n₁, where n₁, n₂, n₃ are the mean motions of three planets. The relative frequencies satisfy the near-integer commensurability: $$\frac{n_2-n_1}{n_3-n_1}=\frac{p}{q},\text{ with }p,q\in \mathcal{Z}.$$(2)

The period T of the system, during which m₂ revolves p times while m₃ revolves q times both around m₁, can be easily derived as $$T=\frac{T_1}{1-T_1/T_2}p=\frac{T_1}{1-T_1/T_3}q.$$(3)

In practice, such periodicity requires that the planets m₁, m₂, and m₃ repeat their configuration after time T. Therefore, the period T should be related to the least common multiple (LCM) of T₁, T₂, T₃, and this type of periodic condition can be rewritten in a more general form $$T/T_i=\mathrm{LCM}(k_1,k_2,k_3)/k_i,\text{ for }i=1,2,3.$$(4)

For the 1:2:3 resonant chain in the planetary system Kepler-51, Eq. (4) gives LCM(1, 2, 3) = 6, and T = 6T₁ = 3T₂ = 2T₃, which is equivalent to p = 3, q = 4 and T = 6T₁ (Antoniadou & Voyatzis 2022). The periodicity defined in Eq. (4) can be generalised to resonant chain of N (N > 3) planets (e.g. a four-planet resonant system Kepler-223 can be found in Appendix B) as $$T/T_i=\mathrm{LCM}(k_1,k_2,\cdots ,k_N)/k_i,\text{ for }i=1,2,\cdots ,N.$$(5)

These periodic orbits can be found in the exact (nominal) location of the MMRs. In addition, according to different orbital configurations, these periodic orbits can be divided into two groups, namely the symmetric periodic orbits and the asymmetric ones. A periodic orbit is symmetric with respect to the GX axis of the rotating frame if it remains invariant under the fundamental symmetry (Hénon 2003) ∑: (t, X, Y) → (−t, X, −Y) $$\begin{array}{rl}& X_i(-t)=X_i(t),Y_i(-t)=-Y_i(t),\\ & {\dot{X}}_i(-t)=-{\dot{X}}_i(t),{\dot{Y}}_i(-t)={\dot{Y}}_i(t),i=1,2,3.\\ & \end{array}$$(6)

This indicates that all the two-body MMRs should take the symmetric configuration. Otherwise, if one of the symmetry condition breaks, i.e. Eq. (6) is not satisfied for any i, then the configuration is referred as asymmetric. The relationship between symmetric condition and the choice of the resonant angles can also be found in e.g. Voyatzis (2008).

2.2 Hamiltonian description

The Hamiltonian of the planar three-planet (from the inside out m₁, m₂, m₃) system takes the form (see e.g. Laskar 1991) $$\begin{array}{rl}\mathcal{H}& ={\mathcal{H}}_0+{\mathcal{H}}_1\\ {\mathcal{H}}_0& =\sum _{i=1}^3(\frac{{\mathit{p}}_i^2}{2M}-\mathcal{G}\frac{M{m}_i}{\left|{\mathit{r}}_i\right|}),\\ {\mathcal{H}}_1& =\sum _{1\le i<j\le 3}\frac{{\mathit{p}}_i\cdot {\mathit{p}}_j}{M}-\frac{\mathcal{G}{m}_i{m}_j}{|{\mathit{r}}_i-{\mathit{r}}_j|},\end{array}$$(7)

where $\mathcal{G}$ is the gravitational constant, r_i is the position vector of the i-th planet, and p_i is the canonically conjugated momentum. The $\mathcal{H}$₀ is the sum of Keplerian part of each planet, while $\mathcal{H}$₁ describes the interaction between planets. The Poincaré canoncial variables in astrocentric coordinate are $$\{\begin{array}{ll}{\mathrm{\Lambda }}_i={\beta }_i\sqrt{{\mu }_i{a}_i},& {\lambda }_i\\ {\mathrm{\Gamma }}_i={\mathrm{\Lambda }}_i(1-\sqrt{1-e_i^2}),& {\gamma }_i=-{\varpi }_i,\end{array}$$(8)

where μ_i = $\mathcal{G}$(M + m_i), β_i = m_iM/(m_i + M), and a_i, e_i, ϖ_i, λ_i are the semi-major axis, eccentricity, longitude of periastron and the mean longitude of the i-th planet m_i, respectively. The action variable Λ_i is the circular angular momentum and Γ_i is often referred as the AMD (see e.g. Petit et al. 2017).

For two planets m_i, m_j with mean motions n_i, n_j, we denote by k_ji/k_ij = n_i/n_j = T_j/T_i the relative period ratio, and q_ij = q_ji = |k_ji − k_ij| the order of the resonance, where k_ij and k_ji are integers. The two-body resonant angles θ_ij, Δϖ_ij are defined as $${\theta }_{ij}=k_{ji}{\lambda }_j-k_{ij}{\lambda }_i-(k_{ji}-k_{ij}){\varpi }_i,\mathrm{\Delta }{\varpi }_{ij}={\varpi }_i-{\varpi }_j.$$(9)

For planetary systems near MMRs, Delisle (2017) introduces the canonical resonant variables for arbitrary resonant chain (with any number of planets, in resonances of any degree). Particularly, the resonant variables of the 1:2:3 resonant chain are $$\{\begin{array}{ll}{\sigma }_1=3{\lambda }_3-2{\lambda }_2-{\varpi }_1,& {\mathrm{\Gamma }}_1,\\ {\sigma }_2=3{\lambda }_3-2{\lambda }_2-{\varpi }_2,& {\mathrm{\Gamma }}_2,\\ {\sigma }_3=3{\lambda }_3-2{\lambda }_2-{\varpi }_3,& {\mathrm{\Gamma }}_3,\\ {\phi }_0=3{\lambda }_3-4{\lambda }_2+{\lambda }_1,& L_0={\mathrm{\Lambda }}_1,\\ {\phi }_G=3{\lambda }_3-2{\lambda }_2,& G=\sum _{i=1}^3({\mathrm{\Lambda }}_i-{\mathrm{\Gamma }}_i,)\\ Q=Q_{23}={\lambda }_2-{\lambda }_3,& J=6{\mathrm{\Lambda }}_1+3{\mathrm{\Lambda }}_2+2{\mathrm{\Lambda }}_3.\end{array}$$(10)

The first three angles in Eq. (10) are resonant angles conjugated to each planet’s AMD, and φ₀ is the angle of Laplace resonance (0-th) between three planets. Since φ_G does not appear explicitly in the Hamiltonian (see details e.g. in Delisle 2017), the corresponding conjugated momentum G is constant, implying the preservation of angular momentum. Thus, the problem is in fact of five degree-of-freedom (DoF). The last angle Q is the fast angle with respect to the first four resonant angles (σ₁, σ₂, σ₃,φ₀), thus it can be removed by numerical averaging (see e.g. Gauss-Legendre formula, Press et al. 1996). The averaging process guarantees that J is another constant, which is similar to the ‘spacing parameter’ in the planar two-body MMR. Moreover, the integral constant J in triple-planet resonant chain T₁: T₂: T₃ = k₁: k₂: k₃ can be rewritten in a more general form L (hereafter we refer L as the spacing parameter in a resonant chain) $$L=\sum _{i=1}^3\frac{{\mathrm{\Lambda }}_i}{k_i}.$$(11)

For the k₁: k₂: k₃ = 1: 2: 3 resonant chain, during the period T (see the definition in Eq. (4)), the angle Q evolves $$Q^T=2\pi \times \mathrm{LCM}(k_1,k_2,k_3)(\frac{1}{k_2}-\frac{1}{k_3})=2\pi .$$(12)

Therefore, the averaged Hamiltonian $\overline{\mathcal{H}}$ can be defined as $$\overline{\mathcal{H}}=\frac{1}{Q^T}{\int }_0^{Q^T}\mathcal{H}\mathrm{d}Q.$$(13)

The resonant periodic orbits will appear as stationary solution of the averaged Hamiltonian Eq. (13). After the numerical averaging process, the final Hamiltonion $\overline{\mathcal{H}}$ is of 4 DoF in canonical variables (σ₁, σ₂, σ₃, φ₀, Γ_1; Γ₂, Γ₃, L₀) parametrised by the two constants (L, G).

Seemingly, the resonance between m₂ and m₃ plays an essential role in Eq. (10), but it is only due to the arbitrary choice of the canonical variables. For any different choice of resonant angles, the expression of the synodic angle Q would be different. For example, the choices Q = Q₁₂ = λ₁ − λ₂, or Q = Q₁₃ = λ₁ − λ₃, the averaged Hamiltonian should be respectively calculated by: $$\begin{array}{rl}\overline{\mathcal{H}}& =\frac{1}{Q_{12}^T}{\int }_0^{Q_{12}^T}\mathcal{H}\mathrm{d}{Q}_{12}=\frac{1}{Q_{13}^T}{\int }_0^{Q_{13}^T}\mathcal{H}\mathrm{d}{Q}_{13},\\ Q_{12}^T& =2\pi \times \mathrm{LCM}(k_1,k_2,k_3)(\frac{1}{k_1}-\frac{1}{k_2})=6\pi ,\\ Q_{13}^T& =2\pi \times \mathrm{LCM}(k_1,k_2,k_3)(\frac{1}{k_1}-\frac{1}{k_3})=8\pi .\end{array}$$(14)

For a better characterisation of the resonant configuration, we perform the linear transformation to angular variables defined in Eq. (10), and the resonant angles (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) between consecutive planets can be obtained as: $$\{\begin{array}{ll}{\theta }_{12}=2{\lambda }_2-{\lambda }_1-{\varpi }_1,& \mathrm{\Delta }{\varpi }_{12}={\varpi }_1-{\varpi }_2,\\ {\theta }_{23}=3{\lambda }_3-2{\lambda }_2-{\varpi }_2,& \mathrm{\Delta }{\varpi }_{23}={\varpi }_2-{\varpi }_3.\end{array}$$(15)

The first two angles θ₁₂, Δϖ₁₂ are for the MMR between the inner two planets m₁, m₂ while θ₂₃, Δϖ₂₃ are for the outer planets m₂, m₃ respectively, just as defined in Eq. (9).

For isolated two-body MMRs, the resonant angles (θ_ij, Δϖ_ij) are all 2π-periodic, i.e., the geometric configuration of MMR characterised by (θ_ij, Δϖ_ij) is identical to the one by (θ_ij + 2k′π,Δϖ_ij + 2k″π), with k′,k″ ∈ $\mathcal{Z}$. For the 1:2:3 resonant chain, the averaged Hamiltonian Eq. (13) is 2π-periodically dependent on all the angular variables defined in Eq. (15). According to the symmetry in the co-rotational frame GXY with respect to m₁, the configuration, in which the resonant angles (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) take either 0 or π, is symmetric, while other configurations are asymmetric. It should be noted that the above conditions may not always guarantee the proper periodicity and symmetric configuration for a general three-planet resonant chain if we choose the resonant angles between adjacent planets (rather θ_i,i+1, Δϖ_i,i+1 than θ_i,i+2, Δϖ_i,i+2) to characterise the system. We leave a detailed discussion for this problem to Appendix A.

We note that the aforementioned procedure can be simply extended to systems containing N (N > 3) planets, in which the two constants take the form: $$\begin{array}{rl}G& =\sum _{i=1}^N{\beta }_i\sqrt{{\mu }_i{a}_i(1-e_i^2)},\\ L& =\sum _{i=1}^N\frac{{\beta }_i\sqrt{{\mu }_i{a}_i}}{k_i}.\end{array}$$(16)

And the resonant configuration for an N-planet system can be represented by (θ₁₂, Δϖ₁₂, ⋯, θ_n−1,n, Δϖ_n−1,n).

Now, searching for the resonant periodic orbits is equivalent to finding the stationary solution of the averaged Hamiltonian with two free constants L and G.

Fig. 1 Topology of the Hamiltonian on different representative planes, (a) The topology on period ratio plane (n1/n2, n2/n3) with fixed eccentricities. (b) The topology on eccentricity plane (e1,e2). The stationary solution located at (n1/n2, n2/n3) = (2, 1.5) and (e1, e2, e3) = (0.00458, 0.00938, 0.00564) is marked by a red dot on both panels. The corresponding critical angles are (θ12, Δϖ12, θ23, Δϖ23) = (π, 0, π, π). The colour indicates the value of Hamiltonian, light yellow for high and dark blue for low. The same colour code is adopted in all figures in this paper, unless otherwise stated.

3 Topology of Hamiltonian and stationary solutions

The averaged Hamiltonian Eq. (13) has been reduced to 4 DoF. Due to the complex structure in high dimensional phase space (eight-dimensional phase space for aplanar three-planet system), searching for stationary solutions in the entire phase space is very difficult. Fortunately, the topology of the Hamiltonian $\overline{\mathcal{H}}$ can be understood from several representative planes (section of phase space), namely, the period ratio plane (n₁/n₂, n₂/n₃), the eccentricity plane (e₁, e₂) or (e₂, e₃), and the angular planes (θ₁₂, Δϖ₁₂) and (θ₂₃, Δϖ₂₃).

Hereafter, as an example, we focus on a 1:2:3 resonant chain with the mass parameters of M = 1.0, m₁ = 2.1 × 10⁻⁶, m₂ = 4.0 × 10⁻⁶, m₃ = 7.6 × 10⁻⁶. The mass ratios between the planet pairs are ρ_i = m₁/m₂ = 0.525 and ρ_o = m₂/m₃ = 0.526 (subscripts “i” for inner and “o” for outer), which in fact is similar to the ones in the Kepler-51 system (Masuda 2014). We set the gravitational constant $\mathcal{G}$ = 1 and fix the semi-major axes of m₂ as a₂ = 1 AU. As a result, the orbital period of m₂ is 1 yr.

3.1 Topology on period ratio plane (n₁/n₂, n₂/n₃)

In a three-planet system, the planetary semi-major axes (a₁, a₂, a₃) are constrained by the spacing parameter L, so that they can be uniquely determined by the period ratio pair (n₁/n₂, n₂/n₃) when L is given. Therefore, we can study the topology of the Hamiltonian on the period ratio plane (n₁/n₂, n₂/n₃) with e₁, e₂, e₃, θ₁₂, Δϖ₁₂, θϖ₂₃ being appropriately fixed.

An example of the topology on (n₁/n₂, n₂/n₃) plane is shown in Fig. la, for which the resonant angles (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (π, 0, π, π) and eccentricities (e₁, e₂, e₃) are fixed to (0.00458, 0.00938, 0.00564). The local extremum (maximum) is indicated by a red dot, and it is very close to the nominal value (n_1c/n_2c, n_2c/n_3c) = (2, 1.5). We note that such structure always holds for other parameters and even in other resonance chains, thus we can fix the period ratios at the nominal value (2, 1.5) as a first approximation and further discuss the topology on the eccentricity plane.

Fig. 2 Different topologies on the (e1, e2) plane for different symmetric configurations (a) (θ12, Δϖ12, θ23, Δϖ23) = (0, π, 0, π) and (b) (π, π, 0, 0). The local extrema are (e1, e2, e3) = (0.119,0.1352,0.119) (red solid circle) and (0.300, 0.0662, 0.0671) (red cross), respectively.

3.2 Topology on eccentricity plane (e₁, e₂)

Once the semi-major axes a₁, a₂, and a₃ have been determined, the conservation of total angular momentum G restricts the eccentricities e₁, e₂, and e₃ to a two-dimensional plane, and therefore we can consider the topology on a representative eccentricity plane, for example the (e₁, e₂) plane. An example is illustrated in Fig. 1b. From this plot, a node-like stationary solution located at (e₁ e₂) = (0.00458, 0.00938) can be found easily. We note that such stationary solutions are just ACRs.

In fact, for different resonant configurations, i.e. different values of (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃), as well as the total angular momentum G, the Hamiltonian has different structures on the (e₁, e₂) plane.

1. For resonant configuration (π, 0, π, π) as shown in Fig. 1, the stationary solution is the local maximum in the (n₁/n₂, n₂/n₃) plane, while is the minimum in the (e₁, e₂) plane, as indicated by colour of the contours. Therefore, the solution is a saddle-like point in the four-dimensional phase space (n₁/n₂, n₂ /n₃, e₁, e₂) and should be considered as unstable.

2. For resonant configuration (0, π, 0, π), with appropriate eccentricities, the topology on (n₁/n₂, n₂/n₃) plane is similar to Fig. 1a, but the topology on (e₁, e₂) plane, as shown by an example in Fig. 2a, has different features (from Fig. 1b). The stationary solution located at (e₁, e₂, e₃) = (0.119, 0.135, 0.119) now appears as the local maximum and is surrounded by robust energy level curves. The stability of this type ACR will be determined by the topological behaviour on the representative angular plane(s).

3. Except for the node-like structure on eccentricity plane, we also observe the saddle-like ACR for resonant configuration (π,π, 0,0) as shown in Fig. 2b. It is worth noting that such an ACR appears only for high eccentricities e₁ ~ 0.3. Since the saddle-like structure always reveals the orbital instability, these ACRs will not be taken into our consideration.

3.3 Symmetric ACRs

In a symmetric ACR, the critical angle could take either 0 or π. For the 1:2:3 resonant chain, the four critical angles make totally 2⁴ = 16 symmetric configurations. With fixed symmetric resonant angles, period ratios (n₁/n₂, n₂/n₃) and the total angular momentum G, we find and locate all the “node-like” (local minimum as in Fig. 1b and local maximum as in Fig. 2a) stationary solution (i.e. ACR) on the eccentricity plane. By varying G, we obtain the symmetric ACR families. The stability of an ACR solution can be determined by analysing the Hamiltonian structure in the two angular planes (θ₁₂, Δϖ₁₂) and (θ₂₃, Δϖ₂₃).

In Fig. 3, on the representative angular planes, we show examples of the Hamiltonian structure for small eccentricity region, where only the symmetric ACRs exist. The saddlelike resonant configuration (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (π, 0, π, π) marked by the red cross in Fig. 3a is unstable, while the node-like configuration (0, π, 0, π) marked as blue dots is stable.

We calculate the symmetric ACRs for small eccentricity and summarise them in Fig. 4. We note that similar results have been obtained by Antoniadou & Voyatzis (2022) using different methods. For the convenience of comparison, we adopt the same notations of the families and use the following resonant angles: $$\{\begin{array}{l}{\theta }_1=2{\lambda }_2-{\lambda }_1-{\varpi }_1,\\ {\theta }_2=2{\lambda }_2-{\lambda }_1-{\varpi }_2,\\ {\theta }_3=3{\lambda }_3-2{\lambda }_2-{\varpi }_2,\\ {\theta }_4=3{\lambda }_3-2{\lambda }_2-{\varpi }_3.\end{array}$$(17)

To distinguish the ACR families belonging to different configurations, we present these symmetric families on the (e₁ cos θ₁, e₂ cos θ₂) and (e₂ cos θ₂, e₃ cos θ₃) plane. The positive or negative value on the plane indicates that the corresponding resonant angle takes the value of 0 or π.

The ACRs we obtained in Fig. 4 are just the periodic orbits in the rotating frame. Starting from the circular orbit configuration, Antoniadou & Voyatzis (2022) find six symmetric families of periodic orbits, and they denote them by S_i, i = 1, ⋯, 6 (“S” for symmetric). We find four of them, namely S₂, S₃, S₅ and S₆, of which the critical angles are detailed in Table 1.

As for the two missing families, S₁ and S₄ in Antoniadou & Voyatzis (2022), their stable segments appear only very shortly when the eccentricities are very small (typically e_i ~ 10⁻⁶), which makes them less significant since a real planetary system is unlikely to be trapped in such a small-eccentricity configuration and survive for long time span. We failed to detect these two families probably because our semi-analytical Hamiltonian only contains the perturbation up to the first order in the planet-star mass ratio, and the contribution from high-order terms are neglected.

In addition, the symmetric family S₅ of configuration (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (π, 0, π, π) labelled as stable in Antoniadou & Voyatzis (2022) is topologically unstable in our calculation. In fact, the S₅ ACR is the local maximum on (n₁/n₂, n₂/n₃) plane but local minimum on (e₁, e₂) plane (e.g. Fig. 1), and such structure holds for the entire family. And we did not find any example that reaches such an orbital configuration in our numerical simulations of the convergent migration (see Sect. 5).

Fig. 3 Topology on representative angular planes. The eccentricities of the symmetric ACR are (e1, e2, e3) = (0.00458, 0.00938, 0.00564). (a) Hamiltonian level curves on (θ12, Δϖ12) plane with (θ23, Δϖ23) fixed to (π, π), the saddle-like point at (θ12, Δϖ12) = (π, 0) is marked by a red cross. (b) Topology on (θ23, Δϖ23) plane with (θ12, Δϖ12) fixed to (π, 0). The node-like point at (θ23, Δϖ23) = (π, π) is marked by a red dot. In addition, the stable symmetric configuration (θ12, Δϖ12, θ23, Δϖ23) = (0, π, 0, π) is shown as blue dots.

Fig. 4 Symmetric ACR families S2, S3, S5, and S6 in the 1:2:3 resonant chain for small-eccentricity configurations on (a) the (e1 cos θ1, e2 cos θ2) plane and (b) the (e2 cos θ2, e3 cos θ3) plane. The stable symmetric ACRs are in blue while the unstable ones in red. The notations of the families are the same as in Antoniadou & Voyatzis (2022), and the resonant angles of different configurations are presented in Table 1.

3.4 Bifurcation and asymmetric ACRs

As shown in Fig. 4, only the S₆ family is stable, and it loses the stability as the eccentricity increases. Along with the increasing eccentricity, the topology around the S₆ family varies in a complicated way. An example is illustrated in Fig. 5.

On the (θ₁₂, Δϖ₁₂) plane with fixed (θ₂₃, Δϖ₂₃) = (0, π), the node-like extreme point at (0, π) changes to a saddle-like point while two additional asymmetric ACRs, namely A₆ (“A” for asymmetric), arise as the eccentricity increases (Fig. 5a). Meantime, we note that the feature of the node-like structure around (0, π) on the (θ₂₃, Δϖ₂₃) plane with fixed (θ₁₂, Δϖ₁₂) = (0, π) remains the same (Fig. 5b).

The stable asymmetric ACRs are the local extrema of the full phase space (n₁/n₂, n₂/n₃, e₁, e₂, θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃). We adopt the geometrical approach (Michtchenko et al. 2006) to search these asymmetric ACRs in the eight-dimensional space.

All the ACRs are computed under two fixed constants L and G. Since the spacing parameter L only sets the global scale of the system and does not influence the resonant dynamics (except for changing the scales of distance, time, and energy), the evolution of ACRs is in fact parametrised by the constant G, or equivalently, the scaled AMD δ defined as $$\delta =\frac{{\mathrm{\Gamma }}_1+{\mathrm{\Gamma }}_2+{\mathrm{\Gamma }}_3}{L}.$$(18)

We summarise in Fig. 6 both the symmetric ACR S₆ and the asymmetric ACR A₆, in which the evolutions of eccentricities (e₁, e₂, e₃) and resonant angles along these ACR solutions are illustrated.

Along the A₆ family, the resonant angles θ₁₂, Δϖ₁₂ between the inner adjacent planets m₁, m₂ shift significantly from the (unstable) symmetric ACR, while the configuration of the outer planet pair (θ₂₃, Δϖ₂₃) experiences only small variations. Particularly, the angle Δϖ₂₃ remains tightly close to π as δ < 0.01.

Starting from near circular orbits (e_i ~ 0), the periodic solution extends to high eccentricity region along the stable branch, blue solid lines first and black dashed lines later in Figs. 6a, b. In a real system, planets in the resonant chain may not revolve exactly along these periodic solutions to high eccentricity, instead, they are expected to be around (but not exactly in) these solutions.

In practice, the numerically averaged Hamiltonian $\overline{\mathcal{H}}$ works for arbitrary planetary eccentricities because it does not depend on the literal expansion. Consequently, all the stable ACRs should appear as local extrema in the entire eight-dimensional phase space, and they can be numerically determined by the downhill simplex method under the constraints of fixed L and G. For the 1:2:3 resonant chain, we find more than ten families of stable ACRs at relatively high eccentricity region and they form very complicated structures on the (e₁, e₂) and (e₂, e₃) planes. However, further investigations reveal that the topology and the dynamical portraits around these ACR families are similar to each other. Therefore, hereafter we simply focus on the ACR families starting from near circular orbits, namely the S₆, A₆ family.

Fig. 5 Topology on representative angular planes for symmetric ACR S6 located at (e1e2,e3) = (0.0569, 0.0846, 0.0531). (a) Topology on (θ12, Δϖ12) plane with fixed (θ23, Δϖ23) = (0, π). The saddle-like point at (θ12, Δϖ12) = (0, π) marked by a cross indicates the instability of this solution, while two asymmetric ACRs bifurcated from this unstable symmetric one are indicated by black dots. (b) The topology on (θ23, Δϖ23) plane with (θ12, Δϖ12) fixed to (0, π). The node-like point at (θ23, Δϖ23) = (0, π) is marked by a pink dot.

Fig. 6 Locations and resonant configurations of the symmetric ACR family (S6) and the asymmetric ACR family (A6). In the upper two panels (a) and (b), the stable symmetric, unstable symmetric, and stable asymmetric ACRs are indicated by blue solid, red solid, and black dashed lines, respectively. Three solid circles on these lines indicate the initial conditions that will be analysed in detail later (see text). The lower two panels (c) and (d) show the evolutions of resonant angles along the asymmetric ACR as the function of δ (scaled AMD, see text). The resonant angles θ12, θ23 and the secular angles Δϖ12, Δϖ23 are labelled by the corresponding lines.

3.5 Law of structure

As we have shown, for the 1:2:3 resonant chain, the S₆−A₆ family is the only family that starts from circular orbits and then evolves to highly eccentric configuration. During the evolution from low to high eccentricity along the ACR solution, the mean motion ratios (n₁/n₂, n₂/n₃) of the adjacent planet pairs may change. Since these ACR solutions are found by searching for the extrema in the eight-dimensional full space, n₁/n₂ and n₂/n₃ are obtained automatically in our calculations. The evolution of these ratios as function of the planets’ eccentricities are summarised in Fig. 7.

As shown in Fig. 7, for the ACR solution, the ratios n_i/n_i+1 (i = 1, 2) are always above the nominal values (2 and 1.5 in the 1:2:3 resonant chain). Therefore, the planet pairs near the stationary resonant configuration, i.e. in resonance with small oscillations, always have orbital period ratio at one side of the nominal value. The deviation from the nominal value is large for near circular orbits and it decreases for larger eccentricities (thus larger AMD).

The mean motion ratio of the ACR solution depends on planetary masses. Except for the planetary masses as in the Kepler-51 system, we also test two additional systems, in which the planets are 10 and 100 times heavier than the ones in Kepler-51, i.e., keep the m_i/m_j (i, j = 1, 2, 3) values but set m_i/M 10 and 100 times larger than in the Kepler-51 system. Assume they are in the same 1:2:3 resonance chain, we calculate the mean motion ratios, and plot the results in Fig. 7, too. For near circular orbits, the stable ACR solution is symmetric (blue segments in Fig. 7), and from a critical e_i it is continued by asymmetric ACR (black segments). Our calculations demonstrate that the critical points for three cases of different total planet masses are the same (e₁, e₂, e₃) = (0.0185, 0.0336, 0.0186), implying the independence of the resonant configuration on the total planetary mass. Indeed, only the mass ratio between adjacent planets influences the resonant configuration (Antoniadou & Voyatzis 2022).

Apparently, for a resonant chain with given mass ratios between planets, the larger the total mass of planets the stronger the short-term perturbations the system would experience. Thus, the resonant chain with higher total planetary mass is more likely to be unstable. The main dynamical mechanism inducing such orbital instability is the 1:1 secondary resonance, which occurs when one of the libration period is almost equal to the period of synodic angle (Pichierri & Morbidelli 2020; Goldberg et al. 2022). This secondary resonance excites the orbital eccentricity and finally destabilises the orbits. In other words, the short term perturbation is too large to be averaged, if the total mass of planets is large enough.

Fig. 7 Law of structure for the 1:2:3 resonant chain. The left (right) column of panels show the relation between the mean motion ratio n1/n2 (n2/n3) and the eccentricities (e1, e2, e3 from top to bottom). Blue lines are for the stable symmetric ACR families (S6) while the black lines for the stable asymmetric family (A6). On each panel, the lowest line is for the planets in the Kepler-51 system, while the middle and the top lines are for the cases of 10 and 100 times more massive planets, respectively (see text). The critical points, at which the S6 family becomes unstable and the asymmetric solution appears, are denoted by red crosses.

4 Dynamics around stable ACRs

For the 1:2:3 resonant chain, the S₆−A₆ is the only stable family that could evolve from circular orbits to regime of high eccentricity. It is possible that a planetary system may evolve along a “path” guided by the S₆−A₆ family. In this part, by constructing dynamical maps on the appropriate representative planes, we explore the motion in the vicinity of this evolution path along the S₆−A₆ ACRs.

4.1 Initial conditions

All the ACR solutions were obtained in the averaged Hamiltonian system, but to construct the dynamical maps we have to use the osculating orbital elements to numerically integrate the orbits. Therefore, we should choose the initial conditions carefully to build a proper association between the averaged system and the real system. Generally, the osculating elements are related to averaged elements by the canonical transformation (see e.g. Ramos et al. 2015; Charalambous et al. 2018, for two-planet and three-planet cases, respectively). But in this paper, we use a simple transformation based on the orbital configuration, particularly the symmetric orbits, as follows.

Typically, the period ratios (n₁/n₂, n₂/n₃) of planets inside a resonant chain have only very small oscillation around their nominal value (n_1c/n_2c, n_2c/n_3c), otherwise the period commensurability and resonance between planets would break. Therefore, the period ratios, thus the osculating semi-major axes a_i, can be fixed at the values of the corresponding stable ACRs ā_i, that is, a_i = ā_i (i = 1, 2, 3). The same treatment can be applied to eccentricities and resonant angles, e_i = ē_i (i = 1, 2, 3), ${\theta }_{i,i+1}={\overline{\theta }}_{i,i+1},\mathrm{\Delta }{\varpi }_{i,i+1}={\overline{\mathrm{\Delta }\varpi }}_{i,i+1}$ (i = 1, 2). In fact, we construct dynamical maps on the (e₁, e₂) plane with fixed constants L and G defined in Eq. (10). For each point on the (e₁, e₂) plane, the conversation of L is automatically satified when the semi-major axes a₁, a₂, a₃ were given (in following maps, they were fixed to the values of the exact ACRs), while the eccentricity of the outer planet e₃ can be determined through the total angular momentum G. After a₁, a₂, a₃, e₁, e₂, e₃, θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃ have been given as the initial conditions, only one angular variable is left free. We choose the synodic angle Q as the last variable to be determined.

Let $(X_i,Y_i,{\dot{X}}_i,{\dot{Y}}_i)$ be the position and velocity of the i-th body m_i in the rotational frame. An orbit starting from the GX axis at time t₀ = 0 and satisfying X_i(t₀) ≠ 0, ${\dot{X}}_i\left(t_0\right)$ = 0, Y_i(t₀) = 0, ${\dot{Y}}_i\left(t_0\right)$ ≠ 0 (i = 1, 2, 3) is symmetric with respect to the GX axis, because the symmetric condition Eq. (6) is apparently fulfilled. The algebraic manipulation reveals that the symmetry in the rotational frame leads to $$\begin{array}{rl}& |{\mathit{r}}_i(t)-{\mathit{r}}_j(t)|=|{\mathit{r}}_i(-t)-{\mathit{r}}_j(-t)|,\\ & {\mathit{p}}_i(-t)\cdot {\mathit{p}}_j(-t)={\mathit{p}}_i(t)\cdot {\mathit{p}}_j(t).\end{array}$$(19)

Thus, the Hamiltonian $\mathcal{H}$ should follow the same symmetry $\mathcal{H}$(−t) = $\mathcal{H}$(t). Since $\mathcal{H}$ depends explicitly on rather Q than t, while Q is a linear function of t, up to the 1st order of the planet-star mass ratio (m₁ + m₂ + m₃)/M, we have $\mathcal{H}$(−Q) = $\mathcal{H}$(Q). Therefore, the Hamiltonian must be symmetric with respect to Q₀ = Q(t = 0) and $\mathcal{H}$(Q₀) must be an extremum. To determine the proper osculating orbital elements, we need to find the appropriate Q, and now it is equivalent to find the Q at which $\mathcal{H}$(Q) reaches the extrema.

We plot in Fig. 8 the variation of the Hamiltonian with respect to the synodic angle Q for both symmetric and asymmetric resonant configurations. It is clear that the Hamiltonian is symmetric with respect to the red vertical line in Fig. 8a, which corresponds to the global minimum. The above calculations can be applied to all the symmetric ACRs, and all the angles (ϖ_i, M_i) (i = 1, 2, 3) at the symmetric solutions have been summarised in Table 1.

For asymmetric configurations, as shown in Fig. 8b, the symmetry condition cannot be fulfilled. The Q value, at which the Hamiltonian attains the minimum, is calculated and adopted as the initial configuration.

Fig. 8 Variation of the Hamiltonian in two synodic periods $2Q_{12}^T=12\pi$ for (a) the symmetric configuration and (b) the asymmetric configuration. The vertical line denotes the position when the Hamiltonian attains the minimum.

4.2 Max (Δ e) indicator and χ² criterion

For each point on the (e₁, e₂) plane, we numerically integrate the equations of motion up to 10⁵ yr from the corresponding initial conditions (osculating elements). We note that the “year” here is defined as the orbital period of the middle planet m₂. During the integration, we keep tracking the maximal and minimal values attained by the planetary eccentricities and adopt the variation range (Δe₁, Δe₂, Δe₃) as the indicator to characterise the orbital stability. As mentioned in Ramos et al. (2015), although the Δe is not a rigorous measure of chaos of the motion, it’s a useful tool for the estimation of both the possible fixed points and the dynamical separatrix of the resonant region. In practise, we may consider the orbits with small Δe as “stable” and the ones with large Δe as “chaotic” (close to the dynamical separatrix).

In a two-body MMR, after the numerical averaging over the synodic angle Q, the system is reduced to a Hamiltonian of only 2 DoF. The existence of the resonant invariant torus, or dynamical separatrix (Alves et al. 2016), forbids the transition between the resonant and the non-resonant region. Inside the resonant zone, the resonant angles librate around the ACR solutions, while they circulate outside the resonant region. Thus, different types of motion can be easily distinguished from each other. Unfortunately, in the three-planet resonant system, the resonant angles behave in a much more complicated way. Due to the perturbation from the third planet, a two-body resonant angle may librate irregularly with a large amplitude and appear like in a complex of both libration and circulation. As a result, it’s difficult to determine accurately whether the system is in a resonance (or dominated by the resonance). Following the idea by Gallardo (2014), another indicator, χ² of the resonant angles, is used in this paper.

The indicator χ² is defined to characterise different effects of the resonance in a three-planet system: deep in resonance (resonant angles librating), at the boundary of resonance (resonant angles circulating or librating with large amplitudes), and outside resonance.

The basic idea is to analyse the resonant angle over a period of time T*. The statistical difference χ² between the distributions of a resonant angle and a uniform distribution in [−π, π] is adopted as an indicator to determine whether the system is in a resonance: the larger the χ², the more likely the resonance angle is in libration, and vice versa.

The semi-quantitative criterion will be given through a standard pendulum model with Hamiltonian: $${\mathcal{H}}_{\text{pend }}=\frac{p_{\phi }^2}{2}-\cos \phi .$$(20)

From each point of a 100 × 100 grid on the phase space [p_φ, φ] ∈ [−4, 4] × [−π, π], we integrate the corresponding Hamiltonian equations and output N = 2¹³ = 8192 points ${\left\{{\phi }_i\right\}}_{i=1}^{N=8192}$ evenly in a time step Δt = 5. The interval [−π, π] is divided into 180 equal bins, and the χ² of such an orbit is computed as $${\chi }^2=\sum _{k=1}^{180}\frac{{(O_k-E_k)}^2}{E_k},$$(21)

where O_k is the number of φ_i observed in the k-th bin, and E_k is the expected number in the same bin if the 8192 values are in a uniform distribution. We show in Fig. 9 the phase space of the standard pendulum and the χ² map.

Apparently, all the points located inside the libration region in Fig. 9 have χ² > 1 (log χ² > 0), while those points with χ² < 0.1 are surely in the circulation region. For initial conditions in the close vicinity of the separatrix, their orbital periods are so long that the total integration time could not cover several times of their periods and therefore, some orbits in the circulation region might be misunderstood as “libration”. Empirically, the resonant behaviours can be estimated by the value of χ²:

• χ² > 1, the resonant angle is in libration and the system is deeply inside the resonance;

• χ² < 0.1, the resonant angle is in circulation and the system is outside the resonance;

• 0.1 < χ² < 1, the system is affected by the resonance, but the resonant angles may experience both the libration and circulation in a long time span. We classify this behaviour as the motion near the separatrix.

For the MMR between planets m_i, m_j, the resonance is characterised by the libration of at least one of the resonant angles $${\theta }_{ij}^l=k_{ji}{\lambda }_j-k_{ij}{\lambda }_i-l{\varpi }_i-(q_{ij}-l){\varpi }_j,l=0,1,\cdots ,q_{ij},$$(22)

where q_ij = |k_ij − k_ji| is the order of the resonance, and the librating angle is referred as the true resonant angle (Alves et al. 2016). Therefore, the criterion χ² of an MMR is defined as the largest χ² value of all resonant angles, and we may adopt the same empirical criterion obtained in the pendulum model.

Fig. 9 χ2 map for the standard pendulum model, (a) Phase structure (energy level) on (φ, pφ) plane. The broad red curve marks the separatrix between the libration and circulation region, (b) χ2 value in logarithm indicated by colour.

Fig. 10 Topology (contour line) and dynamical map (colour map) on the (e1, e2) plane for Case A, with constants L and G determined by the initial condition on the stable ACR solution (e1, e2, e3, θ12, Δϖ12, θ23, Δϖ23) = (0.0124, 0.0237, 0.0128, 0, π, 0, π), which is indicated by a red dot on all maps. The eccentricity variations Δei are used as the stability indicator as labelled in each panel.

4.3 Motion around symmetric ACRs at low eccentricity

Along the S₆−S₆ family, the symmetric ACR solution is stable at low eccentricity (Fig. 6). Arbitrarily, we select on the ACR one initial condition (e₁ e₂, e₃) = (0.0124, 0.0237, 0.0129), which is labelled by ‘A’ in Fig. 6. The resonant angles of this solution are (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (0, π, 0, π), and the period ratios are set as the nominal values in the 1:2:3 resonant chain. Keeping the constant L and G calculated from the above a_i, e_i(i = 1, 2, 3), we generate a set of initial conditions both for the averaged Hamiltonian system and for the corresponding osculating system. We denote the case with these initial conditions by “Case A”, and investigate the motion.

We present in Fig. 10 on the (e₁, e₂) plane both the topologies (contour curves) of the corresponding averaged Hamiltonian and the dynamical maps (coloured map) obtained from numerical integration of the osculating orbital elements. The ACR solution appears as the local maximum and is surrounded by ellipse-like energy curves, implying its stability. The maximum of the eccentricity variation Δe_i is adopted as the indicator of stability as labelled on each panel in Fig. 10. And particularly, the maximum of (Δe_i, Δe₂, Δe₃) is used in the last panel.

As clearly shown in Fig. 10, all the four maps reveal that the motion is very stable around the ACR with very small eccentricity variations. Further analyses show that these stable orbits are characterised by the libration around the stable ACR, performing the so called quasi-periodic motion. As a matter of fact, such a regular region can occupy a large area on the (e₁, e₂) plane.

Away from the stable ACR, the dynamical maps of different indicators may show different patterns. The region of small Δe₁ and Δe₂ extends roughly along either the e₁ (horizontally) or the e₂ direction (vertically). At about (e₁, e₂) = (0.025, 0.02) in Fig. 10a, the planet m₁ would suffer relatively large eccentricity excitation (Δe₁ > 0.02) and the corresponding orbit is much less regular. In addition, starting from about (e₁, e₂) = (0.032, 0.006) we also notice a very narrow strip with relative large value of Δe₁.

The dynamical features in Figs. 10b,c for Δe₂ and Δe₃ are similar. At the bottom of the (e₁, e₂) plane, Δe₂ attains its maximum value (~0.04), implying (relatively) the least stability of motion there. And in all three panels a, b and c, deviated from the ACR, regular region can also be found at the right corner of (e₁, e₂) plane, where motions are distinct from the quasi-libration around the ACR.

Finally, we find that all these dynamical features can be revealed by the maximum of the three eccentricity variations, i.e. max(Δe)= max(Δe₁, Δe₂, Δe₃), as shown in Fig. 10d. Particularly, we see in Fig. 10d that the most stable motions occur in three regions, around the ACR, around (e₁, e₂) = (0.032, 0.017), and around (0.04, 0.01), respectively. Below in this paper, we adopt max(Δe) as the stability indicator and construct the dynamical maps.

The χ² maps might provide more details on the resonant dynamics. In Fig. 11, we depict the χ² maps for all the MMRs inside the 1:2:3 resonant chain. For all MMRs, the closer to the stable ACR, the larger the χ², implying that in the close vicinity of the stable ACR solution, all the resonant angles are librating, as well as the three-body Laplace resonant angle. The regular motion, characterised by the quasi-periodic libration around the ACR, is common around the ACR solution.

The χ² maps in Fig. 11 reveal that, as e₁ increases (while e₂ remains the value of ACR), the 2:3 two-body resonance (MMR23) in the outer planet pair persists and seems to play an essential role to the resonant stability, while other resonances in the resonant chain are weaken or even destroyed. The region with small Δe₂ in Fig. 10b coincides with the region of χ² > 1 in Fig. 11b. When e₂ deviates from the ACR value (with e₁ almost being frozen), the inner planet pair’s 1:2 resonance (MMR 12), the 1:3 resonance (MMR 13) between m₁ and m₃, and the three-planet Laplace resonance, persist in a range of e₂, as shown by the vertical strip-like structures of large χ² in Figs. 11a, c, d. Therefore, away from the ACR solution along the direction of e₁, although the MMR23 is stable, the strength of MMR12, MMR13, and Laplace resonance becomes weaker, implying that the corresponding resonant angles may librate with large amplitudes and eventually the Laplace angle φ₀ will be in circulation.

A nearly vertical narrow strip of small χ² that divides the stable region into two parts, namely Region 1 and Region 2, can be seen in Figs. 11a, c, the separatrix-like structure corresponds to the broad region with large Δe₁ on the dynamical map (Fig. 10a). We check carefully the motion in these two regions and find that the orbits in Region 1 are librating around the stable ACR (Fig. 12), while in Region 2, although all the two-body MMRs (i.e. MMR 12, MMR23, and MMR 13) occur, the Laplace angle φ₀ circulates (Fig. 13), as indicated by the small χ² value in the area of Region 2 in Fig. 11d.

A typical example of orbits in Region 1 is displayed in Fig. 12. The mean motion ratios n₁/n₂, n₂/n₃ librate with very small amplitudes around the nominal values (2 and 1.5), implying that the semi-major axes are almost constants. The eccentricities e₁, e₂, e₃ oscillate around the periodic solution (symmetric ACR, indicated by horizontal lines in the corresponding panels). And the tiny libration of critical resonant angles (including the Laplace angle) around 0° or 180° reveals that the motion is deeply in the resonant chain, even we know that the initial condition is not exactly on the ACR solution.

But for orbits in Region 2, as shown in Fig. 13, the resonant angles of MMR12 and MMR23, i.e. θ₁₂ and θ₂₃, librate around 0° with very small amplitude, as well as the apsidal difference Δϖ₂₃ librates around 180°. As a result, the resonant angle of MMR13, ${\theta }_{13}^1={\theta }_{12}+{\theta }_{23}+\mathrm{\Delta }{\varpi }_{23}=3{\lambda }_3-{\lambda }_1-{\varpi }_1-{\varpi }_3$ (see Eq. (22) for the definition) librates and the system is surely in the MMR between the non-consecutive planets m₁ and m₃. However, due to the circulation of Δϖ₁₂, the Laplace angle φ₀ = θ₁₂ − θ₂₃ − Δϖ₁₂ eventually circulates, indicating that the system is not in the three-body MMR. The eccentricities of the stable ACR solution are also plotted in Fig. 13, in which it can be seen that all e₁, e₂ and e₃ are far away from these ACR values. This implies that the typical orbit in Region 2 is no longer librating around the stable ACR. The motion in Region 2 in fact follows a new stable mode in the resonant chain, in which the dynamics is actually dominated by the individual two-body MMRs. Since in such a configuration each pair of planets is in the corresponding two-body MMR but the three-planet Laplace resonance does not happen, we refer to this type of configuration as “joint two-body resonance”.

It is worth noting that a separatrix-like structure indicated by the red arrow in Fig. 11a further splits Region 2 into two subregions. But the orbits from both sides of this structure are very similar to each other, characterised by libration of θ₁₂, θ₂₃, θ¹₁₃ and circulation of φ₀. The only difference between orbits from two sides is the oscillation centres of eccentricities, (e₁, e₂, e₃) ≈ (0.032, 0.017, 0.009) for the left side (as the one shown in Fig. 13) and (e₁ e₂, e₃) ≈ (0.040, 0.012, 0.007) for the right side. Most probably, such a fine structure in Region 2 is attributed to the nonlinear terms in the perturbation (Rath et al. 2022). We also note that the corresponding structure can also be found in the map of maximum of eccentricity variation (Fig. 10).

Finally, at the bottom of the dynamical map (Fig. 10 at e₂ ~ 0), we observe a thin region with a large eccentricity variation. The eccentricities (e₂, e₃) in this thin region are excited mainly because the outer two-body MMR (MMR23) is close to its boundary, as indicated by the small χ² value in the bottom area in Fig. 11b.

Fig. 11 Resonant structure indicated by χ2 on the (e1, e2) plane for Case A. Four panels are for the two-body MMRs between (a) m1 and m2, (b) m2 and m3, (c) m1 and m3, and (d) the three-planet Laplace resonance with resonant angle φ0. The large value of χ2 (light colour) indicates that the resonant angle is in libration, while the small χ2 (dark tones) indicates circulation. The narrow strip-like structure indicated by a black arrow in (a) divides the map into two parts, namely Region 1 and Region 2. Moreover, another narrow strip (indicated by the red arrow) can be seen in Region 2.

Fig. 12 Typical behaviour of the orbits in Region 1 (Fig. 11a). The temporal variations of n1/n2, n2/n3, e1, e2, e3 are in left panels, while right panels are for resonant angles θ12, Δϖ12, θ23, Δϖ23, φ0. The orbits are librating quasi-periodically around the stable symmetric ACR solution S6. The eccentricities of S6 are indicated by horizontal dashed lines.

Fig. 13 Same as Fig. 12 but for Region 2 in Fig. 11a.

4.4 Motion around asymmetric ACRs at moderate eccentricity

The symmetric ACR solution loses its stability at a specific point, and it bifurcates into two asymmetric ACR solutions. We select again arbitrarily an initial condition from the asymmetric ACR branch, and perform the similar calculations as for the symmetric case. The eccentricities at this initial point, indicated by the black solid circle labelled “B” in Fig. 6, are (e₁, e₂, e₃) = (0.0602, 0.0903, 0.0440). We call this case “Case B”.

The resonant angles at this initial point of asymmetric ACR are (θ₁₂, Δϖ₁₂, θ_23, Δϖ₂₃) = (26.7°,−119.5°, 1.16°, 179.5°). From these initial values, we obtain the corresponding constants L and G. For comparison we also calculate the unstable symmetric ACR solution with the same L and G, of which the eccentricities and resonant angles are (0.0569, 0.0846, 0.0531, 0, 180°, 0, 180°), i.e. the red solid circle in Fig. 6. Then, for both of these two initial conditions, we calculate the topology and dynamical maps, and present them in Fig. 14. In this figure, the max(Δe) is chosen as the stability indicator.

Comparing the two panels in Fig. 14, we find that the topology of the unstable symmetric ACR is similar to the one of the stable asymmetric ACR, but their dynamical maps show different patterns. The non-zero values of max(Δe) ~ 0.04 at the close proximity of the unstable symmetric ACR in Fig. 14a implies the orbits there deviate from the unstable stationary solution in the evolution. On the contrary, the stable asymmetric ACR in Fig. 14b is surrounded by regular orbits as indicated by the max(Δe) ~ 0.

A closer look at the orbits around the unstable symmetric ACR reveals that they are still in resonance, librating around both the asymmetric ACRs and the symmetric ACR with large amplitude, in a way similar to the horseshoe orbit around the Lagrange points L₄ − L₃ − L₅ in the circular restricted three-body problem (see e.g. Murray & Dermott 2000). One typical orbit in that region is shown in Fig. 15. Although it behaves chaotically due to the instability of symmetric ACR, the libration of the resonant angles indicate that the motion is still inside the resonant chain after the end of the numerical integration. Such horseshoe-like orbits can surely survive in the resonant chain in the long term.

On the other hand, the orbits near the stable asymmetric ACR solution are in a regular area (Fig. 14b), and they are protected by the resonant chain. We present in Fig. 16 the χ² map on the representative plane where all the initial resonant angles are fixed as the ones of the asymmetric ACR solution. As can be seen in Fig. 16, the close proximity around the stable ACR is deeply involved in the MMRs, as indicated by the large χ² value (yellow colour). Consistent with the dynamical map in Fig. 14b, such regular region occupies large phase space, with the eccentricities e₁, e₂ extending to e₁ → 0.1, e₂ → 0.1. For orbits inside the yellow region (χ² > 10) in Fig. 16, the libration amplitudes of resonant angles θ₁₂, θ₂₃ are smaller than 10°. As a consequence, some three-body resonant angles Ψ_[M,N] = Mθ₁₂ + Nθ₂₃ librate. We note that the libration of Ψ in this case does not indicate a “real” three-planet resonance but is the necessary consequence of the libration of two-body resonance angles.

Similar to the resonant structure near the stable symmetric ACR (Fig. 11), the librating region of MMR23 extends in the e₁ (but not e₂) direction while the one of MMR12 extends along the e₂ direction. The MMR between planets m₁, m₁₃ (MMR13) and the Laplace resonance appear as the intersection of MMR12 and MMR23. With increasing displacement from the ACR, we notice a broad stripe with small χ² in Fig. 16a that divides the libration zone into two parts, namely Region 1 and Region 2 again. This type of “chaotic border” is similar to the one in Fig. 11a for Case A. And it has apparent correspondence in the dynamical map in Fig. 14b for Case B.

Judging from the χ² value, the orbits in Region 1 are involved in all three two-body MMRs and the three-body Laplace resonance. The motion is nearly the same as the one in Case A shown in Fig. 12, only except that the resonance centre is at the asymmetric values. Close to the ACR solution, the motion is regular, and departing away from the ACR the irregularity increases, as shown by the max(Δe) value in Fig. 14b. However, in Region 2, which is far away from the ACR solution, the small max(Δe) value indicates that the motion is regular, although the χ² value implies that some resonances have broken. To show the orbital behaviour in Region 2, we illustrate a typical example in Fig. 17.

As shown in Fig. 17, the period ratios n₁/n₂, n₂/n₃ still stay close to the nominal values (2,1.5). The eccentricities oscillate regularly around the centre $(e_1^{\ast },e_2^{\ast },e_3^{\ast })=(0.19,0.072,0.035)$ that is quite far away from the ACR solution. Instead of librating around the ACRs (the unstable symmetric one and the two stable asymmetric ones), the resonant angles behave in a similar way as the ones in Case A illustrated in Fig. 13, that is, all the two-body MMRs occur (at least one of the corresponding resonant angles librates) while the three-body Laplace resonance does not (φ₀ circulates). Therefore, the “joint two-body resonance” happens in this region of (relatively) high eccentricity.

The orbits of “joint two-body resonance” found in both Cases A and B (Figs. 13 and 17) reveal that this configuration can occur in a wide range of initial conditions. We note that this type of mode of stable “joint two-body resonance” is missed in our averaged model, i.e. the motion is not the resonant periodic orbits, but it is a common outcome from planetary migration and tidal dissipation, as shown in Morrison et al. (2020).

Fig. 14 Same as Fig. 10 but for Case B. The constants L and G are determined by the initial condition on the stable asymmetric ACR solution (black solid circle in Fig. 6). (a) Around the unstable symmetric ACR solution with the same L, G. The corresponding eccentricities e1, e2, e3 are 0.0569, 0.0846, 0.0531 (red solid circle in Fig. 6). (b) Around the stable asymmetric ACR solution. The colour in both panels indicates the value of max(Δe).

Fig. 15 Horseshoe-like orbit near the unstable symmetric ACR. The unstable symmetric ACR and the stable asymmetric ACRs are indicated by the red and black horizontal lines, respectively.

Fig. 16 Same as Fig. 11 but for Case B. The map is divided into two parts (Region 1 and Region 2) by a stripe-like structure of relatively smaller χ2 value, as indicated by the arrow in panel a.

Fig. 17 Same as Fig. 13 but the initial conditions are from the high eccentricity region of Fig. 14 (see text).

4.5 Stable motion inside the 1:2:3 resonant chain

Basically, in the close neighbourhood of the ACR solutions as we have investigated in this paper, two types of stable motion of the 1:2:3 resonant chain can be found. One is the stable ACR (resonant periodic orbits) and the associated quasi-periodic motion, where all the resonant angles, including the two-body MMRs and the Laplace resonance angle, are librating (χ² > 1). The other one is the “joint two-body MMR”, where θ₁₂, θ₂₃ and ${\theta }_{13}^1$ are librating while the Laplace angle φ₀ circulates. Judging from the dynamical maps and our long-term numerical simulations, the difference in stability between these two types of motion is negligible. Therefore, in the close vicinity of the nominal resonance, we might conclude that the individual two-body MMRs dominate the stability of the 1:2:3 resonant chain, while the libration of the three-body Laplace angle is not a necessary condition of the resonant chain but a necessary outcome of the two-body resonant angles’ librating.

This feature can be explained through the Hamiltonian model. The perturbed part $\mathcal{H}$₁ of Hamiltonian in Eq. (7) only contains the interactions between each two planets. After the numerical averaging process (eliminating all the non-resonant term), up to the 1st order in planet-star mass ratio, the final averaged Hamiltonian $\overline{\mathcal{H}}$ in Eq. (13) only contains the terms induced by the two-body MMRs but not the three-body resonance. Consequently, for a planetary system close to a resonant chain, the dynamics is mainly dominated by the two-body MMRs.

Therefore, if we can verify that a planetary system locates inside all the two-body MMRs (at least one critical resonant angle librates for each MMR), no matter how the Laplace angle φ₀ behaves, we will know for sure that the resonant chain is formed. In the practice of dynamical studies, the χ² indicator can distinguish libration from circulation and thus provide the priori estimation of the stability in short-time numerical integration (only up to tens of resonant periods, empirically 10⁴−10⁵ synodic periods), which is in fact easy to perform.

In addition, the three-body Laplace resonant angle φ₀ can be easily detected in transit data. Once the orbital periods and semi-major axes of the planets are determined, the behaviour of φ₀ is often adopted as an indicator to probe whether the planetary system is involved in a resonant chain (e.g. Jontof-Hutter et al. 2016). According to our analyses of the 1:2:3 resonant chain, the libration of the three-body resonant angle can almost ensure that the system is librating around a certain ACR. In this case, the dynamics and the possible behaviours of the orbital elements can be explained and predicted by studying the Hamiltonian topology and the dynamical maps on appropriate representative planes.

If a system is not in the close vicinity of the nominal resonances, i.e. if the offset from the exact period commensurability is considerable, the pure three-body MMRs (with absence of any two-body MMR) may play its role. Recently, Rath et al. (2022) propose an analytical model (perturbed pendulum) of the chaotic dynamics near the resonant chain, and they demonstrate that the pure three-body resonances should appear as the secondary resonance with the critical angle defined as Ψ_[M,N] = M_ϕ + Nψ, where ϕ and ψ are the critical angles of individual two-body MMRs: $$\varphi =k_{21}{\lambda }_2-k_{12}{\lambda }_1,\psi =k_{23}{\lambda }_2-k_{32}{\lambda }_3.$$(23)

It is important to note that the ϖ terms have been neglected. Obviously, the [M,N]=[1,1] type is the largest secondary resonance (for the example of the 1st order three-body resonance, see Cerioni & Beaugé 2023). For the 1:2:3 resonant chain, it’s just the Laplace resonance φ₀ = λ₁ − 4λ₂ + 3λ₃. Specifically for the 1:2:3 resonant chain, Antoniadou & Voyatzis (2022) find the stable region associated with this zeroth-order three-body resonance. However, we did not detect these configuration, perhaps because our initial condition are always chosen at exact nominal resonance. More efforts are needed to investigate how the pure three-body resonance affects the resonant configuration and the long term stability.

5 Resonant configuration via convergent migration

So far, we have shown both by the Hamiltonian method and by numerical simulations what the orbital configurations might be in a planetary resonant chain. Generally, an MMR between two planets can be easily attained via planetary convergent migration, which can also lead multiple planets to form a resonant chain and then to achieve an orbital configuration of high eccentricity (see e.g. Voyatzis 2016). In this section, we study whether the orbital configurations we have analysed above can be attained through convergent migration of planets.

Due to the interaction between the planet and the circumstellar disk (gaseous disk or planetesimal disk), the migration of planet is quite common. For simplicity, in this paper the migration is simulated by adding an artificial dissipative force on the planet. Specifically, the torque acting on the i-th planet is assumed as (see e.g. Delisle et al. 2015): $${\ddot{\mathit{r}}}_i=-[\frac{{\mathit{v}}_i}{2{\tau }_{a,i}}+\frac{2({\mathit{v}}_i\cdot {\mathit{r}}_i){\mathit{r}}_i}{r_i^2{\tau }_{e,i}}],$$(24)

where r_i, υ_i are the position and velocity vector, and τ_a,i,τ_e,i are timescales of orbital migration and eccentricity damping, for the ith planet. Up to the first order in eccentricity, the effects of the dissipative force on the planetary semi-major axis and eccentricity read: $$a_i(t)=a_i(t=0)\exp (-\frac{t}{{\tau }_{a,i}}),e_i(t)=e_i(t=0)\exp (-\frac{t}{{\tau }_{e,i}}).$$(25)

Constant timescales were adopted in our simulations, and in order to fulfil the global convergent condition (e.g. Beaugé & Cerioni 2022; Wong & Lee 2024), different migration parameters will be chosen for different resonant chains.

5.1 1:2:3 resonant chain

To reach the 1:2:3 resonant chain, following Beaugé & Cerioni (2022), we assumed that all three planets migrated inwards and the migration timescales τ_a,i were set to be 0.5, 0.2, 0.1 Gyr for planets m₁, m₂ and m₃. The eccentricity damping timescales τ_e,i were set in such a way that for each planet κ = τ_a,i/τ_e,i = 50. As we did for the dynamical maps, the initial semi-major axis of m₂ is set to be a₂ = 1 AU. All planets start from coplanar and circular orbits and the initial period ratios T₂/T₁ and T₃/T₂ are set in a 15 × 15 grid in the ranges of [2.01, 2.02] and [1.51, 1.53]. The initial angles ϖ_i, λ_i were all set to be 0°. This choice corresponds to a collinear configuration, in which the mutual distances between planets are minimum. All the orbits starting from the 225 initial conditions were then integrated for a total time span of T = 10⁷ yr.

During the numerical integrations, we keep tracking the behaviour of the resonant offset Δ_ij = n_i/n_j − k_j/k_i between the adjacent planets. Once all the resonant offsets are smaller than a given level ${\mathrm{\Delta }}_{ij}^c$, we then monitor the behaviour of resonant angles by checking their χ² and classify the motion. It should be noted that the critical ${\mathrm{\Delta }}_{ij}^c$ can be chosen empirically, and after some tests, we set ${\mathrm{\Delta }}_{ij}^c$ = 20 × (m₁ + m₂ + m₃)/M × k_j/k_i.

For the 1:2:3 resonant chain, we found that only two types of resonant configuration are achievable in our numerical simulations. The first one is that the planetary system attains the resonant chain and eventually librates around the stable ACR solutions (S₆, A₆). One such example is shown in Fig. 18. As the planets migrate, both period ratios n₁/n₂ and n₂/n₃ decrease, while the eccentricities remain close to zero. After the system reaches the 1:2:3 resonant chain at t ≈ 10⁶ yr, the period ratios are subsequently locked into the resonance and the eccentricities are excited from then on. Meanwhile, all the angles begin to librate around the values of the symmetric ACR, (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (0, π, 0, π), i.e. the S₆ family. Once the bifurcation point is reached at t ≈ 2.4 × 10⁶ yr, the system follows the asymmetric family A₆ and finally librates around the asymmetric ACR at about $(e_1^{\ast },e_2^{\ast },e_3^{\ast },{\theta }_{12}^{\ast },\mathrm{\Delta }{\varpi }_{12}^{\ast },{\theta }_{23}^{\ast },\mathrm{\Delta }{\varpi }_{23}^{\ast })$ = (0.0371, 0.0623, 0.0345, 24.6°, −133.5°, 0.5°, 178.9°).

The second one is the “joint two-body MMR”, of which an example is shown in Fig. 19. Similar patterns as in Fig. 18 can be observed for (n₁/n₂, n₂/n₃, e₁, e₂, e₃): once the system encounters the 1:2:3 resonant chain, the eccentricities monotonously increase while the period ratios remain at the nominal values. However, the resonant angles seem not be affected by the ACRs of the resonant chain, and the three-body Laplace resonant angle (φ₀) always circulates but not librates. All these regular dynamics inside the 1:2:3 resonant chain are just as predicted by our previous analyses.

Comparing the orbital evolution in Figs. 18 and 19, we find that in the latter case the period ratio n₂/n₃ has been trapped to 1.5 (3:2) at t₁ ≈ 1.2 × 10⁶ yr. At this moment, n₁/n₂ has not reached the nominal value of 2.0, but a change of the relative migration rate can be seen clearly in the top left panel of Fig. 19. The 1:2 commensurability between the inner two planets however seems to occur later at t₂ ≈ 1.8 × 10⁶ yr. Its influence on the resonant chain is reflected by the small “spikes” in the e₂ and e₃ variation in two bottom left panels.

The evolution in Fig. 19 seemingly implies that the dynamical configuration a resonant chain finally takes might be influenced by the forming sequence of individual two-body MMRs. But in fact, which configuration the resonant chain will occupy depends on the masses of planets, the initial conditions and the migration rates in a very complicated way. We check carefully all the 225 planetary systems in our simulations, and find that totally 180 out of 225 systems are trapped in the motion around the ACR as in Fig. 18 and the rest 45 systems are finally trapped in the joint two-body MMR. Unfortunately, we were unable to find out a rule regarding which type of motion a system will finally take in the resonant chain, libration around ACR or joint two-body MMR.

Aiming to check whether the 1:2:3 resonant chain can be formed through planets’ migration and whether the periodic solutions can also be obtained in this way, we have set the initial period ratios very close to the target resonant chain in our simulations. We performed some extra simulations from initial period ratios far from the resonant chain (n₁/n₂ > 2.1, n₂/n₃ > 1.6), and found that the capture into the resonance chain can also succeed, and all the captured systems are in either stable ACR or the joint two-body MMRs. The eccentricities may also affect the capture and the subsequent evolution. We note that the eccentricities will be inhibited by the eccentricity damping effect in the numerical model. However, for cases where the system encounters the 1:2:3 resonant chain in eccentric orbits (typically e ∈ [0.05, 0.1]), the resonant chain can also be achieved but the system may temporarily stay around the resonant chain with relatively large libration amplitudes. Such a large eccentricity variation may drive the system out of the resonant chain and eventually leaves only individual two-body MMRs. Therefore, in the following investigations, we only focus on the capture starting near the resonant chain and from near circular orbits.

Fig. 18 Formation and evolution of the 1:2:3 resonant chain via convergent planetary migration. The left panels show the temporal evolutions of period ratios and eccentricities, and the critical angles are shown in the right panels. The planetary system is captured in symmetric configuration S6 firstly, subsequently evolves along the A6 family and finally librates around the asymmetric configuration.

Fig. 19 Same as Fig. 18 but for a ‘joint two-body MMR’ case (see text). In the final resonance configuration, Δϖ12 and φ0 circulate while θ12, θ23, and Δϖ23 librate. The vertical lines in the left panels indicate two moments t1 = 1.2 Myr and t2 = 1.8 Myr.

5.2 Other resonant chains

So far, our investigations have been on the 1:2:3 resonant chain. The consistence between the averaged Hamiltonian and the numerical simulations motivates us to extend our semi-analytical model to other resonant chains. As far as we know, previous works, for example, Morrison et al. (2020) and Kajtazi et al. (2023), mainly focus on the 1st order resonant chains with (k₂₁ − k₁₂ = k₃₂ − k₂₃ = 1). To show how our methods may work for general resonant chains, below we present our analyses on some other resonant chains. We also note that starting from circular orbits, the ACR solutions of a resonant chain may emerge as symmetric configurations with critical resonant angles being 0, π, or as asymmetric ones.

5.2.1 Cases starting from symmetric configuration

We will show two cases, a 2:3:5 resonant chain and a 3:5:7 resonant chain. The former is a resonant chain of mixed orders consisting of the 1st and 2nd order MMRs between adjacent planets, while the latter is a resonant chain with two second order MMRs. Although no exoplanetary system has been confirmed in these resonant configuration, we present these two examples to demonstrate the availability of methods introduced in this paper for higher-order resonances. Meanwhile the existence of such configuration but absence of such real systems also implies that the stability of such configuration is weak.

For the 2:3:5 resonant chain, we consider the fictitious planetary system with central star of mass M = 1 and three planets m₁ = m₂ = m₃ = m = 10⁻⁶, and as before, the semi-major axis of m₂ is set as a₂ = 1 AU. Using the same technique we have applied to the 1:2:3 resonance chain, we obtain the symmetric and asymmetric ACR solutions and their stabilities for the 2:3:5 resonant chain. Following the definition of critical angles in Eq. (9), we define $$\{\begin{array}{ll}{\theta }_{12}=3{\lambda }_2-2{\lambda }_1-{\varpi }_1,& \mathrm{\Delta }{\varpi }_{12}={\varpi }_1-{\varpi }_2,\\ {\theta }_{23}=5{\lambda }_3-3{\lambda }_2-2{\varpi }_2,& \mathrm{\Delta }{\varpi }_{23}={\varpi }_2-{\varpi }_3,\end{array}$$(26)

and we find that the symmetric family with (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (0, π, π, π) is the unique stable ACR arising from the circular configuration. Along this family, the planets may evolve to the configuration of high eccentricities when the angular momentum G decreases. This stable symmetric family finally halts at (e₁, e₂, e₃) = (0.311, 0.275, 0.294) where the planetary system may experience close encounters. Our detailed analyses on the dynamical maps reveal that all the regular motions we found in this 2:3:5 resonant chain are characterised by quasi-periodic libration around the ACR.

To investigate the formation and evolution of this resonant chain, we perform numerical simulations of the convergent migration and resonance capture again. The initial period ratio T₂/T₁ is uniformly distributed in the range [1.51, 1.53], and T₃/T₂ in [1.68, 1.71]. All planets’ initial orbits are assumed to be coplanar and circular, and the angles ϖ_i = 0, λ_i = 0. The migration timescales for three planets τ_a,1, τ_a,2, τ_a,3 are set as 0.3, 0.2, 0.1 Gyr, respectively, and the eccentricity damping time scale τ_e,i is set to satisfy κ = τ_a,i/τ_e,i = 100.

We find that all of 225 fictitious systems are captured into the 2:3:5 resonant chain and evolve along the stable symmetric family (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (0, π, π, π). Instead of plotting the evolutions of orbital elements and the resonant angles, we depict on eccentricities planes the migrating path along this ACR family in Fig. 20. It’s clear that this stable ACR family guides the migration of the planetary system. The evolution eventually stops at $(e_1^{\ast },e_2^{\ast },e_3^{\ast })=(0.0460,0.0612,0.0122)$ in this example, but surely the final status depends on the migrating parameters. If all the torques on planets disappear at this moment, the system will oscillate around the stable ACR at $(e_1^{\ast },e_2^{\ast },e_3^{\ast })$ in the long term.

For the case of the 3:5:7 resonant chain, the critical angles are defined as in Eq. (9) $$\{\begin{array}{ll}{\theta }_{12}=5{\lambda }_2-3{\lambda }_1-2{\varpi }_1,& \mathrm{\Delta }{\varpi }_{12}={\varpi }_1-{\varpi }_2,\\ {\theta }_{23}=7{\lambda }_3-5{\lambda }_2-2{\varpi }_2,& \mathrm{\Delta }{\varpi }_{23}={\varpi }_2-{\varpi }_3,\end{array}$$(27)

and the stable symmetric ACR of (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (π, π, π, π) was found to dominate the regular motion when the eccentricities are small. This symmetric ACR solution loses its stability at (e₁, e₂, e₃) = (0.196, 0.238, 0.0845) and it bifurcates into two branches of asymmetric families in a similar way as in the 1:2:3 resonant chain. Along the asymmetric family, the eccentricities might finally attain (e₁, e₂, e₃) = (0.267, 0.568, 0.617). Since both the 3:5 and 5:7 MMRs are relatively high order resonances, the resonant capture can be observed only if the migration is very slow. We show an example of the formation of 3:5:7 resonant chain in Fig. 21, where the migration timescales τ_a,1, τ_a,2, τ_a,3 for three planets are 5, 2, and 1 Gyr, respectively. In fact, if a faster migration rate is applied, most of the systems will skip the 3:5:7 resonance chain and fall into a more compact configuration. We note that such a phenomenon has been reported in Kajtazi et al. (2023).

To produce the 3:5:7 resonant chain as in Fig. 21, we started our simulations from initial (n₁/n₂, n₂/n₃) = (1.668, 1.402), which are very close to the nominal values in the resonant chain. And an eccentricity damping effect with κ = 100 was imposed in our simulations. The system evolves along the stable symmetric ACR (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (π, π, π, π) and the eccentricity damping effect finally halts the system around the configuration $(e_1^{\ast },e_2^{\ast },e_3^{\ast })=(0.0272,0.0627,0.0356)$.

Fig. 20 Resonant capture and evolution of the 2:3:5 resonant chain on the eccentricities plane. The stable (θ12, Δϖ12, θ23, Δϖ23) = (0, π, π, π) ACR is indicated by blue dashed lines while the grey dots show the migrating path of the system evolving to higher eccentricities.

Fig. 21 Same as Fig. 20 but for the 3:5:7 resonant chain. The blue dashed lines indicate the location of the stable ACR of (θ12, Δϖ12, θ23, Δϖ23) = (π, π, π, π).

5.2.2 Cases starting from asymmetric configuration

For all the aforementioned cases (1:2:3, 2:3:5, and 3:5:7 resonant chains), the stable ACR solutions arising from the circular configuration are always symmetric ones. As a result, the three-body Laplace angle φ₀ is always librating around the symmetric value, 0 or π. In this sense, if the φ₀ is found to librate around 0 or π, the system must be in such a resonant chain. However, this rule breaks for the 2:3:4 resonant chain. Siegel & Fabrycky (2021) find that the unreduced three-body critical angle 2_φ0 (they adopt 2_φ0 = 2λ₁ + 4λ₃ − 6λ₂ instead of φ₀) for the 2:3:4 resonant chain shift to asymmetric values ~ 162° and ~ 198°, instead of the symmetric value 2π.

Motivated by this interesting result, we apply our model to the 2:3:4 resonant chain with mass parameters M = 1, m₁ = m₂ = m₃ = 10⁻⁵, similar to the ones in Siegel & Fabrycky (2021). The resonant angles are $$\{\begin{array}{ll}{\theta }_{12}=3{\lambda }_2-2{\lambda }_1-{\varpi }_1,& \mathrm{\Delta }{\varpi }_{12}={\varpi }_1-{\varpi }_2,\\ {\theta }_{23}=4{\lambda }_3-3{\lambda }_2-{\varpi }_2,& \mathrm{\Delta }{\varpi }_{23}={\varpi }_2-{\varpi }_3.\end{array}$$(28)

Different from the 1:2:3 case, the final averaged Hamiltonian $\mathcal{H}$ for the 2:3:4 resonant chain is not 2π-periodically dependent on the angles (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃). Following the procedure in Appendix A, we find that the degeneracy of these four resonant angles are (2, 2, 2, 1), respectively. Instead of the symmetric ACRs, we found two stable asymmetric ACRs in the near-circular configuration. We present these ACRs by black lines in Fig. 22. For better view, only one of the two asymmetric families is plotted in Fig. 22, the other family simply takes the same values in terms of (e₁, e₂, e₃) and the symmetrically opposite values of the resonant angles.

Starting from circular orbits and initial period ratios T₂/T₁ uniformly in [1.51,1.53] and T₃/T₂ in [1.34,1.36], three planets m₁, m₂, m₃ in 225 fictitious systems with different initial period ratios are forced to migrate with migration timescales τ_a,i = 1, 0.5, 0.2 Gyr, while the eccentricity damping effect satisfies κ = τ_a,i/τ_e,i = 50. All these systems are eventually captured into either one of the two asymmetric ACR families. And a typical example of resonant capture and evolution is overlaid on the stable ACRs in Fig. 22.

As shown in Figs. 22c, d, once the system reaches the 2:3:4 resonant chain just after the migration begins when the orbits are still near circular, the resonant angles are immediately locked into the asymmetric value (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) = (27.0°, 144.8°, 9.9°, 162.7°) and the Laplace angle librates around 2_φ0 = 162.7° (Mod[360°]). The system evolves subsequently along the asymmetric ACR family and eventually librates around the ACR located at $(e_1^{\ast },e_2^{\ast },e_3^{\ast })=(0.0357,0.0512,0.0265)$ and $({\theta }_{12}^{\ast },\mathrm{\Delta }{\varpi }_{12}^{\ast },{\theta }_{23}^{\ast },\mathrm{\Delta }{\varpi }_{23}^{\ast })=({33.4}^{\circ },{138.4}^{\circ },{15.0}^{\circ },{160.0}^{\circ })$.

Therefore, for the 2:3:4 resonant chain, our numerically averaged Hamiltonian successfully predicts the librating centres of resonant angles, as well as the adiabatic evolution induced by the convergent migration. In fact, starting from circular orbits, asymmetric ACR solutions are also found in the 3:4:6 and the 12:15:20 resonant chains (Siegel & Fabrycky 2021), and our calculations detected three and four stable ACR families starting from the circular configuration, respectively. For reference, we list in Table 2 the basic properties of the ACR solutions in some low-order three-planet resonant chains.

From Table 2, we see the number of stable ACR solutions is equal to the highest degeneracy number of the corresponding resonant angles. This number, and the symmetry properties (symmetric or asymmetric) of these stable ACR families, might depend not only on the resonant orders, but also on the degeneracies of the involved two-body MMRs. A general model for this question is still needed.

Fig. 22 Resonant capture and evolution of the 2:3:4 resonant chain. The black lines indicate the location of the asymmetric ACR solution arising from the circular configuration. The grey dots stand for the capture and evolution of planets by the resonant chain in the numerically simulated convergent migration.

5.3 Remarks on final resonant configurations

As we have shown above, starting from circular configuration (e₁ = e₂ = e₃ = 0), the stable families of periodic orbits (ACR solutions, either symmetric or asymmetric) always exist in the resonant chains. These solutions guide the formation and adiabatic evolution of the resonant chains during the planetary migration. Depending on the migration parameters (τ_a,i, τ_e,i) and initial conditions, a triple-planet system finally attains different configurations of resonant chain, including the quasi-periodic motion around the ACR solutions and the joint two-body MMRs.

After the three planets have been trapped in a resonant chain, they may acquire high eccentricities as the system evolves along the path defined by the ACR solutions as the convergent migration continues. There are also some ACR solutions in the high eccentricity region that do not originate from the circular configuration (Voyatzis 2016). Since our methods of finding the ACR solutions work for arbitrary eccentricities, these high-eccentricity configurations can be calculated, but they cannot be attained through the adiabatic process from circular orbits as described above. Considering the fact that the typical eccentricities in the resonant chain are small (see e.g. Sun et al. 2017), we skip these configurations in this paper.

Another factor that can affect the final orbital configuration in a resonant chain is the planetary masses. For massive planets (e.g. Jupiter mass), the short-term perturbations will manifest their significances, and thus the formation and evolution of resonances in the system may not be smooth (adiabatic). One numerical example of the 1:2:4 resonant chain of massive planets can be found in Voyatzis (2016). The system can be captured into resonances through chaotic (irregular) processes and finally attain configurations of high eccentricities. Last but not least, in a compact planetary system, the tidal effect may modify the planets semimajor axes and eccentricities in an adiabatic way, thus it might play an important role in shaping the configuration of a resonant chain. As a matter of fact, the most well-known resonant chain, the literal ‘Laplace resonance’ in the three Jovian moons is believed to have formed through the tidal evolution (see e.g. Malhotra 1991).

6 Conclusion

In this paper, we generalise the semi-analytical model of MMR between two planets to the resonant chain in the triple-planet systems. By constructing and appropriately averaging the Hamiltonian that describes the dynamics of three planets in a resonant chain, we numerically computed the periodic orbits in the resonant chain. The periodic orbits in the resonance (also known as ACR solutions) are the stationary solutions to the averaged Hamiltonian equations and correspond to the extrema of the Hamiltonian function, which can be searched out numerically in the phase space.

Adopting the 1:2:3 resonant chain in the Kepler-51 system as an example, we calculated these ACR solutions and then analysed their stabilities by checking the topologies of the Hamiltonian on different representative planes. The stable symmetric ACRs we found are in good agreement with the stable symmetric periodic orbits obtained previously through the continuation method (Antoniadou & Voyatzis 2022). The symmetric ACR solutions lose their stabilities as the eccentricities increase, and they might bifurcate to asymmetric ACRs, which create stable configurations in the high-eccentricity region.

To investigate the orbital dynamics of the resonant chain near the stable ACRs, we constructed dynamical maps on the representative (e₁, e₂) plane and associate the orbital features with the Hamiltonian topology on the same plane. The dynamical maps show that the domain near the stable ACRs is always characterised by regular motion of quasi-periodic libration. The consistence between the topology and dynamical map simplifies the calculation of the stability of planetary systems in the resonant chain.

The behaviour of resonant angles in the resonant chain might be so complicated that we cannot always determine the resonance state simply by checking whether the corresponding critical angle is librating or not. In some cases, the short-term perturbations muddle the libration of resonance angles even when the system is dominated by the resonance. Thus, we calculated the deviation (χ²) of a resonant angle from a uniform distribution, and empirically defined the χ² criterion to distinguish the resonance.

In the vicinity of the ACRs, we found two types of stable motion in the resonant chain. One is the periodic or quasi-periodic motion at or around the ACR solution, which is characterised by concurrent libration of the two-body resonant angles and the three-body Laplace angle. The other is the quasi-periodic motion, in which only the two-body resonant angles but not the Laplace angle librate. We denote the latter motion as joint two-body MMRs. Another difference between these two types of motion is that the centre (in terms of eccentricities and critical angles) of libration for the former is just the ACR solution while for the latter the ACR is not even ‘reached by the libration. According to our experiences, in both cases, the two-body MMRs between planet pairs are the essential mechanism that stabilises the planetary resonant chain, and the pure three-body resonance contributes little to the stability. Last but not least, we note that the two-body MMRs in a resonant chain are not the same as the usual two-body resonance that is isolated from other planets in a planetary system.

Just like the usual two-body MMR, the resonant chain could also be formed via the convergent migration of planets. Setting the initial conditions and migration velocities for each planet carefully, we verified the formation and evolution of the resonant chain. Starting from near circular orbits, the planets can be trapped to form the resonant chain. And all these systems are either in the quasi-periodic orbits around the ACR solution or in the joint two-body MMRs, as predicted by our model. As the migration continues, the system evolves to a high-eccentricity region along the ACR solutions.

In the high-eccentricity region, there exist some ACR solutions that do not arise from near circular orbits. These solutions can also be determined using the methods introduced in this paper, but these resonant configurations generally cannot be achieved through the convergent migration. In addition, such a resonant chain of high eccentricity has a relatively poor stability.

The configuration and stability of a resonant chain may be affected by the planetary mass. Not only do the locations of the ACRs change, but also the stability region shrinks as the individual planet mass increases. In fact, the more massive the planets are, the stronger the disturbances between planets. Thus the resonant chains, especially those compact ones, are more likely to occur in planetary systems of sub-Neptunes or even lighter planets.

Our methods introduced in this paper can be applied to any resonant chains. We have briefly shown our calculations of the ACR solutions and formation for the 2:3:5, 3:5:7, 2:3:4, 3:4:6, and 12:15:20 resonant chains. Diversified orbital configurations are found in different systems. Particularly, both symmetric ACR solutions (e.g. in the 1:2:3, 2:3:5, and 3:5:7 resonant chains) and asymmetric ones (e.g. in 2:3:4, 2:3:6, 3:4:6, and 12:15:20) emerge from circular orbits. Our methods can also be applied to longer resonant chains consisting of more than three planets. As an example, in Appendix B, we calculated the periodic solutions of the four-planet 3:4:6:8 resonant chain in Kepler-223, and show the resonant configurations as the function of scaled AMD δ.

Acknowledgements

This work has been supported by the National Key R&D Program of China (2019YFA0706601) and National Natural Science Foundation of China (NSFC, Grants No.12373081, No.12150009 & No.11933001). We also thank the supports from the science research grant from the China Manned Space Project with No. CMS-CSST-2021-B08.

Appendix A Degeneracy of the resonance

The literal expansion of perturbation part $\mathcal{H}$₁ of the Hamiltonian in Eq. (7) takes the form $$\begin{array}{rl}{\mathcal{H}}_1& =\sum _{1\le i<j\le n}\frac{p_i{p}_j}{M}-\mathcal{G}\frac{m_i{m}_j}{r_{ij}}=\sum _{1\le i<j\le n}{H}_{1,ij}({\theta }_{ij},\mathrm{\Delta }{\varpi }_{ij},Q_{ij})\\ & =\sum _{1\le i<j\le n}\sum _{l,m,k}{H}_{1,ij}^{lmk}\cos (l{\theta }_{ij}+m\mathrm{\Delta }{\varpi }_{ij}+k{Q}_{ij}),\end{array}$$(A.1)

where l, m, k ∈ Z, θ_ij = k_jiλ_j − k_ijλ_i − (k_ji − k_ij)ϖ_i, Δϖ_ij = ϖ_i − ϖ_j, Q_ij = λ_i − λ_j. After the numerically averaging process, terms with k ≠ 0 in Eq. (A.1) will not appear in the averaged Hamiltonian $\overline{\mathcal{H}}$. Thus, the literal expansion of ${\overline{\mathcal{H}}}_1$₁ can be reduced to $$\begin{array}{rl}{\overline{\mathcal{H}}}_1& =\sum _{1\le i<j\le n}\frac{1}{T_{Q_{ij}}}{\int }_0^{T_{Q_{ij}}}(\frac{p_i{p}_j}{M}-\mathcal{G}\frac{m_i{m}_j}{r_{ij}})\mathrm{d}{Q}_{ij}\\ & =\sum _{1\le i<j\le n}{H}_{1,ij}({\theta }_{ij},\mathrm{\Delta }{\varpi }_{ij})\\ & =\sum _{1\le i<j\le n}\sum _{l,m}{H}_{1,ij}^{lm}\cos (l{\theta }_{ij}+m\mathrm{\Delta }{\varpi }_{ij}).\end{array}$$(A.2)

There are n(n − 1) resonant angles (θ_ij, Δϖ_ij), 1 ≤ i < j ≤ n appearing in the final expression at first glance. But in fact, these angles are linearly dependent. If we choose the two-body resonant angles between the adjacent planet pair θ_i,i+1, Δϖ_i,i+1 to characterise the configuration, the resonance angles between non-adjacent planets shall be linear combinations of angles θ_i,i+1, Δϖ_i,i+1 with the coefficients being fractions rather than integers. As a result, the Hamiltonian may not be 2π-periodic for all the resonant angles, e.g. the resonant configuration characterised by (θ₁₂, Δϖ₁₂, ⋯, θ_n−1,n, Δϖ_n−1,n) may not be identical to the one by (θ₁₂ + 2π, Δϖ₁₂, ⋯, θ_n−1,n, Δϖ_n−1,n). For instance, the resonant angles for the 2:3:4 resonant chain are $$\{\begin{array}{ll}{\theta }_{12}=3{\lambda }_2-2{\lambda }_1-{\varpi }_1,& \mathrm{\Delta }{\varpi }_{12}={\varpi }_1-{\varpi }_2,\\ {\theta }_{23}=4{\lambda }_3-3{\lambda }_2-{\varpi }_2,& \mathrm{\Delta }{\varpi }_{23}={\varpi }_2-{\varpi }_3.\end{array}$$(A.3)

And the resonant angles between planets m₁, m₃ $${\theta }_{13}=2{\lambda }_3-{\lambda }_1-{\varpi }_1,\mathrm{\Delta }{\varpi }_{13}={\varpi }_1-{\varpi }_3$$(A.4)

can be linearly expressed by θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃ as $${\theta }_{13}=\frac{1}{2}({\theta }_{12}+{\theta }_{23}-\mathrm{\Delta }{\varpi }_{12}),\mathrm{\Delta }{\varpi }_{13}=\mathrm{\Delta }{\sigma }_{12}+\mathrm{\Delta }{\varpi }_{23}.$$(A.5)

In this case, θ₁₂, θ₂₃, and Δϖ₁₂ are not 2π-periodic angles, but 4π-periodic.

Generally, we define the degeneracy of an angle ϕ (ϕ stands for any one of θ_ij, Δϖ_ij) as the smallest positive integer q which satisfies the following condition: $$\begin{array}{rl}& {\overline{\mathcal{H}}}_1({\theta }_{12},\cdots ,\varphi ,\cdots ,{\theta }_{n-1,n},\mathrm{\Delta }{\varpi }_{n-1,n})\\ & ={\overline{\mathcal{H}}}_1({\theta }_{12},\cdots ,\varphi +2q\pi ,\cdots ,{\theta }_{n-1,n},\mathrm{\Delta }{\varpi }_{n-1,n}).\end{array}$$(A.6)

Therefore, for the 2:3:4 resonant chain, the degeneracy of the angles (θ₁₂, Δϖ₁₂, θ₂₃, Δϖ₂₃) are (2, 2, 2, 1). The topology of the averaged Hamiltonian with e₁ = e₂ = e₃ = 10⁻⁴ on angular plane is shown in Fig. A.1, in which the periodicity confirms the degeneracy of resonant angles. From the plots in Fig. A.1, we also note that no stable symmetric configuration exists for near circular orbit.

According to the form of averaged Hamiltonian, the asymmetric ACRs always appear in pair. If an ACR locates at (θ₁₂, Δϖ₁₂, ⋯, θ_n−1,n, Δϖ_n−1,n), one must find an ACR at (−θ₁₂, −Δϖ₁₂, ⋯, −θ_n−1,n, −Δϖ_n−1,n).

Fig. A.1 Topology of the averaged Hamiltonian for the 2:3:4 resonant chain on (a) the (θ12, Δϖ12) plane and (b) the (θ23, Δϖ23) plane. The solid black and red circles indicate the location of the stable and unstable ACRs, respectively.

Appendix B Extension to four-planet resonant chain: Kepler-223 as an example

The methods introduced in this paper can be extended to resonant chains of more planets. As an example, we briefly present our analyses on the four-planet system Kepler-223, in which the orbital periods are close to T₁: T₂: T₃: T₄ = 3:4:6:8. The following resonant angles are adopted to characterise the resonant configuration. $$\{\begin{array}{ll}{\theta }_{12}=4{\lambda }_2-3{\lambda }_1-{\varpi }_1,& \mathrm{\Delta }{\varpi }_{12}={\varpi }_1-{\varpi }_2,\\ {\theta }_{23}=3{\lambda }_3-2{\lambda }_2-{\varpi }_2,& \mathrm{\Delta }{\varpi }_{23}={\varpi }_2-{\varpi }_3,\\ {\theta }_{34}=4{\lambda }_4-3{\lambda }_3-{\varpi }_3,& \mathrm{\Delta }{\varpi }_{34}={\varpi }_3-{\varpi }_4.\end{array}$$(B.1)

After some algebra, we obtain the degeneracy of these resonant angles as (3, 3, 6, 2, 2, 1), which means that the searching for resonant configuration will be done in $\mathcal{T}$⁶ = [0, 6π] × [0, 6π] × [0, 12π] × [0, 4π] × [0, 4π] × [0, 2π]. The numerical averaging over the synodic angle Q reduces the final Hamiltonian to 6 DoF with two free parameters, the spacing parameter L and the total angular momentum G, $$L=\sum _{i=1}^4\frac{{\mathrm{\Lambda }}_i}{k_i},G=\sum _{i=1}^4{\mathrm{\Lambda }}_i-{\mathrm{\Gamma }}_i,$$(B.2)

where Λ_i and Γ_i are the Poincaré variables as in Eq. (8).

We only consider here the stable ACRs that correspond to libration in the resonance. Taking the mass parameters of Kepler-223 from Mills et al. (2016), and following the geometric method described in the main text, we found 6 (= 2 × 3) families of stable ACR. In Fig. B.1 we summarise the locations and configurations of these ACR families as the function of the scaled AMD, defined in the same way as Eq. (18), $\delta ={\sum }_{i=1}^4{\mathrm{\Gamma }}_i/L$. We also plot the Laplace angles φ₁₂₃ = λ₁ + λ₃ − 2λ₂, φ₂₃₄ = λ₂ + 2λ₄ − 3λ₃ for adjacent three planets. It’s worth noting that the six asymmetric solutions we obtained actually take only three configurations, and each configuration contains two branches that are symmetric with respect to each other. That’s why only three lines in each panel for eccentricities in Fig. B.1 can be seen.

Our results obtained from the averaged Hamiltonian are consistent with the analytical results (Delisle 2017). Since the analytical model in Delisle (2017) is truncated at the first order of eccentricity, it gives constant critical angles at the ACR solutions. However, our method provides more precise results. Some resonant angles θ_ij, Δϖ_ij are found to experience large variations along the ACR solutions, while the three-planet Laplace angle φ₁₂₃ remains almost constant.

Fig. B.1 Locations and resonant configurations of six families of stable ACR for the four-planet resonant chain in Kepler-223. The variations of the eccentricities, the resonant angles, and the three-planet Laplace angles, are plotted against the scaled AMD (δ). Each family together with its symmetric opposite (see text) is represented by a specific colour in all panels.

References

All Tables

All Figures

Fig. 1 Topology of the Hamiltonian on different representative planes, (a) The topology on period ratio plane (n1/n2, n2/n3) with fixed eccentricities. (b) The topology on eccentricity plane (e1,e2). The stationary solution located at (n1/n2, n2/n3) = (2, 1.5) and (e1, e2, e3) = (0.00458, 0.00938, 0.00564) is marked by a red dot on both panels. The corresponding critical angles are (θ12, Δϖ12, θ23, Δϖ23) = (π, 0, π, π). The colour indicates the value of Hamiltonian, light yellow for high and dark blue for low. The same colour code is adopted in all figures in this paper, unless otherwise stated.

In the text

	Fig. 2 Different topologies on the (e1, e2) plane for different symmetric configurations (a) (θ12, Δϖ12, θ23, Δϖ23) = (0, π, 0, π) and (b) (π, π, 0, 0). The local extrema are (e1, e2, e3) = (0.119,0.1352,0.119) (red solid circle) and (0.300, 0.0662, 0.0671) (red cross), respectively.
In the text

Fig. 3 Topology on representative angular planes. The eccentricities of the symmetric ACR are (e1, e2, e3) = (0.00458, 0.00938, 0.00564). (a) Hamiltonian level curves on (θ12, Δϖ12) plane with (θ23, Δϖ23) fixed to (π, π), the saddle-like point at (θ12, Δϖ12) = (π, 0) is marked by a red cross. (b) Topology on (θ23, Δϖ23) plane with (θ12, Δϖ12) fixed to (π, 0). The node-like point at (θ23, Δϖ23) = (π, π) is marked by a red dot. In addition, the stable symmetric configuration (θ12, Δϖ12, θ23, Δϖ23) = (0, π, 0, π) is shown as blue dots.

In the text

Fig. 4 Symmetric ACR families S2, S3, S5, and S6 in the 1:2:3 resonant chain for small-eccentricity configurations on (a) the (e1 cos θ1, e2 cos θ2) plane and (b) the (e2 cos θ2, e3 cos θ3) plane. The stable symmetric ACRs are in blue while the unstable ones in red. The notations of the families are the same as in Antoniadou & Voyatzis (2022), and the resonant angles of different configurations are presented in Table 1.

In the text

Fig. 5 Topology on representative angular planes for symmetric ACR S6 located at (e1e2,e3) = (0.0569, 0.0846, 0.0531). (a) Topology on (θ12, Δϖ12) plane with fixed (θ23, Δϖ23) = (0, π). The saddle-like point at (θ12, Δϖ12) = (0, π) marked by a cross indicates the instability of this solution, while two asymmetric ACRs bifurcated from this unstable symmetric one are indicated by black dots. (b) The topology on (θ23, Δϖ23) plane with (θ12, Δϖ12) fixed to (0, π). The node-like point at (θ23, Δϖ23) = (0, π) is marked by a pink dot.

In the text

Fig. 6 Locations and resonant configurations of the symmetric ACR family (S6) and the asymmetric ACR family (A6). In the upper two panels (a) and (b), the stable symmetric, unstable symmetric, and stable asymmetric ACRs are indicated by blue solid, red solid, and black dashed lines, respectively. Three solid circles on these lines indicate the initial conditions that will be analysed in detail later (see text). The lower two panels (c) and (d) show the evolutions of resonant angles along the asymmetric ACR as the function of δ (scaled AMD, see text). The resonant angles θ12, θ23 and the secular angles Δϖ12, Δϖ23 are labelled by the corresponding lines.

In the text

Fig. 7 Law of structure for the 1:2:3 resonant chain. The left (right) column of panels show the relation between the mean motion ratio n1/n2 (n2/n3) and the eccentricities (e1, e2, e3 from top to bottom). Blue lines are for the stable symmetric ACR families (S6) while the black lines for the stable asymmetric family (A6). On each panel, the lowest line is for the planets in the Kepler-51 system, while the middle and the top lines are for the cases of 10 and 100 times more massive planets, respectively (see text). The critical points, at which the S6 family becomes unstable and the asymmetric solution appears, are denoted by red crosses.

In the text

	Fig. 8 Variation of the Hamiltonian in two synodic periods $2Q_{12}^T=12\pi$ for (a) the symmetric configuration and (b) the asymmetric configuration. The vertical line denotes the position when the Hamiltonian attains the minimum.
In the text

	Fig. 9 χ2 map for the standard pendulum model, (a) Phase structure (energy level) on (φ, pφ) plane. The broad red curve marks the separatrix between the libration and circulation region, (b) χ2 value in logarithm indicated by colour.
In the text

	Fig. 10 Topology (contour line) and dynamical map (colour map) on the (e1, e2) plane for Case A, with constants L and G determined by the initial condition on the stable ACR solution (e1, e2, e3, θ12, Δϖ12, θ23, Δϖ23) = (0.0124, 0.0237, 0.0128, 0, π, 0, π), which is indicated by a red dot on all maps. The eccentricity variations Δei are used as the stability indicator as labelled in each panel.
In the text

Fig. 11 Resonant structure indicated by χ2 on the (e1, e2) plane for Case A. Four panels are for the two-body MMRs between (a) m1 and m2, (b) m2 and m3, (c) m1 and m3, and (d) the three-planet Laplace resonance with resonant angle φ0. The large value of χ2 (light colour) indicates that the resonant angle is in libration, while the small χ2 (dark tones) indicates circulation. The narrow strip-like structure indicated by a black arrow in (a) divides the map into two parts, namely Region 1 and Region 2. Moreover, another narrow strip (indicated by the red arrow) can be seen in Region 2.

In the text

	Fig. 12 Typical behaviour of the orbits in Region 1 (Fig. 11a). The temporal variations of n1/n2, n2/n3, e1, e2, e3 are in left panels, while right panels are for resonant angles θ12, Δϖ12, θ23, Δϖ23, φ0. The orbits are librating quasi-periodically around the stable symmetric ACR solution S6. The eccentricities of S6 are indicated by horizontal dashed lines.
In the text

	Fig. 13 Same as Fig. 12 but for Region 2 in Fig. 11a.
In the text

Fig. 14 Same as Fig. 10 but for Case B. The constants L and G are determined by the initial condition on the stable asymmetric ACR solution (black solid circle in Fig. 6). (a) Around the unstable symmetric ACR solution with the same L, G. The corresponding eccentricities e1, e2, e3 are 0.0569, 0.0846, 0.0531 (red solid circle in Fig. 6). (b) Around the stable asymmetric ACR solution. The colour in both panels indicates the value of max(Δe).

In the text

	Fig. 15 Horseshoe-like orbit near the unstable symmetric ACR. The unstable symmetric ACR and the stable asymmetric ACRs are indicated by the red and black horizontal lines, respectively.
In the text

	Fig. 16 Same as Fig. 11 but for Case B. The map is divided into two parts (Region 1 and Region 2) by a stripe-like structure of relatively smaller χ2 value, as indicated by the arrow in panel a.
In the text

	Fig. 17 Same as Fig. 13 but the initial conditions are from the high eccentricity region of Fig. 14 (see text).
In the text

Fig. 18 Formation and evolution of the 1:2:3 resonant chain via convergent planetary migration. The left panels show the temporal evolutions of period ratios and eccentricities, and the critical angles are shown in the right panels. The planetary system is captured in symmetric configuration S6 firstly, subsequently evolves along the A6 family and finally librates around the asymmetric configuration.

In the text

	Fig. 19 Same as Fig. 18 but for a ‘joint two-body MMR’ case (see text). In the final resonance configuration, Δϖ12 and φ0 circulate while θ12, θ23, and Δϖ23 librate. The vertical lines in the left panels indicate two moments t1 = 1.2 Myr and t2 = 1.8 Myr.
In the text

	Fig. 20 Resonant capture and evolution of the 2:3:5 resonant chain on the eccentricities plane. The stable (θ12, Δϖ12, θ23, Δϖ23) = (0, π, π, π) ACR is indicated by blue dashed lines while the grey dots show the migrating path of the system evolving to higher eccentricities.
In the text

	Fig. 21 Same as Fig. 20 but for the 3:5:7 resonant chain. The blue dashed lines indicate the location of the stable ACR of (θ12, Δϖ12, θ23, Δϖ23) = (π, π, π, π).
In the text

	Fig. 22 Resonant capture and evolution of the 2:3:4 resonant chain. The black lines indicate the location of the asymmetric ACR solution arising from the circular configuration. The grey dots stand for the capture and evolution of planets by the resonant chain in the numerically simulated convergent migration.
In the text

	Fig. A.1 Topology of the averaged Hamiltonian for the 2:3:4 resonant chain on (a) the (θ12, Δϖ12) plane and (b) the (θ23, Δϖ23) plane. The solid black and red circles indicate the location of the stable and unstable ACRs, respectively.
In the text

	Fig. B.1 Locations and resonant configurations of six families of stable ACR for the four-planet resonant chain in Kepler-223. The variations of the eccentricities, the resonant angles, and the three-planet Laplace angles, are plotted against the scaled AMD (δ). Each family together with its symmetric opposite (see text) is represented by a specific colour in all panels.
In the text