© The Authors 2024

1 Introduction

More than 5000 extrasolar planets have been observed so far, with hundreds of multi-planet systems listed in the catalogues. However, these discoveries alone do not provide us with a deeper knowledge of planetary systems. The current limitations of the detection techniques in particular hinder our complete understanding of the formation and evolution of planetary systems, as some of the orbital parameters are still unknown in the vast majority of cases. We have focused our study on systems detected via the radial velocity (RV) method. This detection method does not provide any information on the inclinations of a multi-planet system, since neither the orbital inclinations, i, of the planets (i.e. the inclination with respect to the line of sight) nor the mutual inclination, i_mut, between their orbital planes are determined. These unknowns entail uncertainties for the actual masses of the planets (for which only a minimal value can be inferred) and the potentially three-dimensional (3D) structure of the systems. There are a few rare cases in which inclinations were determined by combining different detection methods, and there is proof that exoplanetary systems can show a highly mutually inclined architecture. The most famous example is υ Andromedae, with a mutual inclination between the orbital planes of planets c and d of ~30º (McArthur et al. 2010). There are also Kepler-108, with its two giants that are mutually inclined by 24º (Mills & Fabrycky 2017), and K2 – 266b, with two different pairs, at i_mut = 12º and 15º (Rodriguez et al. 2018).

The dynamics of mutually inclined systems is strongly affected by the Lidov–Kozai (LK) resonance (von Zeipel 1909; Lidov 1962; Kozai 1962). This secular resonance can serve as a phase-protection mechanism for such systems, even though the eccentricity and inclination of both planetary orbits vary strongly. The LK resonance has been investigated for a wide variety of problems involving different mass ratios among the bodies (for instance, see the review by Libert 2022 on the LK resonance at different scales) and in the contexts of multi-planet systems (e.g. Michtchenko et al. 2006; Libert & Henrard 2007; Migaszewski & Goździewski 2009) and multi-star systems (e.g. Innanen et al. 1997; Naoz et al. 2013). In addition to being a phase-protection mechanism, its relevance for the formation of extrasolar systems has also been highlighted. In particular, the migration through LK cycles, combined with the tidal friction, has been widely invoked to explain the formation of hot Jupiters in binary star systems (e.g. Wu & Murray 2003; Fabrycky & Tremaine 2007; Naoz et al. 2012; Petrovich 2015; Anderson et al. 2016; Ngo et al. 2016; Vick et al. 2019), while alternative scenarios based on secular planet-planet interactions have also been put forward (e.g. Petrovich & Tremaine 2016).

In this work we focus on the long-term stability of RV-detected two-planet exosystems, in particular the stability of possible 3D configurations of these systems; we do not address the question of their formation. Following the preliminary investigations of Libert & Tsiganis (2009), Volpi et al. (2019) recently studied the dynamics of RV-detected two-planet systems by focusing on their long-term stability for different values in the parameter space (i, i_mut), in order to identify 3D system parametrisations compatible with the observations and thereby constrain the unknown system parameters. The stability of the systems was studied not only by means of a chaos indicator to identify the regular regions of the parameter space, but also by using an analytical model to identify the subsets of the parameter space that are influenced by the LK resonance. They find that all ten of the RV-detected exosystems considered in their work have orbital parameters compatible with the LK resonance (in particular the eccentricity values) when considering highly inclined system configurations. In particular, long-term regular evolutions were identified either at a low mutual inclination of the two planetary orbits or at high mutual inclination, preferentially for LK-resonant evolutions, due to the significant extent of chaos along the separatrix of the resonance for many systems.

However, the Volpi et al. (2019) study was limited to systems whose planets are distant enough from the central star that the action of general relativity (GR) could be disregarded without affecting the reliability of the results. In the present work, we aim to extend this previous study to RV-detected two-planet systems with close-in planets and place constraints on the orbital inclinations, i, and mutual inclination, i_mut, that would guarantee the long-term stability of these systems. A correct description of the evolution of the system must then include GR, whose major effect on the motion of the planet is the advance of the precession of its pericenter. Given that the LK resonance acts on the pericenter as well, it is relevant to study how these two phenomena compete over the long-term evolution of the two-planet systems. As in Volpi et al. (2019), we reduced the dimensions of the problem using an analytical approach. We considered a high-order Hamiltonian power series expansion in the eccentricities and inclinations, which is reduced to four degrees by adopting the Laplace invariant plane and accurately describes the secular dynamics of two-planet systems (Libert & Henrard 2007; Libert & Sansottera 2013; Volpi et al. 2019).

Several previous works have studied the effects of GR on extrasolar systems using different approaches and for different contexts. Veras & Ford (2010) carried out extensive long-term N-body integrations to explore the dynamics of five hierarchical two-planet systems (i.e. with small semi-major axis ratios) with close-in planets, studying the full range of possible line-of-sight and relative inclinations. They concluded that GR affects the dynamics of systems with high relative inclinations. However, unlike in analytical works, providing a general picture of the dynamics for the whole phase space is particularly tricky. By analytically expanding the perturbing Hamiltonian in the ratio of the semi-major axes, Migaszewski & Goździewski (2011) investigated the dynamics of hierarchical systems, such as cir-cumbinary planets, and showed that the dynamics is qualitatively different when considering GR. By using the same approximation but limited to the octupole level, Naoz et al. (2013) discussed how GR can suppress or excite the eccentricities in triple-star systems. More recently, employing linear secular theory, Marzari & Nagasawa (2020) showed that, for coplanar two-planet systems with a close-in planet, a significant damping of the eccentricities of the inner and outer planets occurs when GR is included, which implies that systems with a chaotic behaviour caused by mutual interactions between the planets could actually have a stable evolution thanks to GR effects.

We should remark that we do not consider here the tidal effects from the host star, which act on planets extremely close to the star (a ≲ 0.08 AU typically), for several reasons. Firstly, the computation of the tidal effects relies on parameters (such as the radius of the planet and the tidal Love number) that are currently not known for most of the RV-detected exoplanets. Secondly, with these parameters fixed to arbitrary values, previous works (e.g. Migaszewski & Goździewski 2009; Veras & Ford 2010) found that the tidal effects are generally negligible with respect to the relativistic corrections. Thirdly, the perturbation due to the dynamical flattening of the star and/or of the inner planet depends on the orbital inclination of the equatorial plane of the star (see e.g. Migaszewski & Goździewski 2009), which implies that the Laplace plane is no longer a constant reference plane and a new analytical approach to reduce the number of degrees of freedom is required.

The paper is organised as follows: in Sect. 2 we briefly describe the analytical secular model used here, as well as the properties of the LK resonance. In Sect. 3, we set the parametric study and the criteria for the selection of the RV-detected systems analysed here. An in-depth study of our results regarding the extent of the LK region for the selected systems is given in Sect. 4. We draw our conclusions in Sect. 5.

2 Analytical secular approximation

The analytical expansion used in this study is an extension of the secular approximation from Volpi et al. (2019), which is briefly described in Sect. 2.1. The inclusion of GR terms is detailed in Sect. 2.2. Finally, the LK resonance is presented in Sect. 2.3.

2.1 Secular expansion of the three-body problem

We studied the dynamics of a system formed by a central star (referred to by the index 0) and two planets (indexes 1 and 2 for the inner and outer planet, respectively) that are not close to a mean-motion resonance. As a reference plane, we chose the invariant Laplace plane, which is perpendicular to the constant total angular momentum of the system. This allowed us to perform the Jacobi reduction of the nodes (Jacobi 1842). In this plane, the Hamiltonian formulation of the problem no longer depends on the longitudes of the nodes Ω₁ and Ω₂, only on their constant difference: ΔΩ = Ω₁ - Ω₂ = π. The problem is then reduced to four degrees of freedom, and for a fixed total angular momentum, the inclinations can be expressed as functions of the eccentricities. We considered the four Poincaré variables, $$\begin{array}{l}{\text{Λ}}_i={\beta }_i\sqrt{{\mu }_i{a}_i},{\xi }_i=\sqrt{2{\text{Λ}}_i}\sqrt{1-\sqrt{1-e_i^2}}\cos {\omega }_i,\\ {\lambda }_i=M_i+{\omega }_i,{\eta }_i=-\sqrt{2{\text{Λ}}_i}\sqrt{1-\sqrt{1-e_i^2}}\sin {\omega }_i,\end{array}$$(1)

where a, e, ω, and M are the semi-major axis, eccentricity, argument of the pericenter, and mean anomaly, respectively, and having set $${\mu }_i=\mathcal{G}\left(m_0+m_i\right),{\beta }_i=\frac{m_0{m}_i}{m_0+m_i}$$(2)

for i = 1,2. The parameter D₂, defined in Robutel (1995) as $$D_2=\frac{{\left({\text{Λ}}_1+{\text{Λ}}_2\right)}^2-C^2}{{\text{Λ}}_1{\text{Λ}}_2}\text{, }$$(3)

with C the norm of the total angular momentum of the system, is a variation of the angular momentum deficit (e.g. Laskar 1997) and is by definition quadratic in terms of eccentricities and inclinations.

We expanded the Hamiltonian of the three-body problem into power series of the variables L, ξ, and η and the parameter D₂, and into Fourier series of λ, at first order in the masses. As our goal is to study the secular dynamics of the system, we averaged the Hamiltonian over the fast angles, λ, and obtained the following expression: $$H\left(D_2,\xi ,\eta \right)={\displaystyle \sum _{j=0}^6{C}_{j,m,n}}{D}_2^j{\displaystyle \sum _{m+n=0}^{12-j}{\xi }^m}{\eta }^n,$$(4)

where we have fixed the maximal order in the eccentricities considered here to 12. We highlight that the Hamiltonian, Eq. (4), is expanded with respect to the eccentricities and inclinations, and not with respect to the semi-major axis ratio, α, as is the case for the quadrupole and octupole expressions¹. As a result, this analytical approach is valid for a wide range of orbital separations between the two planets, not just for hierarchical planetary systems. Moreover, it has been shown that the secular evolution of extrasolar systems that are far from a mean-motion resonance is accurately described by this Hamiltonian approach (e.g. Libert & Henrard 2007; Libert & Sansottera 2013). A further validation for highly inclined system configurations is performed in Sect. 3. Within the two-degrees-of-freedom secular Hamiltonian formulation, Eq. (4), the semi-major axes of the two planets are constant throughout the evolution. We refer to Volpi et al. (2018, 2019) for more details.

2.2 GR correction

The main effect induced by GR on the planetary motion is an advance of the pericenter. It is common practice in the secular approximation framework to consider the effects of GR on the inner planet only (e.g. Naoz 2016) when studying the LK resonance, since the argument of the pericenter of the inner planet, ω₁, is a critical angle for the LK resonance (see Sect. 2.3). As we focused on the long-term evolution of the system, we added to the purely gravitational three-body problem described by the Hamiltonian formulation, Eq. (4), the well-known expression of secular perturbation caused by GR on the pericenter of the inner planet after an averaging over the fast angles (e.g. Migaszewski & Goździewski 2009): $$H_{\text{GR}}=-\frac{3{\mu }_1^4{\beta }_1^5}{c^2{L}_1^3{G}_1},$$(5)

where c is the speed of light, and L₁ = Λ₁ and $G_1=L_1\sqrt{1-e_1^2}$ are the Delaunay variables. Since the coordinate relative to the conjugate momentum, G₁, is the argument of the pericenter of the inner planet, ω₂, it results in the following precession induced by GR on the inner planet: $${\dot{\omega }}_{\text{GR}}=\frac{\partial H_{\text{GR}}}{\partial G_1}=\frac{3{\mu }_1^4{\beta }_1^5}{c^2{L}_1^3{G}_1^2}=\frac{3{\mu }_1^{3/2}}{c^2{a}^{5/2}\left(1-e_1^2\right)}.$$(6)

The secular Hamiltonian that includes the GR correction can be rewritten in the heliocentric Poincaré variables as: $$H\left(D_2,\xi ,\eta \right)={\displaystyle \sum _{j=0}^6{C}_{j,m,n}}{D}_2^j{\displaystyle \sum _{m+n=0}^{12-j}{\xi }^m}{\eta }^n+H_{\text{GR}}$$(7)

with $$H_{\text{GR}}=-3{\mu }_1^4{\beta }_1^5{\left[c^2{\text{Λ}}_1^3\left({\text{Λ}}_1-\frac{{\xi }_1^2+{\eta }_1^2}{2}\right)\right]}^{-1}.$$(8)

2.3 LK resonance

The LK resonance has been widely discussed since the discovery of exoplanets in eccentric orbits. This secular resonance acting on 3D configurations of planetary systems is a highly effective phase-protection mechanism that ensures their stability in the long term, even though both orbits can vary widely in terms of both eccentricity and inclination. As previously said, in the invariant Laplace plane reference frame, the planetary three-body problem averaged over the short periods can be reduced to two degrees of freedom. The main dynamical features depend on the location and stability of the equilibria of the problem. As shown by Libert & Henrard (2007), there exists an equilibrium at which both eccentricities are zero, which is stable for small mutual inclinations between the two orbital planes. However, for larger mutual inclinations, the equilibrium becomes unstable and induces a large chaotic zone around it, but also generates via bifurcation two new stable equilibria around which two regular regions exist, the so-called LK resonance regions. The LK resonance is identified by the coupled variation in the eccentricity and the inclination of the inner planet as well as the libration of the argument of the pericenter of the same planet around 90° or 270° (in the Laplace plane reference frame; see e.g. Libert & Tsiganis 2009). Additional details on the LK resonance can be found in the reviews by Naoz (2016) and Libert (2021).

The critical mutual inclination corresponding to the change in stability of the equilibrium at zero eccentricities (and thus the apparition of the LK equilibria) depends on the mass and semi-major axis ratios and is typically around 40°–45° for mass ratios between 0.5 and 2 (Libert & Henrard 2007). The linear stability of the central equilibrium can be studied by following the approach adopted in, for example, Henrard & Lemaître (2005) and Libert & Henrard (2007). Around an equilibrium point (ξ^*, η^*), the quadratic approximation of the Hamiltonian, Eqs. (4) or (7), can be written as $$2Q=a_{11}{s}_1^2+2a_{12}{s}_1^2{s}_2^2+a_{22}{s}_2^2+b_{11}{S}_1^2+2b_{12}{S}_1{S}_2+b_{22}{S}_2^2,$$(9)

where $$a_{i,j}={\frac{{\partial }^{2}Q}{\partial {\xi }_i\partial {\xi }_j}|}_{\left({\xi }^{*},{\eta }^{*}\right)},b_{i,j}={\frac{{\partial }^2\mathcal{Q}}{\partial {\eta }_i\partial {\eta }_j}|}_{\left({\xi }^{*},{\eta }^{*}\right)}(i,j=1,2),$$(10)

and (s_i, S_i) are the increments relative to the variables (ξ_i, η_i). Applying the untangling transformation described in detail in Henrard & Lemaître (2005), we obtain a new formulation of the Hamiltonian, which can now be written as a linear combination of squares: $$2Q=\left(c_{11}{\tilde{s}}_1^2+d_{11}{\tilde{S}}_1^2\right)+\left(c_{22}{\tilde{s}}_2^2+d_{22}{\tilde{S}}_2^2\right),$$(11)

where $\left({\tilde{s}}_i,{\tilde{S}}_i\right)$ for i, j = 1,2 are the new variables. By studying the sign of the products c_iid_ii, we can determine the linear stability: when it is positive, the corresponding degree of freedom is linearly stable; otherwise, it is linearly unstable. The linear frequencies can be easily computed when c_iid_ii > 0 by introducing a new set of action-angle variables: $${\tilde{s}}_i=\sqrt{2R_i{R}_i^{*}}\mathrm{\hspace{0.17em}sin\hspace{0.17em}}{r}_i,{\tilde{S}}_i=\sqrt{2R_i/R_i^{*}}\mathrm{\hspace{0.17em}cos\hspace{0.17em}}{r}_i,$$(12)

where $R_i^{*}=\sqrt{d_{ii}/c_{ii}}$. We finally obtain the following simplified form of the Hamiltonian: $$Q=\sqrt{c_{11}{d}_{11}}{R}_1+\sqrt{c_{22}{d}_{22}}{R}_2.$$(13)

Thus, the evolution of the ratio between the two linear frequencies, $\sqrt{c_{11}{d}_{11}}\text{ and }\sqrt{c_{22}{d}_{22}}$, associated with the equilibrium at zero eccentricities (also called central equilibrium), when increasing the value of the mutual inclination, provides the information needed about the change in stability of the central equilibrium.

3 Parametric study

In this section, we first present the criteria used to select the systems analysed in this work as well as the parametric study we carried out. A validation of the analytical approach is then performed.

3.1 Methodology

In order to select two-planet systems with close-in planets around a single star (no binary companion) investigated here, we used the following criteria: (a) the orbital period of the inner planet is shorter than 45 days (close-in inner planet), (b) the semi-major axis of the outer planet is smaller than 10 AU (systems with significant planet-planet interactions), (c) the system is not close to a mean-motion resonance, (d) the planetary eccentricities are less than 0.65, and (e) the masses of the planets are lower than 10 M_J. The first condition comes from the observation that, for the systems considered in Volpi et al. (2019; i.e. with inner orbital periods longer than 45 days), the GR effects do not influence the dynamics. All the two-planet systems referenced on exoplanet.eu that meet these criteria are listed in Table 1 along with their orbital parameters and references. As of today, 11 extrasolar systems fulfilling all the criteria have been reported.

For the parametric study, we used the same approach as in Volpi et al. (2019). We assumed that the two orbital planes are inclined by the same angle to the line of sight; therefore, i₁ = i₂ = i. We varied both the inclination, i, of the orbital planes and the mutual inclination between the orbital planes, i_mut, since these two quantities are unknown. We set the grid for the mutual inclination, i_mut, to range from 0° to 80° in steps of 0.5° and for the orbital inclination, i, from 5° to 90° in steps of 5°. These variations affect other quantities. Firstly, as the observed planetary mass is m sin i, for each value of i considered, we had to change the mass accordingly by multiplying it by a factor of 1 / sin i. Secondly, in the general reference frame, the following relation holds, $$\cos i_{\text{mut }}=\cos i_1\cos i_2+\sin i_1\sin i_2\cos \text{\hspace{0.17em}ΔΩ}$$(14)

and it determines boundaries for the compatible values of i_mut for a given value of i, namely i_mut ≤ 2i. Having fixed the values of i and i_mut, we used Eq. (14) to determine the values of the longitudes of the nodes by setting Ω₁ = ΔΩ and Ω₂ = 0. By doing so, we obtained the complete set of initial conditions, all the other parameters being fixed to their observational values. We finally performed a change of coordinates to the Laplace plane by using the following relations: $${\text{Λ}}_1\sqrt{1-e_1^2}\mathrm{\hspace{0.17em}cos\hspace{0.17em}}{i}_{\text{L}1}+{\text{Λ}}_2\sqrt{1-e_2^2}\mathrm{\hspace{0.17em}cos\hspace{0.17em}}{i}_{\text{L}2}=C\text{, }$$(15) $${\text{Λ}}_1\sqrt{1-e_1^2}\mathrm{\hspace{0.17em}sin\hspace{0.17em}}{i}_{\text{L}1}+{\text{Λ}}_2\sqrt{1-e_2^2}\mathrm{\hspace{0.17em}sin\hspace{0.17em}}{i}_{\text{L}2}=0$$(16)

where i_L1 and i_L2 denote the orbital inclinations in the Laplace plane reference frame.

Concerning the parameters relative to the numerical integration, the fact that some of the selected systems have planets whose orbital periods are only a few days has to be taken into account. The time step must therefore be chosen such that we can perform integrations that correctly reproduce the dynamics. In order to achieve an optimal balance between precision and computational cost, we set the time step as either 20 times the orbital period of the inner planet or 1 yr, whichever was shorter. We fixed the integration time to 10⁶ yr. On one hand, said integration time is long enough to evaluate the preservation of the energy and, if any, the libration of the angle ω₁. On the other, it does not demand a high computational cost: this is a crucial factor, considering the extent of the parametric study we propose. We chose a fourth-order Runge-Kutta integrator to integrate the Hamiltonian equations associated with Eq. (7).

3.2 Validation of the approach

To validate our analytical approach, we studied the convergence of the Hamiltonian, Eq. (4), for the selected systems. The goal was to ensure the accuracy of the analytical approximation, in particular for high mutual inclinations, i_mut. In Table 2, we show the convergence au sens des astronomes (i.e. the numerical convergence of the expansion, instead of the mathematical one) computed for i = i_mut = 50º (see Poincaré 1893; Libert & Henrard 2005). Having set the values for both inclinations, the parameter D₂ can be evaluated and therefore the Hamiltonian takes the form $$H(\xi ,\eta )={\displaystyle \sum _{m+n=0}^6{\overline{\mathcal{C}}}_{m,n}}{\xi }^m{\eta }^n.$$(17)

The quantities shown in Table 2 are defined as the contributions to the Hamiltonian value of the terms (in absolute values) from order 2-12 in eccentricities and inclinations. More precisely, $$H_{2i}={\displaystyle \sum _{m+n=i}\left|{\overline{C}}_{m,n}{\xi }^m{\eta }^n\right|}.$$(18)

From the last column of Table 2, which shows the ratio between H₁₂ and H₂, we see that the numerical convergence of the expansion at high mutual inclinations is obvious for the majority of the systems. For cases that show a less clear decrease in the value of H_2i, the results concerning higher mutual inclinations should be considered with caution.

4 Results

The goal of the present work is to study how GR affects the long-term evolution of systems that harbour close-in planets, in particular when considering 3D configurations of the systems. The influence of GR on the extent of the LK region is discussed in Sect. 4.1 for the 11 selected RV-detected systems with close-in planets, while analytical estimates for the importance of the GR effects in these systems are presented in Sect. 4.2.

4.1 Influence of GR on the extent of the LK region

First, we focus here on the evolution of two specific systems we chose as examples, as they depict the two main behaviours we have encountered. In Fig. 1, we show the results obtained for the HD 147018 system, whose orbital parameters are listed in Table 1. The top panels show the long-term evolution of the system for i = 50^º and i_mut = 30^º given by our analytical approach without and with relativistic corrections (left and right panels, respectively). Similarly, the bottom panels refer to the evolutions when i = 50^º and i_mut = 50^º. At low mutual inclination, namely i_mut = 30º, the argument of the pericenter of the inner planet, ω₁, circulates (red curve); thus, the system, as expected, is not in a LK resonance. On the contrary, when i_mut = 50^º, we observe the libration of ω₁, as the system is in a LK-resonant state. When the relativistic corrections are considered (right panels), we do not observe any significant difference for HD 147018, neither concerning the LK resonance nor the general behaviour.

In the top panels of Fig. 2, we show, for the whole (i_mut, i) parameter space for which Eq. (14) can be solved, the maximal eccentricity of the inner planet, defined as $$e_{1_{\max }}=\underset{t}{\max }{e}_1(t).$$(19)

The bottom panels show instead the libration amplitude of the angle ω₁, which we computed as $${\text{libr}}_{\text{\_}}\text{ampl}\left({\omega }_1\right)\text{\hspace{0.17em}=\hspace{0.17em}}\underset{t}{\max }{\omega }_{\text{1}}(t)\text{\hspace{0.17em}−\hspace{0.17em}}\underset{t}{\min }{\omega }_{\text{1}}(t)\text{. }$$(20)

We note that the planetary masses are fixed along a horizontal line in Fig. 2 and increase when moving down along a vertical line, which means that the coefficients C_j,_m,_n of Eq. (7) have to be recomputed for each horizontal line. In the left panels, we show the results obtained with the purely gravitational formulation, Eq. (4). We remark the presence of a LK resonance region where the argument of pericenter, ω₁, librates, appearing from i_mut ≃ 40° (dark blue in the bottom panel). On the right, we present the plots for which relativistic corrections are included in the Hamiltonian formulation, Eq. (7). It is straightforward to note that there is no significant difference between the two cases; in particular, the LK resonance region is unperturbed by the additional effect.

The negligible effect of GR for HD 147018 can also be observed when studying the linear stability of the central equilibrium point, as previously detailed in Sect. 2.3. In Fig. 3 (top panels), we show the evolution of the ratio of the linear frequencies $\sqrt{c_{22}{d}_{22}/c_{11}{d}_{11}}$ relative to the central equilibrium point, having fixed i = 50°. In the left panel, we display the results when taking only gravitational forces into account, and in the right panel when also including the relativistic corrections. This shows that the critical value of the mutual inclination (i.e. when the ratio changes sign, namely i_mut ≃ 40^°) is not affected by the introduction of the relativistic effects.

We now focus on the GJ 649 system, whose orbital parameters are listed in Table 1. In Fig. 4 we show the long-term evolution of the system for i = 50° and i_mut = 30° (top panels) and i = 50^° and i_mut = 50^° (bottom panels). Again, in the left panels, we show the results obtained considering only gravitational interactions, and in the right panels those related to the model that also includes the relativistic corrections. Unlike the previous case, the differences between the left and right plots are obvious. For i_mut = 30^° (top panels), we observe a drastic change in the secular period of the evolution: for example, the secular period of the eccentricity of the inner planet decreases from ~1.1 × 10⁵ yr down to ~1.5 × 10⁴ yr. Moreover, we see that the eccentricity excitation of the inner planet is reduced with the GR contribution (see e.g. Migaszewski & Goździewski 2009; Sansottera et al. 2014; Marzari & Nagasawa 2020). For i_mut = 50^° (bottom panels), we note an additional consequence: the libration of the angle ω1 (bottom left) vanishes when considering relativistic corrections (bottom right). The introduction of GR has therefore inhibited the LK resonance. Moreover, we observe effects of GR on the eccentricity of the outer planet as well. In Fig. 5, we focus on the evolution of e₂ when considering the purely gravitational model (red line) and when including GR (blue line), for i_mut = 30° (left panel) and i_mut = 50° (rightpanel). An important variation decrease is observed when introducing the GR corrections (as previously observed in e.g. Marzari & Nagasawa 2020).

The effects of GR are even clearer when we plot the maximal eccentricity of the inner planet and the libration amplitude of ω₁ for the whole (i_mut, i) parameter space (Fig. 6). In the top-right panel, we see a stabilising effect for high mutual inclinations: the maximal eccentricity of the inner planet is much closer to the initial value (e₁ = 0.2) since the excitation is damped by GR. The LK resonance region vanishes when GR effects are included: all the initial conditions show the circulation of the angle ω₁ (bottom-right panel).

This can also be observed when studying the linear stability of the central equilibrium point. In Fig. 3 (bottom panels), we show the evolution of the ratio of the linear frequencies $\sqrt{c_{22}{d}_{22}/c_{11}{d}_{11}}$, associated with i = 50°, without considering the GR effects (left panel) and when including the relativistic corrections (right panel). In the last case, we observe no change in the linear stability of the central equilibrium around i_mut = 40^°, but the destabilisation, if any, would take place at a very high mutual inclination. This explains why no LK resonance region is observed in the bottom-right panel of Fig. 6.

In Figs. 7 and 8, we present the results for the 11 systems with close-in planets considered here, namely the maximal eccentricity and the libration amplitude of the argument of the pericenter for the inner planet, respectively, in the purely gravitational case. In Figs. 9 and 10, we display the corresponding results when taking the relativistic corrections into account. We observe that many systems experience a damping of the excitation of the inner eccentricity at high mutual inclinations, similar to the one observed for GJ 649. This is in line with the results obtained by Veras & Ford (2010), who show that the inclusion of GR in the simulations flattens the evolution of the eccentricities. Our results show that this flattening corresponds to the disappearance of the LK resonance region. As can be observed from Fig. 10, the LK region remains in place for only four systems when the relativistic effects are considered. It should be noted that the systems HD 187123, HD 217107, and HD 47186 present, in the purely gravitational case, extremely long secular periods (see Sect. 4.2). This explains the strange shape of their LK resonance region in Fig. 8, whose actual extent could be revealed by integrations over a longer time span.

Fig. 1 Dynamical evolution of the HD 147018 system. In the left panels, the evolution is given by the secular Hamiltonian, Eq. (4), and in the right panels by the Hamiltonian that includes relativistic corrections, Eq. (7). The inclination of the orbital plane is fixed to i = 50° and the mutual inclination between the planets to imut = 30° (top panels) and 50° (bottom panels).

Fig. 2 Long-term evolution of the HD 147018 system when varying the mutual inclination, imut (x-axis) and the inclination of the orbital plane, i (y-axis), both expressed in degrees. Top panels: maximal eccentricity of the inner planet, as defined by Eq. (19), without (left) and with (right) relativistic corrections. Bottom panels: libration amplitude of the argument of the pericenter, ω1 (in degrees), as defined by Eq. (20), without (left) and with (right) relativistic corrections.

Fig. 3 Evolution of the ratio of the linear frequencies $\sqrt{c_{11}{d}_{11}}\text{ and }\sqrt{c_{22}{d}_{22}}$, as defined by Eq. (13), relative to the central equilibrium point for the systems HD 147018 (top panels) and GJ 649 (bottom panels), with i fixed to 50°. In the left panels, the results only take gravitational forces into account. In the right panels, relativistic corrections are also included.

Fig. 4 Same as Fig. 1 but for the GJ 649 system.

Fig. 5 Evolution of the eccentricity of the outer planet of GJ 649 for ί = 50°, with and without relativistic effects. Left panel: imut = 30°. Right panel: imut = 50°.

4.2 Analytical estimates for the importance of the GR effects

As highlighted in the previous section, the relevance of GR in the dynamics of the selected systems varies from case to case. We were therefore interested in evaluating the relative contributions of the different effects acting on the pericenter of the inner planet. To investigate the importance of the GR effects on the inner planet, we computed two different ratios of the precession rates, following Veras & Ford (2010), namely ${\chi }_{\text{sec}}={\dot{\omega }}_{\text{sec}}/{\dot{\omega }}_{\text{GR}}$ and ${\chi }_{\text{LK}}={\dot{\omega }}_{\text{LK}}/{\dot{\omega }}_{\text{GR}}$, with ${\dot{\omega }}_{\text{sec}}$ the pericenter precession induced by the outer planet and ${\dot{\omega }}_{\text{LK}}$ the one corresponding to the LK oscillations.

The expression of the pericenter precession due to GR was previously derived in Eq. (5). The precession caused by the secular interaction of the inner planet with the outer planet can be derived using the Laplace-Lagrange secular theory (approximation to the second order in the eccentricities) and can be written as (e.g. Zhou & Sun 2003) $${\dot{\omega }}_{\text{sec}}\simeq \sqrt{{\left(y_1-y_2\right)}^2+4y_1{y}_2\left(\frac{b_{3/2}^{(2)}(\alpha )}{b_{3/2}^{(1)}(\alpha )}\right)},$$(21)

where $$\begin{array}{cc}& y_1=\frac{1}{4}\frac{m_2{a}_1^{-3/2}}{\sqrt{m_0+m_1}}{\alpha }^2{b}_{3/2}^{(1)}(\alpha )\\ & y_2=\frac{1}{4}\frac{m_1{a}_2^{-3/2}}{\sqrt{m_0+m_2}}\alpha b_{3/2}^{(1)}(\alpha ),\end{array},$$(22)

$b_{3/2}^{(j)}(\alpha )$ are the Laplace coefficients, and α = a₁/a₂ is the semi-major axis ratio.

The ratio of these two precession rates, $${\chi }_{\text{sec}}=\frac{{\dot{\omega }}_{\text{sec}}}{{\dot{\omega }}_{\text{GR}}}\text{,}$$(23)

gives us an estimation of the importance of GR for the long-term evolution of the system. The smaller the value of χ_sec, the more relevant the GR effects. In Table 3 (third column), we list the values of χ_sec for the 11 exosystems considered with i=90^○.The result shown in Figs. 7–10, are perfectly coherent with the values of χ_sec. All the systems for which χ_sec ≪ 1 are clearly affected by the introduction of the relativistic corrections, namely GJ 649, GJ 832, HD 187123, HD 190360, HD 217107, and HD 47186. This is particularly visible in Fig. 9, where the maximal eccentricity of the inner planet in these systems stays almost constant at a low level for all i_mut values. The importance of the GR corrections for the GJ 649 system is also apparent in Fig. 4. For all the other systems, the influence of the relativistic effects are more limited. No variation in the long-term evolution is reported when χ_sec ≫ 1, as previously shown in Fig. 1 for the HD 147018 system.

Regarding the highly mutually inclined systems, we can also estimate the influence of the relativistic effects on the LK resonance. The period of the LK oscillations is given by (Kiseleva etal. 1998) $$P_{\text{LK}}=\frac{2P_2^2}{3\pi P_1}\frac{m_0+m_1+m_2}{m_2}{\left(1-e_1^2\right)}^{3/2}.$$(24)

Consequently, the ratio of the two precession rates can be written as $${\chi }_{\text{LK}}=\frac{{\dot{\omega }}_{\text{LK}}}{{\dot{\omega }}_{\text{GR}}}=\frac{2\pi }{P_{\text{LK}}{\dot{\omega }}_{\text{GR}}}.$$(25)

The smaller the value of χ_LK, the more relevant the effects of GR on the long-term evolution. Table 3 (fourth column) reports the values of χ_LK for the different systems. The systems for which the LK region is clearly maintained in Fig. 10, namely GJ 682, HD 11964, HD 147018, and HD 9446, correspond well to the higher values of χ_LK in Table 3.

It should be noted that both χ_sec and χ_LK depend on the masses of the three bodies. Therefore, they should be computed for each value of the orbital plane inclination, i, as its variation determines a change in the masses of the planets as well. In Table 3, we list the precession ratios corresponding to different i values for GJ 649, which is highly influenced by GR, and HD 147018, which is not significantly affected by GR. We see that the change in the precession ratios is small, and for the two systems considered here, we observe that the balance between the effects is independent of the value of the inclination of the orbital plane.

Fig. 6 Same as Fig. 2 but for the GJ 649 system.

Fig. 7 Maximal eccentricity of the inner planet of the 11 systems listed in Table 1 when no relativistic corrections are taken into account.

Fig. 8 Libration amplitude of ω1 for the 11 systems listed in Table 1 when no relativistic corrections are taken into account.

Fig. 9 Same as Fig. 7, but with relativistic corrections included.

Fig. 10 Same as Fig. 8, but with relativistic corrections included.

5 Conclusions

We have studied the possibility that 11 RV-detected exoplanetary systems with close-in planets have a stable evolution in a 3D configuration. Given the proximity to the host star, we considered the relativistic corrections for the innermost planet and determined how the GR effects influence the purely gravitational evolution. To do so, we adopted an analytical approach in a Hamiltonian formalism, where a term for the GR effects on the pericenter was added to the classical expansion in the eccentricities and inclinations. We find that most of the considered systems show LK resonance regions when we only take the gravitational interactions between the bodies into account, in line with the results obtained in previous studies (see Libert & Tsiganis 2009; Volpi et al. 2019). However, in the majority of the cases, since the LK resonance and GR both act on the evolution of the argument of the pericenter of the inner planet, the LK resonance region disappears as we introduce the relativistic corrections. These observations are in agreement with analytical estimates for the pericenter precession ratios, highlighting the validity of our secular Hamiltonian approach.

Thanks to this study, we can now address the question of the possible architectures of RV-detected exosystems with close-in planets. Our study indicates that there are two categories of systems. On the one hand, the long-term stability of the GJ 649, GJ 832, HD 187123, HD 190360, HD 217107, and HD 47186 systems is possible whatever the mutual inclination of the two planetary orbits, since GR strongly affects the dynamical evolution of the planets and dominates over the LK resonance. On the other hand, the GJ 682, HD 11964, HD 147018, and HD 9446 exosystems are not compatible with all the mutual inclinations between the two orbital planes, only either at a low mutual inclination of the orbits or in the LK regime at high mutual inclination. This is due to the rapid destabilisation of highly mutually inclined orbits resulting from the significant chaos that develops around the stability islands of the LK resonance (see Volpi et al. 2019, for more details). The analytical estimates described in this work can easily be used to determine which category a given system belongs to.

This work is open to future developments. Firstly, a careful study of the stability of the systems by means of a chaos indicator would be complementary to the present study. It would provide a more complete panoramic view of the stability of the systems and precisely identify all the initial conditions that guarantee long-term stability. Secondly, taking tidal effects into account, which we have not done here, would provide a more precise description of the dynamics of the systems. However, as previously noted, this is not straightforward and is beyond the scope of this work.

Acknowledgements

M.V. acknowledges financial support from the FRIA fellowship (F.R.S.-FNRS). The work of A.-S.L. is supported by the Fonds de la Recherche Scientifique – FNRS under Grant No. F.4523.20 (DYNAMITE MIS-project). Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11.

References

All Tables

All Figures

	Fig. 1 Dynamical evolution of the HD 147018 system. In the left panels, the evolution is given by the secular Hamiltonian, Eq. (4), and in the right panels by the Hamiltonian that includes relativistic corrections, Eq. (7). The inclination of the orbital plane is fixed to i = 50° and the mutual inclination between the planets to imut = 30° (top panels) and 50° (bottom panels).
In the text

Fig. 2 Long-term evolution of the HD 147018 system when varying the mutual inclination, imut (x-axis) and the inclination of the orbital plane, i (y-axis), both expressed in degrees. Top panels: maximal eccentricity of the inner planet, as defined by Eq. (19), without (left) and with (right) relativistic corrections. Bottom panels: libration amplitude of the argument of the pericenter, ω1 (in degrees), as defined by Eq. (20), without (left) and with (right) relativistic corrections.

In the text

Fig. 3 Evolution of the ratio of the linear frequencies $\sqrt{c_{11}{d}_{11}}\text{ and }\sqrt{c_{22}{d}_{22}}$, as defined by Eq. (13), relative to the central equilibrium point for the systems HD 147018 (top panels) and GJ 649 (bottom panels), with i fixed to 50°. In the left panels, the results only take gravitational forces into account. In the right panels, relativistic corrections are also included.

In the text

	Fig. 4 Same as Fig. 1 but for the GJ 649 system.
In the text

	Fig. 5 Evolution of the eccentricity of the outer planet of GJ 649 for ί = 50°, with and without relativistic effects. Left panel: imut = 30°. Right panel: imut = 50°.
In the text

	Fig. 6 Same as Fig. 2 but for the GJ 649 system.
In the text

	Fig. 7 Maximal eccentricity of the inner planet of the 11 systems listed in Table 1 when no relativistic corrections are taken into account.
In the text

	Fig. 8 Libration amplitude of ω1 for the 11 systems listed in Table 1 when no relativistic corrections are taken into account.
In the text

	Fig. 9 Same as Fig. 7, but with relativistic corrections included.
In the text

	Fig. 10 Same as Fig. 8, but with relativistic corrections included.
In the text