© The Authors 2024

1 Introduction

The long-term dynamical evolution of planetary systems is defined by what is called their secular dynamics. Some planetary systems are occasionally known to be in mean-motion resonance (MMR), which means that their orbital periods are related in a simple way. The secular dynamics of planetary systems within MMR are generally very different from those outside MMR (Beaugé & Michtchenko 2003; Callegari et al. 2004; Batygin & Morbidelli 2013).

Methods for studying the secular dynamics of planetary systems within an MMR must initially contemplate and correctly reproduce the dynamics of the MMR. Some methods for planetary systems in MMR have been developed using analytical methods (Batygin & Morbidelli 2013; Quillen & French 2014), some were pure numerical studies (Haghighipour et al. 2003), and semi-analytical methods were also proposed (Michtchenko et al. 2008; Gallardo et al. 2021).

The secular evolution of planetary systems within MMRs is a more challenging problem (Batygin & Morbidelli 2013). The restricted problem is a simpler case, that is, a particle that is perturbed by an unperturbed planet. In our first paper (Pons & Gallardo 2022), we proposed a method for following the secular dynamics of a particle in deep resonance (i.e., zero libration amplitude) with a planet. In this paper, we generalize the method to the study of a two-planet system in deep MMR, allowing us to obtain the long-term evolution of the planetary system without restrictions on the eccentricity.

This work is organized as follows. Section 2 describes the theoretical framework, starting with the Hamiltonian formalism, the calculation formulae for the equilibrium points, and a brief description of the invariant adiabatic principle that allows the application of the model on secular timescales. Then, Sect. 3 describes the method we used to apply the model in the most generic and comprehensive way. There, we recall the definition and construction of the phase-space graphical representation we call ℌ surface. Finally, we describe some nomenclature that we used to distinguish and classify the apsidal corrota-tion resonances (ACRs). Section 4 presents the most interesting results. This section is divided into two subsections, the first of which presents the generic results found for 1:1, 2:1, 3:1, and 3:2 MMRs, and the second subsection contains the results we achieved by applying the model to two real exoplanetary systems. In Sect. 5, we present our conclusions.

2 Theoretical framework

2.1 Semi-analytical theory

As is commonly used, index 1 is reserved for the inner and index 2 for the outer planet. This means that a₁ ≤ a₂. Index 0 is reserved for the star. The method we describe in this work is valid for coplanar orbits of planets in deep MMR with arbitrary masses. In order to obtain some simplified expressions below, we assume that m_1,2 << m₀.

Following the model developed by Gallardo et al. (2021), we can express the canonical elements of the planets in a coplanar configuration as follows: $$\begin{array}{l}{\lambda }_i;L_i={\beta }_i\sqrt{{\mu }_i{a}_i}\\ {\varpi }_i;{\text{Γ}}_i=L_i\left(\sqrt{1-e_i^2}-1\right),\end{array}$$(1)

where β_i = m₀m_i/(m₀ + m_i) is the reduced mass, µ_i = k²(m₀ + m_i), a_i is the ith planet semi-major axis, e_i is the eccentricity, and k is the gravitational Gauss constant. With these variables, the Hamiltonian takes the next form, $$H=-\frac{{\mu }_1^2{\beta }_1^3}{2L_1^2}-\frac{{\mu }_2^2{\beta }_2^3}{2L_2^2}-R,$$(2)

where R is the disturbing function, which can be expressed in the following way: $$R=\frac{k^2{m}_1{m}_2}{\left|r_2-r_1\right|}-\frac{m_1{m}_2}{m_0}{v}_1\cdot v_2.$$(3)

Here, r_i, v_i are the planetary position and velocities in the Poincare reference system. The nonperturbative Hamiltonian is given by the term ℌ₀ = H + R.

Our next step is to introduce the so-called critical angle σ = k₁λ₁ − k₂ λ₂ + (k₂ − k₁ )ϖ₁. In order to do so, we apply a canonical transformation (see Appendix A) and obtain the following set of canonical variables: $$\begin{array}{l}\sigma ;I_1=L_1/k_1\\ {\lambda }_2;I_2=L_2+k_2{L}_1/k_1\\ \text{Δ}\varpi ;{\mathcal{K}}_1={\text{Γ}}_1-\left(k_2-k_1\right)L_1/k_1\\ {\varpi }_2;{\mathcal{K}}_2={\text{Γ}}_1+{\text{Γ}}_2-\left(k_2-k_1\right)L_1/k_1.\end{array}$$(4)

This transformation also introduces the variable ∆ϖ = ϖ₁ − ϖ₂, which is useful for studying the ACRs in planetary systems (Beaugé et al. 2003; Michtchenko et al. 2008).

To conclude the basis of our theoretical development, we calculate the averaged disturbing function 𝓡 (see Appendix A) following the approach given in Gallardo (2020), that is, we make an average in λ₂ assuming that on the resonant timescale (which is longer than the timescale on which the average is performed) a₁, a₂, e₁, e₂, and ∆ϖ remain constant. This is a good approximation as long as the eccentricities are not too close to 0, because the closer they are to 0, the larger the longitudes of the perihelion change rate.

After the averaging, the Hamiltonian no longer depends on λ₂. Therefore, I₂ becomes a constant of motion, correlating the evolution of the semi-major axis. Additionally, when we apply D’Alembert rules to a generic argument of the form ϕ = k₁λ₁ + k₂λ₂ + k₃ϖ₁ + k₄ϖ₂, it can be demonstrated that 𝓡 depends on σ and ∆ϖ, but not on ϖ₂, implying that 𝓚₂ is another constant of motion. Based on these arguments, the Hamiltonian takes the following form: $$\mathscr{H}={\mathscr{H}}_0\left(I_1;I_2\right)-ℛ\left(\sigma ,I_1,{\mathcal{K}}_1,\text{Δ}\varpi ;I_2,{\mathcal{K}}_2\right),$$(5)

where all the values before the semicolon are variables, and all values after the semicolon are considered fixed parameters. These dependences of the resonant Hamiltonian variables are valid on secular timescales.

2.2 Equilibrium points

Under the assumption of deep MMR condition, the two semi-major axis are constant, that is, L₁ and L₂ are also constant (see Sect. 2.3). In addition, on resonant timescales, as we previously stated, e₁, e₂, and ∆ϖ can also be considered constant, that is, 𝓚₁ and 𝓚₂ are constant. This allows us to write that $$\mathscr{H}={\mathscr{H}}_0\left(I_1;I_2\right)-ℛ\left(\sigma ;I_1=I_{1\text{ nom }},I_2,{\mathcal{K}}_1,{\mathcal{K}}_2,\text{Δ}\pi \right),$$(6)

where I₁ = I_1nom (where nom stands for nominal, which refers to the values for the exact resonance, given by Eq. (8)), and the remaining values after the semicolon are considered fixed parameters, at least on the resonant timescale. We consider I₁ to be a variable in the nonperturbative term, but a constant in the perturbative term. This manipulation does not change the results and simplifies the search for the equilibrium point.

To obtain the equilibrium points under these further assumptions, we used the canonical equations. Therefore, we have that the equilibrium points must satisfy the following conditions: $$\frac{\partial {\mathscr{H}}_0}{\partial I_1}=0;\frac{\partial ℛ}{\partial \sigma }=0.$$(7)

The first equation leads to n₁k₁ = n₂k₂, where n_i describes the planetary mean motions. This is the deep resonant condition, which is equivalent to the following formula: $$a_1=a_2{\left(\frac{m_0+m_1}{m_0+m_2}\right)}^{1/3}{\left(\frac{k_1}{k_2}\right)}^{2/3}\equiv a_{1\text{nom}},$$(8)

where for a given a₂, a₁ is defined by the masses and by the MMR k₂:k₁ . We refer to this value as the nominal semi-major axis of the resonance.

The second equation indicates that all the equilibrium points are located in the extrema of the 𝓡(σ) function. In particular, it can be demonstrated (see Appendix A) that the stable equilibrium points occur in the minima of this function. We assume that a minimum occurs in σ_n. Therefore, we can say that there is an stable equilibrium point in (I₁, σ) = (I_1nom, σ_n). In resonant dynamics, the stable equilibrium points are also known as libra-tion centers. In general, for a given set of (a₁, a₂, e₁, e₂, ∆ϖ), N ≥ 1 libration centers (I_1nom, σ_n) can exist with n = 1…N. For better graphical representations, we occasionally use an alternative for the critical angle, which is defined as follows: $$\theta =k_1{\lambda }_1-k_2{\lambda }_2.$$(9)

A priori, all libration centers are valid. Nevertheless, when a close encounter between the planets occurs, the evolution can become unstable, leading to strong changes in the orbital elements, including a_i. This implies that the MMR is broken. In order to dismiss these stable equilibrium points in the encounter condition, we use Hill’s criteria. Consequently, we disregard a libration center when the following condition is met: $$\left|r_2-r_1\right|<\xi R_{\text{Hill }},$$(10)

where R_Hill is the planetary Hill radius, which is defined as follows (Gladman 1993): $$R_{\text{Hill }}=\frac{a_1+a_2}{2}{\left(\frac{m_1+m_2}{3m_0}\right)}^{1/3}.$$(11)

For practical use, we chose in general a ξ value from 3 to 4, depending on the particular studied case. A more complete analysis of ξ values can be found in Gallardo et al. (2021).

Fig. 1 Averaged resonant Hamiltonian contour curves in plane (a2, σ) in 3:1 MMR for a system with m1 = m2, e1 = 0,47, e2 = 0,38, and ∆ϖ = 84°.

2.3 Secular model

With the purpose of deriving a model that is applicable to secular timescales, we used the adiabatic invariant principle. In the celestial mechanics context, in particular, in planetary dynamics, this principle states that as long as the variations in e_i and ∆ϖ are slow enough, a quantity J exists that is related to faster variables (I₁ and σ). This variable remains constant in time. This quantity is known as the adiabatic invariant and is defined as follows: $$J={\displaystyle \oint I_1}d\sigma$$(12)

This definition can be interpreted as the enclosed area inside one of the curves of the Hamiltonian level shown in Fig. 1. The mathematical details of which exactly mean that the variations in e_i and ∆ϖ are slow enough can be found in Henrard (1993). In simple terms, we assumed that e_i and ∆ϖ change much more slowly (i.e., adiabatically) than I₁ and σ.

For practical reasons, we assumed the ideal case of J = 0, as in Pons & Gallardo (2022). This implies that the area of the enclosed curve in the (I₁ ,σ) plane is 0. In other words, we are located exactly in the libration center. In this situation, the planetary resonant angles librate around the equilibrium value, but with null resonant amplitude, which describes the deep MMR hypothesis we mentioned in Sect. 2.1.

The value of a given libration center σ_n can change as a consequence of e_i and ∆ϖ variations (we recall Eqs. (6) and (7)) that are present in the secular evolution. When no encounter occurs, the adiabatic invariant principle ensures that the system continues in a deep MMR condition.

3 Method

3.1 AMD conservation

Since we have that 𝓚₂ is a constant of motion and L₁ is assumed constant (because of the deep MMR hypothesis), we conclude that −(Γ₁ + Γ₂) is also constant. This is the quantity known as angular momentum deficit (AMD), which is defined as follows (Laskar 1997): $$\text{AMD}={\displaystyle \sum _{i=1}^2{\beta }_i}\sqrt{{\mu }_i}\left(1-\sqrt{1-e_i^2}\right).$$(13)

Petit et al. (2017) concluded that in the presence of MMR, there is no guarantee that AMD is conserved. This is mainly due to the chaos that emerges when an MMR overlap occurs, particularly in high-eccentricity domains. We assumed that there is no MMR overlapping, which is reasonable in the deep-resonance hypothesis. In other words, we assumed that AMD is conserved even in a MMR.

3.2 Eccentricity domain

If the AMD is conserved, then the following quantity is also conserved: $$\mathcal{A}ℳ={\beta }_1\sqrt{{\mu }_1{a}_1\left(1-e_1^2\right)}+{\beta }_2\sqrt{{\mu }_2{a}_2\left(1-e_2^2\right)}.$$(14)

This quantity becomes the angular momentum of the system when m_1,2 << m₀ (Michtchenko et al. 2008). Regardless of the case, 𝒜𝓜 always has the same functional form, which is $$\mathcal{A}ℳ=C_1\sqrt{1-e_1^2}+C_2\sqrt{1-e_2^2},$$(15)

where $C_1={\beta }_1\sqrt{{\mu }_1{a}_1}\text{ and }{C}_2={\beta }_2\sqrt{{\mu }_2{a}_2}$ are constants because they only depend on the masses and semi-major axis (which are constant because of the deep MMR hypothesis). It is convenient to define a normalized 𝒜𝓜 as follows: $$\mathcal{A}{ℳ}_{\text{norm}}=\frac{\mathcal{A}ℳ}{\mathcal{A}{ℳ}_{\max }}=\frac{\sqrt{1-e_1^2}}{1+\eta }+\frac{\sqrt{1-e_2^2}}{1+{\eta }^{-1}},$$(16)

where 𝒜𝓜_max = C₁ + C₂, and η = C₂/C₁ is a new parameter we defined. With some simple operations, we obtained $$\eta =\frac{m_2\sqrt{\left(m_0+m_2\right)a_2}}{m_1\sqrt{\left(m_0+m_1\right)a_1}}\simeq \frac{m_2}{m_1}\sqrt[3]{\frac{k_2}{k_1}},$$(17)

where the last equality is valid only when m_1,2 << m₀. In this condition, the functional form of 𝒜𝓜_norm only depends on the planetary mass ratio and on the MMR in which the masses are locked. Figure 2 shows an example of the 𝒜𝓜_norm function and its contour curves. These curves allow a rapid identification of the eccentricity domains, that is, the possible variation ranges of e₁ and e₂.

In order to inclusively and generically study the secular dynamics of resonant coplanar planets with eccentric orbits, we assumed 𝒜𝓜_norm and η as free parameters. A priori, a double sweep in these two parameters might be performed. Nevertheless, in the context of planetary dynamics, we considered a predefined set of MMRs and a set of mass ratios that determined the values for η. On the one hand, the MMR set comprised the most typical and stronger resonances, which are the 2:1, 3:2, 3:1, and 1:1 MMR. For the ratio m₂/m₁, we considered the values 1/20, 1/5, 1, 5, and 20. From the combination of these two sets, we obtained 20 different values for η that ranged from 0.5 to almost 30. On the other hand, the 𝒜𝓜_norm values were chosen (after η was defined) after an inspection of the contour curves, with the aim of selecting the cases for which a higher eccentricity excursion appeared possible. This condition was usually met for 𝒜𝓜_norm values between 0.7 to 0.95. We present the most interesting cases hereafter studying all those that resulted from our selection criteria.

After 𝒜𝓜_norm was fixed, the value of one eccentricity determined that of the other. For instance, for a value of e₁, e₂ is given by the formula below, $$e_2=\sqrt{1-{\left[\left(1+{\eta }^{-1}\right)\mathcal{A}{ℳ}_{\text{norm }}-{\eta }^{-1}\sqrt{1-e_1^2}\right]}^2}.$$(18)

This equation allows drawing two restrictions. The first restriction emerges just from the fact that the term inside the radicals must be positive. Therefore, $$\mathcal{A}{ℳ}_{\text{norm }}\le \frac{\sqrt{1-e_1^2}+\eta }{1+\eta }.$$(19)

The second restriction arises when we inspect the term inside the straight parentheses in expression (18). This term must be also positive because if it is not, the second term in Eq. (16) would be negative. For this reason, we have $$\mathcal{A}{ℳ}_{\text{norm }}\ge \frac{\sqrt{1-e_1^2}}{1+\eta }.$$(20)

These two limiting conditions for the eccentricity variation ranges are useful in the next section, when we detail the method we used to explore the phase space.

Fig. 2 Normalized angular momentum function 𝒜𝓜norm(e1, e2) for 2:1 MMR with m2 = m1. (a) Three-dimensional plot of 𝒜𝓜norm(e1, e2). (b) 𝒜𝓜norm(e1, e2) contour curves with numerical values above the curves.

3.3 ℌ surfaces

Exploiting the fact that we work with conservative systems, we constructed some phase portraits that are very useful because they contain the constant Hamiltonian contour curves. When a system has more than one degree of freedom, these phase portraits are constructed considering one coordinate and its conjugate momenta as variables and the rest as parameters. In our case, the averaged Hamiltonian has three degrees of freedom a priori, which translates into six variables that essentially are a₁, a₂, e₁, e₂, ϖ₁, and ϖ₂. On the one hand, as a result of the deep MMR condition, a₁ and a₂ were assumed constant, and as we mentioned previously, 𝓡 depends on Am and σ (not on ϖ₁ and ϖ₂ separately). On the other hand, AM conservation allows us to establish a relation between e₁ and e₂. Therefore, we have essentially three variables of interest to understand the secular evolution in the resonant context, which are e_i (where the index i takes the value 1 or 2), ∆ϖ, and σ. For the three variables, we used the graphical representation developed by Pons & Gallardo (2022). This representation allows us to build one single phase portrait that contains all the information of the secular evolution of the system instead of several bidimensional phase portraits with one of the three variables as a parameter. The single phase portrait containing all the information is called 𝓗 surface.

Pons & Gallardo (2022) introduced the 𝓗 surface for the restricted case of the coplanar and resonant three-body problem. In the present research, we need two 𝓗 surfaces, one for each planet. Consequently, one surface is associated with the variables (e₁, ∆ϖ, σ), and the other surface with (e₂, ∆ϖ, σ). We call these phase portraits surface 𝓗₁ and surface 𝓗₂. They are related through the 𝒜𝓜 conservation because this condition relates e₁ with e₂ (we recall expression (15)). As a consequence, the topology of the two surfaces is the same, but they can vary somewhat in size or scale. For this reason, it is usually sufficient to study one of them to understand the system evolution qualitatively.

3.4 Apsidal corotation resonances

In terms of secular evolution, there are two possible behaviors for the angular variable Am. The first behavior is circulation, where ∆ϖ eventually takes on all the possible values for an angular variable and essentially maintains the same sign as the time derivative. Alternatively, ∆ϖ can librate (or oscillate) around a certain value. In this last case, the two planets are in ACR (Beaugé et al. 2003; Zhou et al. 2004). If this libration has null amplitude, the angle between each the line of the apsides of each orbit is fixed and the two eccentricities also remain unchanged. The situation in which each orbit is frozen with respect to the other persists until some external mechanism removes the planets from the exact ACR point.

The ACRs are classified into two groups: symmetric ACRs, and the asymmetric ACRs. The symmetric ACRs occur when ∆ϖ is 0° or 180°. In asymmetric ACRs, ∆ϖ takes the remaining possible angular values. This ACR distinction is very frequently used in the literature. Is reasonable that asymmetric ACR have double multiplicity because of the geometrical symmetry of the problem, as we explain next. When the orientation of one orbit with respect to the other implies that the relative position between them is frozen on secular timescales, then the evolution should be the same when we consider a new configuration of the orientation, making an axial symmetry with respect to the symmetry of the apsis lines. Moreover, the planetary position in the orbit should also be symmetrical to the line of aspsis. The mathematical counterpart of these statements is that when M is the mean anomaly, when an ACR exists in (∆ϖ, M₁) = (∆ϖ_α, M_1α) (we assumed M₂ = 0°), then another ACR exists in (∆ϖ, M₁) = (−∆ϖ_α, −M_1α). With some simple operations, the same holds for the critical angle, that is, when the first ACR has σ = σ_α, then the other ACR has σ = −σ_α. A numerical example would be that when an ACR exists in ∆ϖ = 45° and σ = 60°, then another asymmetric ACR must exist in ∆ϖ= 360 − 45 = 315° and σ = 360 − 60 = 300° (here we summed 360° just to work with positive values). In order to simplify the ACR count, we always refer to a group of two ACRs as one asymmetric ACR.

We introduce another criterion below to classify ACRs in the framework of resonant dynamics. ACR type I or ACR_I are ACRs with a fixed libration center on the secular timescale. Conversely, ACR type II or ACR_II are ACRs with a changing libration center on secular timescale (in all studied cases, this change was always of a circulation type, but there is no reason a priori that prevents it from having a librational-type change in the nominal value of the ACR). In this case, the exact ACR may coincide with the edge of an 𝓗 surface. This implies that a system in an exact ACR_II cannot possibly maintain invariant J because it is too close to the edge of the 𝓗 surface (Pons & Gallardo 2022).

3.5 Numerical integrations

To perform numerical integrations, we made use of the Mer-curius algorithm included in the Rebound Python package, which is a hybrid symplectic integrator that is very similar to Mercury (Chambers 1999). It is basically composed of a WHFast algorithm that switches over smoothly to IAS15 whenever a close encounter occurs.

4 Results and discussion

4.1 General results

Without losing generality, we fixed some variables to explore the phase space. For the rest of the work, we have a₂ = 1 au and a larger planet with mass m_i = 1 M_J.

With the aim of testing and validating the model, we compared every phase portrait presented in this paper with numerical integrations of the exact equations of motion. A small correction was sometimes required to the initial value for a₁ given by Eq. (8) in order to ensure the deep MMR condition. This is mainly due to the short-term perturbations and to the law of structure (Ferraz-Mello 1988), which can slightly modify the resonant nominal semi-major axis. The MMR we studied first is presented below. In some of the integrations, we had to implement these corrections to ensure the zero-amplitude condition in the resonant librations.

4.1.1 MMR 2:1

In panel a of Fig. 3, we show two examples of 𝓗 surfaces, both regarding a planetary system locked in 2:1 MMR with 𝒜𝓜_norm = 0.9. The first system is a planetary system with m₂/m₁ = 0.2. In this case, there are two topologically separated surfaces in the phase portrait. When this occurs, we name them subsurfaces. One of the subsurfaces has a constant libration center (θ = 0°) and an ACR_I in (e₂, ∆ϖ, θ) = (0.88, 180°, 0°). The contour curves of the constant Hamiltonian are plotted in black, and the five numerical integrations run to test the model are shown in red (some initial a₁ values were modified up to 7% from the nominal value to obtain the deep MMR condition). Two integrations were required in the case of apsidal corotation resonance, and three were required with a circulating ∆ϖ.

The other subsurface contains an asymmetric ACRI located near ∆ϖ= θ = 0° (see the red marker in panel b of Fig. 3) and for a lower eccentricity value. In this particular case, the subsurface was filtered to remove areas in which a close encounter might occur. In one of these encounter zones lay another asymmetric ACR_I. However, as expected, numerical integrations with initial conditions in this zone resulted in a chaotic evolution. Similar to the other subsurface, we show in colors five 1 kyr numerical integrations, four of which are in ACR condition (but around the two close ACR points) and the remaining one is in a ∆ϖ circulating condition. This last integration shows a great e₂ excursion of at least 0.6 (green curve).

The second case presented here for the MMR 2:1 is a system with m₂/m₁ = 1. Figure 3b shows the resulting 𝓗₁ surface. Here, four ACR type I are distributed in two subsurfaces. Due to the circularity of the angular variables, it seems to be more than two. Notwithstanding, there are only two topologically separated subsurfaces. This time, nine numerical integrations were performed to compare with the model. Most of them were evolutions around an ACR, except for one (yellow curve), with a circulating ∆ϖ. The assymetric ACR with an integation in green is in a quasi-encounter situation because at this point |r₂ − r₁| ≃ 4R_Hill. This integration is stable, but for a displacement from the ACR initial conditions, the evolution remains stable for 400 yr at most and then shows an irregular behavior.

We recall that three topological types of constant 𝓗 curves are possible in these surfaces: closed ones with a librating ∆ϖ (ACR condition), closed ones with a circulating ∆ϖ, and open curves (Pons & Gallardo 2022). The last type of curves implies a breach of the invariance principle that might hamper the comparison between the model and the numerical experiments and might in some cases lead to chaotic behavior. All the initial conditions of the numerical integrations shown here correspond to the first two groups of curves. The yellow curve is at the limit of being in a closed curve. If the initial e₁ were increased slightly, the evolution changes drastically turning into a chaotic behavior because an open curve family lies on the surface.

We found that the model contour curves and the numerical integrations in all the cases presented so far agree well. This is not a strictly mathematical proof that the model is correct, but it contributes to validating the model.

Fig. 3 Phase portraits for MMR 2:1 and 𝒜𝓜norm = 0.9. (a) 𝓗2 surface for a system with m2/m1 = 0.2. The ten numerical integrations run for 10 kyr are plotted in colors (some areas of the surface were removed for clarity). (b) 𝓗1 surface for a system with m2/m1 = 1. The nine numerical integrations run for 1 kyr are shown.

Fig. 4 𝓗1 (e1, ∆ϖ) contour curves for the 3:2 MMR for two different systems. (a) First system: m2/m1 = 1 and 𝒜𝓜norm = 0.95. (b) Second system: m2/m1 = 5 and 𝒜𝓜norm = 0.9. In both panels, the numerical integrations are drawn in color. The red crosses indicate the initial condition in each case.

4.1.2 MMR 3:2

In this section, we present two examples in MMR 3:2. The first example consists of a planetary system with m₂/m₁ = 1 and 𝒜𝓜_norm = 0.95. We focused on the secular evolution of e₁ and ∆ϖ and consequently disregarded the evolution of the libration center. This was achieved by projecting the 𝓗₁ surface onto the (e₁, ∆ϖ) plane. The result of this operation is shown in Fig. 4. In particular, panel a of the figure shows the evolution of the system. Two symmetric ACRs are visible. One is located in ∆ϖ= 180° and the other in ∆ϖ= 0°. The last ACR is type II and is placed inside the encounter zone, and therefore, the contour curves in the proximity of the ACR are missing. In this case in particular, surface areas where |r₂ − r₁ | < 4.5R_Hill were removed. The evolution in an encounter zone in general is related to abrupt variations in the semi-major axis (which obviously disrupt the resonant condition) and strong changes in eccentricity, which can lead to a collision with the star or to an ejection from the system. In between the two ACRs lies a transition zone where ∆ϖ circulates.

The second system is characterized by m₂/m₁ = 5 and 𝒜𝓜_norm = 0.9. In this case, the dynamical structure is more complex, as shown in Fig. 4b. It shows four symmetric ACRs, three of which are type I, and the remaining ACR is type II. The transition areas between ACRs seem to be distorted because 𝓗₁ overlaps slightly with itself. This constitutes a practical problem when a projection is intended.

The numerical integration and the model contour curves in the two systems again agree well. In these examples, eccentricity excursions up to 0.4 at least can occur.

4.1.3 MMR 3:1

As in the previous MMRs, we present two examples for the 3:1 MMR as well. The phase-space structure is shown in Fig. 5. Panel a shows the phase portrait of a system with two equal-mass planets and 𝒜𝓜_norm = 0.8. The resulting 𝓗₁ surface is similar to the surface shown in Fig. 3a because in both cases, a symmetric ACRI and an asymmetric ACRI lie close to ∆ϖ= σ = 0°. They differ basically in the position of the eccentricity value.

The second example consists of a system with m₂/m₁ = 0.2 and 𝒜𝓜_norm = 0.99. This time, there is an asymmetric ACR_I and, for lower e₁ values, a symmetric ACR_II around which the libration center circulates on secular timescales.

Fig. 5 Phase portraits for MMR 3:1. (a) 𝓗1 surface for a system with m2/m1 = 1 and 𝒜𝓜norm = 0.8. The ten numerical integrations run for 10 kyr are plotted in colors. (b) 𝓗1 surface for a system with m2/m1 = 0.2 and 𝒜𝓜norm = 0.99. The five numerical integrations run for 1 kyr are shown.

4.1.4 MMR 1:1

From all the cases investigated in this MMR, we present the most interesting, which consists of a system with m₂ = m₁ and 𝒜𝓜_norm = 0.9. In this case, there are three subsurfaces with negligible libration center variation, that is, which are essentially planes located at σ ~ 0°, σ ~ 120° and σ ~ 240°. We therefore present just the Hamiltonian contour curves projected onto the (e_i, ∆ϖ) domain. On the one hand, in Fig. 6a, we show the contour curves of the case σ ~ 120^° , where two asymmetric ACRI lie. We only show this subsurface, but the other sub-urface around σ ~ 240°, has a similar aspect, but is flipped symmetrically with respect to the ∆ϖ= 0^° axis. This ACR distribution completely satisfies the explanation given in Sect. 3.4 for the pair multiplicity in all asymmetric ACRs. On the other hand, in Fig. 6b, we show a stable symmetric ACRI (the other one is in encounter condition). These five ACRs correspond to the classical Lagrangian solutions L₃, L₄, and L₅ and to the two points mentioned as anti-Lagrangian points in Giuppone et al. (2010, see Table 1 of that work). There, the authors carried out a very complete research in co-orbital planetary motion, where they explored the entire ACR families for different m₂/m₁ values.

4.2 Real exoplanetary systems

After investigating the secular dynamics of coplanar resonant three-body systems with a generic approach, we explored the exoplanet database¹ to find examples of resonant planet pairs with the aim of applying the model to them. Our intention was to focus on systems that might be in or near the deep MMR condition and with at least one planet with moderate to high eccentricity. In order to meet these criteria, we disregarded systems with a resonant offset (Charalambous et al. 2022) larger than 2%. Additionally, we filtered the database and selected systems with at least one planet with e > 0.2 and with a low uncertainty in their determination (∆e < 0. 05). After we searched the database, we found 22 exoplanet pairs. Most of them were analyzed, and we finally selected two of the more interesting systems. The first system is in a very common resonance, but in a deep MMR condition. The second system is in a very rare MMR and also very close to the deep MMR condition.

For a real exoplanet system locked in a k₂ : k₁ MMR, we calculated the 𝒜𝓜_norm from its planetary masses, semi-major axes, and eccentricities. The procedure of the 𝓗 surface construction was then analogous to the generic case. Finally, a numerical integration of the real system was carried out, and the result was compared directly with the 𝓗 surface contour curves.

For a given system, the initial conditions used for the numerical integration were obtained from the orbital elements for which observational data were fit. In the following sections, the corresponding paper is referenced for each analyzed system. At the moment of setting up the numerical integration, some adjustments to the initial conditions were needed for a more straightforward comparison with our model. On the one hand, our model focuses in the relative orientation of the orbits, that is, it considers ∆ϖ to describe the orientation. We therefore assumed ϖ₂ = 0° and calculated ϖ₁ = ϖ_1i − ϖ_2i, where variables with the subindex i refer to the data obtained from the original paper.

On the other hand, we also redefined the reference zero time by imposing M₂ = 0°. This implies that M₁ = (τ_2i − τ_1i) * (360°/P₁), where P₁ is the orbital period of the inner planet, and τ is the time of periapsis passage. For each system, we finally summarized the orbital elements in a table where we also compare the libration center value σ given by the model with the real value σ_int given by the orbital elements of the system, in order to be quantitatively aware of how far the system is from being in deep MMR.

4.2.1 System HD 73526

This system possess two confirmed exoplanets that lie almost exactly in the 2:1 MMR (Wittenmyer et al. 2014). The resonant offset is precisely 0.343%. Table 1 lists the orbital elements of the system members, and Fig. 7 shows the 𝓗₁ surface compared with a numerical integration (green curve) of the system that uses as initial condition the orbital elements of the real exoplanet.

This example is outstanding because the system is found to be almost in the deep MMR hypothesis assumed by our model.

The real critical angle value is just 4° away from the center of libration (see the footnote of Table 1). To correctness of the model is confirmed by the numerical integration, which matches one of the constant 𝓗 curves from the phase space.

This system is evolving around an ACR_I located in (∆_α,σ_α) = (0°, 0°). This ACR appears to be rather usual in the 2:1 MMR (we recall Fig. 3), and due to the topology of the contour curves, some local maxima appear in the eccentricity excursions in the proximity of ∆ϖ ± 90°.

In this case, e₁ evolves around 0.35 with a secular amplitude of approximately 0.2 (see Fig. A.1). The secular evolution of the center of libration has a certain amplitude (~ 90°), and its evolution its related to ∆ϖ and to the secular evolution of the eccentricities. This is explained, for instance, by the numerical integration in the time domain, which clearly shows that the secular frequency of all the variables is the same. This is coherent with the shape of the contour curve in which the system evolves.

Fig. 6 𝓗2(e2, ∆ϖ) contour curves in the 1:1 MMR for a system with m2/m1 = 1 and 𝒜𝓜norm = 0.9. (a) The first projection shows a near planar subsurface with σ ~ 120°. (b) The second projection shows a near planar subsurface with σ ~ 0° .In both panels, the nine numerical integrations are drawn in colors. The red crosses indicate the initial condition in each case.

4.2.2 System HD 31527

System HD 31527 is composed of three exoplanets whose orbital elements are summarized in Table 2. According to Gallardo et al. (2021), planets c and d are in a very uncommon MMR, the 16:3 MMR. In this case, we show the projection of surface 𝓗₁ in Fig. 8 since the three-dimensional representation is extremely complex.

This is an interesting case because the secular evolution of the system is really close to the a phase-space separatrix that defines whether the system is in the ACR condition. Currently, e₁ = 0.03, but in the future, it way reach values that are higher by one magnitude, regardless of whether it is in the ACR condition (see Fig. A.2).

The secular frequency relation between the orbital elements and the critical angle contrast with the system HD 73526 because in this case, for each secular period of ∆ϖ and the eccentricities, there are seven secular periods of the center of libration. This is closely related to the complexity of the 𝓗 surface (which we do not show).

We clarify that the numerical integration for this system was performed considering planets c and d alone. When we compare them with the numerical integration results shown in Gallardo et al. (2021), who considered planet b, there is a difference in the behaviour of ∆ϖ. Specifically, they obtained that ∆ϖ usually oscillates with some eventual circulations, while we obtained that it circulates in general.

Fig. 7 𝓗1 surface compared to a 10 kyr numerical integration (green curve) of the HD 73526 extrasolar system. The orbital elements are listed in Table 1.

Fig. 8 𝓗1 surface projection compared to a 100 kyr numerical integration (green curve) of the HD 31527 extrasolar system. The orbital elements are listed in Table 2. The red cross indicates the initial conditions.

5 Conclusions

The main conclusion of our work is that we successfully extended the resonant prescription in Gallardo et al. (2021) by applying the method constructed in Pons & Gallardo (2022). In that work, we studied the secular evolution of resonant motions in the framework of an elliptical restricted planar three-body problem. This model allowed us to visualize the phase space in a three-dimensional phase portrait we call 𝓗 surface. This phase portrait provided us with the secular evolution information of e₁, e₂, ∆ϖ, and σ (or θ). This has several advantages, the most important of which is the capability of encountering exhaustively all the families (or subsurfaces) of resonant evolution that can exist. It is also useful to distinguish a resonant motion with a circulating libration center from a nonresonant motion, where the critical angle circulates on resonant timescales.

We divided the study cases into two groups, the general cases, and the real cases. In the first group, we chose the most common MMRs, which are the 2:1, 3:2, 3:1, and 1:1 MMR. This last MMR, also known as co-orbital motion, is interesting because not every model is capable of describing it. We showed the two most interesting results obtained after exploring several generic systems for each MMR. This uncovered the great complexity that can exist in phase space, and also the high-eccentricity excursions that in some conditions can occur on secular timescales.

We presented two real exoplanet systems in different evolving conditions. First, system HD 73526, which is in an ACR_I condition where e₁ presents relatively low excursions, but with a moderately high absolute value. This could imply that some other mechanism has excited its eccentricity to such high values, and then, the system encountered a stable condition to continue its evolution.

The other real system we analyzed is system HD 31527, which is in a veryy rare 16:3 MMR and also very close to the strange condition in which it alternates between being in and out of an ACR. In this case, e₁ evolves between low and moderate values, whereas e₂ remains almost unchanged at a high value.

Because there is no encounter between the planets and they are in a deep MMR condition, the numerical integrations agreed with the model in all cases we studied, which is remarkable. A small correction in initial semi-major axis was sometimes required in order to be in the deep MMR condition. Another important detail that ensured that the deep MMR resonant condition (and therefore, the adiabatic invariant principle) was not broken were the initial conditions for the systems, so that they lay in a closed-contour curve on the 𝓗 surface.

Acknowledgements

The completion of this work was possible not only by the authors’ willingness and determination but also by the support of the CSIC – UdelaR since this was done in the CSIC I+D 2022 “Dinámica Secular y Resonante en Sistemas Planetarios” project context. We also want to mention the valuable corrections and comments given by Jérémy Couturier that help us significantly improve the quality of this article.

Appendix A Mathematical developments

A.1 Canonical transformation

In this section, we outline that the transformation from variables 1 to variables 4 is canonical.

There are several ways to prove whether a given transformation is canonical. We verified this by confirming whether the matrix 𝕄 associated with the transformation satisfied the symplectic condition.

This matrix the Jacobian of the new variables with respect to the old ones. Therefore, it is easy to see that $$\mathbb{M}=\left(\begin{array}{cccccccc}{k}_1& -k_2& k_2-k_1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1/k_1& 0& 0& 0\\ 0& 0& 0& 0& k_2/k_1& 1& 0& 0\\ 0& 0& 0& 0& \left(k_1-k_2\right)/k_1& 0& 1& 0\\ 0& 0& 0& 0& \left(k_1-k_2\right)/k_1& 0& 1& 1.\end{array}\right)$$(A.1)

In the context of Hamiltonian mechanics, the matrix 𝕄 satisfies the symplectic condition if (Goldstein 1987) $$\mathbb{M}\mathbb{J}{\mathbb{M}}^t=\mathbb{J},$$(A.2)

where 𝕁 is an antisymmetric matrix, defined as follows: $$\mathbb{J}=\left(\begin{array}{cc}{\mathbb{O}}_{4x4}& {\mathbb{I}}_{4x4}\\ -{\mathbb{I}}_{4x4}& {\mathbb{O}}_{4x4}.\end{array}\right)$$(A.3)

Here, 𝕆4x4 represents the null matrix, and 𝕀4x4 represents the identity matrix. Both are 4x4 matrices.

Equation A.2 states that the transformation (represented by the matrix 𝕄) preserves the symplectic structure. This means that after the transformation is performed, the canonical Hamil-tonian equations are still valid. For the sake of brevity, we do not show the detailed calculations here that prove the fulfillment of the condition A.2.

A.2 Averaged disturbing function

The third term in equation 2 was averaged to obtain the function 𝓡, the calculation of which was carried out following the approach adopted by Gallardo (2006, 2019, 2020), $$ℛ(\sigma )=\frac{1}{2\pi k_1}{\displaystyle {\int }_0^{2\pi k_1}R}\left({\lambda }_2,{\lambda }_1\left({\lambda }_2,\sigma \right)\right)d{\lambda }_2.$$(A.4)

We assumed that during k₁ revolutions of the outer planet (and therefore, k₂ revolutions of the inner planet), a_i, e_i, and ∆ϖ remain constant, which is reasonable because the variations in these elements tend to be much more slowly than those of σ (Gallardo et al. 2021). For this calculation, we set a₁ to the nominal value given by equation 8. The function R represents the instantaneous disturbing function defined in equation 3.

The relation λ₁(λ₂, σ) was obtained from Equation 9, which represents the resonant condition².

A.3 Stable equilibrium points

We demonstrate that stable solutions³ occur at the local minima of the disturbing function. To do this, we followed the idea of small oscillations that was used, for example, in Gallardo (2020).

The idea is to consider small displacements of the equilibrium points on the resonant timescale. Let $\left(I_1,\sigma \right)=\left(I_{1_0},{\sigma }_0\right)$ be an equilibrium point, that is, values that satisfy the equations 7. We define the small displacements as $S=I_1-I_{1_0}$ and s = σ − σ₀. It is easy to see that $$\begin{array}{l}\frac{dS}{dt}=\frac{d{I}_1}{dt}=-\frac{\partial \mathscr{H}}{\partial \sigma }=-{\mathscr{H}}_{\sigma }\\ \frac{ds}{dt}=\frac{d\sigma }{dt}=\frac{\partial \mathscr{H}}{\partial I_1}={\mathscr{H}}_{I_1},\end{array}$$(A.5)

where in the last equality of each equation, we introduce the compact notation for partial derivatives.

We performed a first-order expansion of the functions 𝓗_σ(I₁,σ) and ${\mathscr{H}}_{I_1}\left(I_1,\sigma \right)$. The result expressed in terms of the variables S and s in matrix form is $$\left(\begin{array}{c}\dot{S}\\ \dot{s}\end{array}\right)=\left(\begin{array}{cc}-{\mathscr{H}}_{\sigma I_1}& -{\mathscr{H}}_{\sigma \sigma }\\ {\mathscr{H}}_{I_1{I}_1}& {\mathscr{H}}_{I_1\sigma }\end{array}\right)\left(\begin{array}{c}S\\ s\end{array}\right)=Q\left(\begin{array}{c}S\\ s\end{array}\right).$$(A.6)

For the equilibrium point $\left(I_{1_0},{\sigma }_0\right)$ to be stable, the eigenvalues of the matrix in equation A.6 must be purely imaginary. Therefore, we proceeded to calculate the characteristic polynomial of Q, $$\det (Q-\lambda \mathbb{I})=\left|\begin{array}{cc}-{\mathscr{H}}_{\sigma I_1}-\lambda & -{\mathscr{H}}_{\sigma \sigma }\\ {\mathscr{H}}_{I_1{I}_1}& {\mathscr{H}}_{I_1\sigma }-\lambda \end{array}\right|=-{\mathscr{H}}_{I_1\sigma }^2+{\lambda }^2+{\mathscr{H}}_{I_1{I}_1}{\mathscr{H}}_{\sigma \sigma }.$$(A.7)

The eigenvalues are given by the roots of the characteristic polynomial. We proceeded to compute them and imposed the stability condition $${\lambda }^2=-{\mathscr{H}}_{I_1{I}_1}{\mathscr{H}}_{\sigma \sigma }+{\mathscr{H}}_{I_1\sigma }^2<0.$$(A.8)

We recall the expression for the Hamiltonian (Eq. (6)) and have $$\begin{array}{l}{\mathscr{H}}_{I_1{I}_1}=-\frac{3{\mu }_1^2{\beta }_1^3}{k_1^2}\frac{1}{I_1^4}\\ {\mathscr{H}}_{\sigma \sigma }=-{ℛ}_{\sigma \sigma }.\end{array}$$(A.9)

Therefore, the only way for λ² < 0 to hold is if 𝓡_σσ > 0. This corresponds to positive concavity in the perturbing function, in other words, the condition for a local minimum. Using this procedure, we can conclude that if 𝓡_σσ < 0 (local maximum), then the evolution of that equilibrium point is unstable.

A.4 Time evolution of real exoplanet examples

In this section we report the results of numerical integrations of the two real exoplanet systems in the time domain.

Fig. A.1 Orbital element evolution in the time domain for exoplanet system HD 73526.

Fig. A.2 Orbital element evolution in the time domain for exoplanet system HD 31527. This integration was computed considering only the planets involved in the MMR. σ librates around a long-term circulating libration center.

References

All Tables

All Figures

	Fig. 1 Averaged resonant Hamiltonian contour curves in plane (a2, σ) in 3:1 MMR for a system with m1 = m2, e1 = 0,47, e2 = 0,38, and ∆ϖ = 84°.
In the text

	Fig. 2 Normalized angular momentum function 𝒜𝓜norm(e1, e2) for 2:1 MMR with m2 = m1. (a) Three-dimensional plot of 𝒜𝓜norm(e1, e2). (b) 𝒜𝓜norm(e1, e2) contour curves with numerical values above the curves.
In the text

	Fig. 3 Phase portraits for MMR 2:1 and 𝒜𝓜norm = 0.9. (a) 𝓗2 surface for a system with m2/m1 = 0.2. The ten numerical integrations run for 10 kyr are plotted in colors (some areas of the surface were removed for clarity). (b) 𝓗1 surface for a system with m2/m1 = 1. The nine numerical integrations run for 1 kyr are shown.
In the text

	Fig. 4 𝓗1 (e1, ∆ϖ) contour curves for the 3:2 MMR for two different systems. (a) First system: m2/m1 = 1 and 𝒜𝓜norm = 0.95. (b) Second system: m2/m1 = 5 and 𝒜𝓜norm = 0.9. In both panels, the numerical integrations are drawn in color. The red crosses indicate the initial condition in each case.
In the text

	Fig. 5 Phase portraits for MMR 3:1. (a) 𝓗1 surface for a system with m2/m1 = 1 and 𝒜𝓜norm = 0.8. The ten numerical integrations run for 10 kyr are plotted in colors. (b) 𝓗1 surface for a system with m2/m1 = 0.2 and 𝒜𝓜norm = 0.99. The five numerical integrations run for 1 kyr are shown.
In the text

	Fig. 6 𝓗2(e2, ∆ϖ) contour curves in the 1:1 MMR for a system with m2/m1 = 1 and 𝒜𝓜norm = 0.9. (a) The first projection shows a near planar subsurface with σ ~ 120°. (b) The second projection shows a near planar subsurface with σ ~ 0° .In both panels, the nine numerical integrations are drawn in colors. The red crosses indicate the initial condition in each case.
In the text

	Fig. 7 𝓗1 surface compared to a 10 kyr numerical integration (green curve) of the HD 73526 extrasolar system. The orbital elements are listed in Table 1.
In the text

	Fig. 8 𝓗1 surface projection compared to a 100 kyr numerical integration (green curve) of the HD 31527 extrasolar system. The orbital elements are listed in Table 2. The red cross indicates the initial conditions.
In the text

	Fig. A.1 Orbital element evolution in the time domain for exoplanet system HD 73526.
In the text

	Fig. A.2 Orbital element evolution in the time domain for exoplanet system HD 31527. This integration was computed considering only the planets involved in the MMR. σ librates around a long-term circulating libration center.
In the text