# Secular Orbit Evolution in Systems with a Strong External Perturber-A Simple and Accurate Model

Eduardo Andrade-Ines, Siegfried Eggl

► To cite this version:

Eduardo Andrade-Ines, Siegfried Eggl. Secular Orbit Evolution in Systems with a Strong External Perturber-A Simple and Accurate Model. The Astronomical Journal, 2017, 153 (4), pp.148. 10.3847/1538-3881/153/4/148 . hal-02202118

HAL Id: hal-02202118

https://hal.science/hal-02202118v1

Submitted on 5 Sep 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Secular Orbit Evolution in Systems with a Strong External Perturber—
A Simple and Accurate Model

Eduardo Andrade-Ines¹ and Siegfried Eggl²

¹Institute de Mécanique Céleste et des Calcul des Éphémérides—Observatoire de Paris,
77 Avenue Denfert Rochereau, F-75014 Paris, France; eandrade.ines@gmail.com

²Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, 91109 Pasadena, CA, USA; siegfried.eggl@jpl.nasa.gov

Received 2016 November 14; revised 2017 January 4; accepted 2017 January 12; published 2017 March 9

## Abstract

We present a semi-analytical correction to the seminal solution for the secular motion of a planet's orbit under gravitational influence of an external perturber derived by Heppenheimer. A comparison between analytical predictions and numerical simulations allows us to determine corrective factors for the secular frequency and forced eccentricity in the coplanar restricted three-body problem. The correction is given in the form of a polynomial function of the system's parameters that can be applied to first-order forced eccentricity and secular frequency estimates. The resulting secular equations are simple, straight forward to use, and improve the fidelity of Heppenheimers solution well beyond higher-order models. The quality and convergence of the corrected secular equations are tested for a wide range of parameters and limits of its applicability are given.

**Key words:** binaries: general – celestial mechanics – planets and satellites: dynamical evolution and stability – planets and satellites: formation

## 1. Introduction

Understanding the gravitational three-body problem is one of the fundamental challenges in celestial mechanics. More than three centuries of study have produced a wide variety of methods to tackle the dynamics of three gravitating point masses (Marchal 1990). Although general analytical solutions do exist, they tend to be unwieldy due to the complexity of the problem (e.g., Sundman 1907). Thus, simplified approaches tailored to specific needs are still popular. In order to describe the long-term evolution of hierarchical three-body systems, for instance, that fast-evolving variables can be eliminated from the equations of motion either through averaging or via canonical transformations of the corresponding Hamiltonian system (Ferraz-Mello 2007). The resulting “secular equations” are less cumbersome to deal with and permit simple analytical solutions to the dynamical evolution of the system. In his seminal work on planet formation in binary stars, Heppenheimer (1978) described the orbital evolution of a satellite influenced by an external perturber, e.g., a planet orbiting one component of a binary star considering all bodies to be point masses.<fnref n="3" /> Assuming the planet started out on a circular orbit in the same orbital plane as the massive bodies, its eccentricity ($e_1$) and longitude of pericenter ($\varpi_1$) evolve as

$$e_1(t) = \epsilon [4 - 4 \cos(\text{gt})]^{1/4}, \quad \varpi_1(t) = \frac{g}{2} t, \quad (1)$$

where $\epsilon$ and $g$ are the so-called *forced eccentricity* and *secular frequency*, respectively. Both $\epsilon$ and $g$ are constants that depend on the system parameters as shown in Section 2. As in Heppenheimer's approximation, the planet and the perturber only exchange angular momentum; specific orbital energies, semimajor axes, and periods stay constant for all bodies. Hence, Equation (1) completely describes the secular evolution of the system.<fnref n="4" /> The

nonlinear and singular aspects of the three-body problem, however, cause simplified solutions such as the one above to sometimes fail to capture the actual evolution of the system (Andrade-Ines et al. 2016).

Today, direct numerical integration of the equations of motion is a viable option. While purely numerical methods have their own shortcomings, they provide relatively accurate solutions on short to intermediate timescales (e.g., Eggl & Dvorak 2010, p. 431). Yet, for large-scale planet formation simulations in particular, the time requirements to individually propagate each particle through numerical integrations are still prohibitive (Thebault 2011; Leiva et al. 2013). A natural way to tackle such issues is to apply more complete analytic models, but those often contain higher-order expansions leading to a substantial increase in complexity (Georgakarakos 2003; Andrade-Ines et al. 2016). Alternatively, one can try to modify lower order theories so as to better fit the behavior of the actual system under investigation (Giuppone et al. 2011).

In this work we take the latter approach and derive an empirical correction to the classical solution of Heppenheimer (1978) for the forced eccentricity and secular frequency. Our corrections to Equation (1) presented in this work improve the fidelity of Heppenheimer's solution substantially while still retaining its simple and elegant form. The rest of the article is structured as follows. In Section 2, we present a more detailed derivation of Heppenheimer's equations and discuss recent developments and extensions to the first-order secular approach. Section 3 contains information on how we derive the corrections to the forced eccentricity and secular frequency from the integration of the exact equations of motion. Results for the corrections stemming from least squares fits for both the secular frequency and the forced eccentricity are presented in Section 4, and in Section 5, we discuss some applications of the new model. Concluding remarks are provided in Section 6.

<footnote n="3">³ That is, neglecting the dynamical effects due to the shapes of the bodies.</footnote>

<footnote n="4">⁴ This is assuming that the planet does not influence the orbit of the perturber and that both the planet's and the perturber's angular momentum vectors are parallel.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Representation of the system consisted of a central body of mass $m_0$ being orbited by a satellite of (negligible) mass $m_1$. Another body of mass $m_2$ revolves around the primary so that its orbit never crosses the path of the satellite.</figcaption></figure>

## 2. Analytical Secular Orbit Evolution Models

Let us consider a system consisting of a main body with mass $m_0$, a planet with mass $m_1$, and a secondary body with mass $m_2$ in a reference frame centered on $m_0$. Let $\mathbf{r}_i$ be the position vector of the $i$th body with respect to $m_0$ (see Figure 1). We assume, furthermore, that all bodies move in the same plane. If $m_1 \ll m_0 \sim m_2$ and $|\mathbf{r}_i| < |\mathbf{r}_2|$ for all time and neglecting gravitational effects due to $m_1$, we arrive at the coplanar restricted three-body problem of S-type (Dvorak 1986). The orbit of $m_2$ around $m_0$ will be regarded as a fixed Keplerian ellipse, whereas the orbit of the planet evolves with time. This section is dedicated to discussing previous analytical results that model the secular behavior of the host-planet-perturber configuration. Models by Heppenheimer (1978), Marchal (1990), and Andrade-Ines et al. (2016) have been selected for this purpose, as they represent various degrees of complexity and fidelity with respect to the numerical reference solutions. The here presented analytic models will serve as a benchmark for the empirical correction derived in Section 4.

### 2.1. Heppenheimer (1978)

Assuming there are no significant mean-motion resonances or high-amplitude short-period oscillations, the dynamical evolution of the planets' orbit ($m_1$) will be dominated by the secular interaction with the perturber ($m_2$). The corresponding disturbing function can be expanded in terms of the semimajor axis ratio by means of Legendre polynomials ($P_i$). While this kind of expansion limits our model to hierarchical systems, it allows us to obtain finite expressions for arbitrary eccentricities (Kaula 1962; Laskar & Boué 2010). By limiting the expansion in Legendre polynomials to $P_3$ (quadrupole problem), truncating the perturbation to low values of the eccentricity of the planet and performing a first-order averaging, Heppenheimer (1978) obtained the following secular disturbing function in orbital elements:

$$ \mathcal{R}_H = \frac{Gm_2}{(1 - e_2^2)^{3/2}} \frac{a_1^2}{a_2^3} \left[ \frac{1}{4} + \frac{3}{8} e_1^2 - \frac{15}{16} \frac{a_1}{a_2} \frac{e_1 e_2}{(1 - e_2^2)} \cos(\Delta\varpi) \right], \quad (2) $$

where $G$ is the gravitational constant, $\Delta\varpi = \varpi_1 - \varpi_2$, $a_i$, $e_i$, and $\varpi_i$ are the semimajor axis, the eccentricity, and the longitude of the pericenter of the $i$th body, respectively. Neglecting all constant terms, as they will not change the system's equations of motion, Equation (2) can be expressed in

a simpler form,

$$ \mathcal{R}_H = n_1 a_1^2 g_H \left[ \frac{e_1^2}{2} - \epsilon_H e_1 \cos \Delta\varpi \right], \quad (3) $$

where

$$ g_H = \frac{3}{4} \mu \alpha^3 \frac{n_1}{(1 - e_2^2)^{3/2}} \quad (4) $$

is the secular frequency and

$$ \epsilon_H = \frac{5}{4} \alpha \frac{e_2}{(1 - e_2^2)} \quad (5) $$

is the forced eccentricity, with $n_1 = \sqrt{Gm_0/a_1^3}$ being the mean motion of the planet, $\alpha = a_1/a_2$, and $\mu = m_2/m_0$. Introducing the variables

$$ \begin{aligned} k &= e_1 \cos \Delta\varpi, \\ h &= e_1 \sin \Delta\varpi, \end{aligned} \quad (6) $$

the Lagrange–Laplace planetary equations up to order $O(e_1^2)$, which determine the satellite's orbit evolution, read (Brouwer & Clemence 1961):

$$ \begin{aligned} \frac{dh}{dt} &= \frac{1}{n_1 a_1^2} \frac{\partial \mathcal{R}}{\partial k} = g_H(k - \epsilon_H), \\ \frac{dk}{dt} &= - \frac{1}{n_1 a_1^2} \frac{\partial \mathcal{R}}{\partial h} = -g_H h. \end{aligned} \quad (7) $$

The system of Equation (7) admits the analytical solution

$$ \begin{aligned} k &= e_p \cos(g_H t + \phi) + \epsilon_H, \\ h &= e_p \sin(g_H t + \phi), \end{aligned} \quad (8) $$

where $e_p$ and $\phi$ are constants of integration, obtained from the initial conditions $e_1(0)$ and $\varpi_1(0)$ as

$$ e_p = \sqrt{(k_0 - \epsilon_H)^2 + h_0^2}, \quad (9) $$

$$ \cos \phi = \frac{k_0 - \epsilon_H}{e_p}, \quad \sin \phi = \frac{h_0}{e_p}. \quad (10) $$

Here, $k_0 = e_1(0) \cos \Delta\varpi(0)$ and $h_0 = e_1(0) \sin \Delta\varpi(0)$, with $\Delta\varpi(0) = \varpi_1(0) - \varpi_2$. Figure 2 shows an example of how the time evolution of the secular variables $k$ and $h$ translates into the elliptic osculating elements $e_1$ and $\Delta\varpi$. Three different sets of initial conditions were chosen for the planet's orbit, $e_1 \in \{0.001, 0.03, 0.12\}$, sharing otherwise similar system parameters $\alpha = 0.1, \mu = 1$, and $e_2 = 0.3$. The evolution of the elliptic elements $e_1$ and $\Delta\varpi$ is qualitatively different for each set of initial conditions. Yet, the very same dynamical system can be described by simple circles in the $(k, h)$ plane. While the circles have different radii depending on the initial conditions, they are all centered on $(k, h) = (\epsilon_H, 0)$. This behavior is due to the form of Equation (8) that define circles in the $(k, h)$ plane. The planet's orbit evolves along those circles with the frequency $g_H$. Additionally, we can see from Equation (4) and (5) that $g_H$ and $\epsilon_H$ are functions of the parameters of the problem $n_1$, $\alpha = a_1/a_2$, $\mu = m_2/m_0$, and $e_2$ only. Those are constants of motion in the secular restricted three-body problem. Hence, for a fixed set of parameters, the secular orbits will always have the same forced eccentricity and the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Secular orbits obtained from Heppenheimer's model for a system with parameters $\mu = 1$, $\alpha = 0.1$, and $e_2 = 0.3$, integrated from orbits initially at the pericenter ($\Delta\varpi = 0$) for three different initial planetary eccentricities: $e_1 = 0.001$ (blue), $e_1 = 0.03$ (red), and $e_1 = 0.12$ (green). Note that for all three systems $\epsilon_H = 0.041$. All circles in the (k, h) plane (bottom-right graph) are centered on $(\epsilon_H, 0)$. Each of the three orbits evolves along their respective circle and complete one revolution in 36.5 years corresponding to a secular frequency of $g_H = 0.172$ rad yr$^{-1}$ regardless of their radius.</figcaption></figure>

same secular frequency, no matter the initial values of $\varpi_1$ and $e_1$.<fnref n="5" /> As the planets eccentricity increases, higher-order terms become more important and the secular frequency may change from the predicted value. However, this is only expected to happen for $e_1 > 0.2$ (Andrade-Ines et al. 2016). If the planet's orbit was initially circular, i.e., $h(0) = k(0) = 0$, so that $e_p = \epsilon$ and $\phi = \pi$, then Equation (8) reduce to Equation (1) presented in the introduction.

## 2.2. Marchal (1990)

Using the Von Zeipel averaging method on the three-body disturbing function, Marchal (1990) obtained a second-order<fnref n="6" /> solution for the secular orbit evolution of an S-type hierarchical triple system. In the following we briefly discuss Marchal's model and the resulting equations. Still considering the planar problem, the disturbing function of Marchal's model is given by

$$ \mathcal{R}_M = Q_1 + Q_2 + Q_3, \quad (11) $$

where

$$ Q_1 = \frac{1}{4} G \frac{m_0 m_2}{m_0 + m_1} \frac{a_1^2}{a_2^3} \left( 1 + \frac{3}{2} e_1^2 \right) \frac{1}{(1 - e_2^2)^{3/2}}, \quad (12) $$

$$ Q_2 = -\frac{15}{16} G \frac{m_0 m_2 (m_0 - m_1) a_1^3}{(m_0 + m_1)^2 a_2^4} \times e_1 \left(1 + \frac{3}{4} e_1^2\right) \frac{e_2}{(1 - e_2^2)^{5/2}} \cos \Delta\varpi \quad (13) $$

originate in the Legendre polynomial expansion up to $P_3$ and

$$ Q_3 = \frac{15}{64} G \frac{m_0 m_2^2}{(m_0 + m_1)^{3/2} (m_0 + m_1 + m_2)^{1/2}} \times e_1^2 (1 - e_1^2)^{1/2} \frac{[5(3 + 2e_2^2) + 3e_2^2 \cos(2\Delta\varpi)]}{(1 - e_2^2)^3} \quad (14) $$

arises from the Von Zeipel's averaging method as a second-order term. Assuming now the restricted problem ($m_1/m_0 \to 0$) and neglecting terms of order $O(e_1^3)$, Equation (11) becomes

$$ \mathcal{R}_M = n_1 a_1^2 g_M \left[ \frac{e_1^2}{2} - \epsilon_M e_1 \cos \Delta\varpi \right], \quad (15) $$

where

$$ g_M = g_H \left[ 1 + \frac{25}{8} \frac{\mu}{\sqrt{1+\mu}} \alpha^{3/2} \frac{3+2e_2^2}{(1-e_2^2)^{3/2}} \right], \quad (16) $$

$$ \epsilon_M = \epsilon_H \left[ 1 + \frac{25}{8} \frac{\mu}{\sqrt{1+\mu}} \alpha^{3/2} \frac{3+2e_2^2}{(1-e_2^2)^{3/2}} \right]^{-1}. \quad (17) $$

The disturbing function $\mathcal{R}_M$ in Equation (15) has a form similar to $\mathcal{R}_H$ in Equation (3) where $g_M$ and $\epsilon_M$ are the quantities corresponding to $g_H$ and $\epsilon_H$, respectively. Therefore, the same general form of solutions given by Equation (8) can be expected to work. Note that, although of similar form, the expressions for the secular frequency and forced eccentricity are not identical in the two models. Equations (16) and (17) relate Marchal's secular quantities with those in Heppenheimer's model. Only if $\mu \ll 1$ (i.e., $m_2 \ll m_0$) do both models share the same analytic expressions.

## 2.3. Andrade-Ines et al. (2016)

Utilizing Hori's averaging theory (Hori 1966; Ferraz-Mello 2007), Andrade-Ines et al. (2016) recently constructed a second-order<fnref n="7" /> model that, too, provided more accurate results than Heppenheimer's solution. Relying on an expansion of the disturbing function in powers of the semimajor axis ratio and eccentricities, however, the model is complex and a complete analytical form is not available at this time. The model developed by Andrade-Ines et al. (2016) works as follows. The second-order Hamiltonian ($\mathcal{H}_A$) and the secular frequencies $g_A$ are calculated from analytical expressions presented in Andrade-Ines et al. (2016) using tabulated numerical coefficients of the disturbing function. The forced eccentricity ($\epsilon_A$) is a root of the Equation

$$ \left. \frac{\partial \mathcal{H}_A}{\partial e_1} \right|_{e_1 = \epsilon_A, \Delta\varpi = 0} = 0. \quad (18) $$

<footnote n="5">⁵ This statement holds for models up to order $O(e_1^2)$.</footnote>

<footnote n="6">⁶ In $\mu$ and $\alpha$, see also Georgakarakos (2003, 2005).</footnote>

<footnote n="7">⁷ in $\alpha^2 \cdot \mu$.</footnote>
Equation (18) can be solved for $\epsilon_A$ numerically by applying, for instance, the geometric method presented in Michtchenko & Malhotra (2004). Finally, the secular orbits for the second-order Hamiltonian can be obtained from Equation (8), replacing $g_H$ and $\epsilon_H$ by the numerical values obtained for $g_A$ and $\epsilon_A$, respectively.

### 3. Extracting Secular Solutions from Direct N-body Simulations

So far, we have seen that the secular problem we are investigating can be reduced to one degree of freedom associated with the constant frequency $g$. As $g$ only depends on fixed parameters of the problem ($\alpha$, $\mu$, and $e_2$), it can be calculated in advance. Returning to Figure 2 we see that the secular evolution of the planet's orbit in the $(k, h)$ plane can be described with circles centered at $(\epsilon, 0)$. The radius of such a circle is equal to $e_p$, which, similar to $\phi$, depends exclusively on the system's initial conditions. In particular for $e_p = 0$ we have a stationary secular solution, with $e_1(t) = \epsilon$ and $\Delta\varpi(t) = 0$. In other words, the planet's orbit ($a_1, e_1, \varpi_1$) does not change its shape and keeps its prior orientation with respect to the perturber's orbit, ($\dot{a}_1 = \dot{e}_1 = \dot{\varpi}_1 = 0$). In the complete<fnref n="8" /> planar problem, however, we have three degrees of freedom, associated with the following three frequencies: the mean motions of the planet ($n_1$) and the perturber ($n_2$) and the secular frequency $g$. In this case, the stationary secular solution corresponds to a quasi-periodic orbit in osculating elements with frequencies $n_1$ and $n_2$ and zero amplitude in all quantities associated with the secular frequency $g$. Therefore, the orbits and the initial conditions of the complete problem can be fundamentally different from those of the secular (averaged) problem, which makes the comparison between the two models a non-trivial task. Methods based on the frequency analysis (Laskar 1990; Michtchenko et al. 2002) of numerically integrated orbits have been used successfully by many authors to determine such quasi-periodic solutions (Noyelles et al. 2008; Couetdic et al. 2010; Andrade-Ines et al. 2016). In this work, we chose to apply the method described in Andrade-Ines et al. (2016) that was developed for the case of stationary secular solutions in the restricted three-body problem. For a fixed set of parameters, this method finds the quasi-periodic secular solution after successive iterations of a digital filter that eliminates the secular component of the time series of the variable $\xi = e_1 e^{i\Delta\varpi}$. Once the initial conditions of the quasi-periodic orbit are determined, the secular frequency $g_E$ is obtained via harmonic decomposition, while the forced eccentricity is calculated as

$$ \epsilon_E = \langle e_1(t) \cos \Delta\varpi(t) \rangle_T, \quad (19) $$

where $T$ is the total time of integration. Unless the secular period $P_E = 2\pi/g_E$ is known exactly, this process yields the true values for $\epsilon_E$ only if the total time of integration tends to infinity. However, good approximations can be found for sufficiently long<fnref n="9" /> times of integration. Fortunately, as the secular component was already eliminated from the quasi-periodic solution, the only periodic components still present are the fast frequencies $n_1$ and $n_2$. Therefore, by considering a time of integration of that is much longer than the Keplerian period of the outer body, the error in the determination of $\epsilon_E$ is going to be small as well. In our calculations, we considered a time of

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Forced eccentricity (top) and secular frequency (bottom) estimates as a function of the semimajor axis ratio for systems sharing the parameters $m_0 = 1 M_\odot$, $m_1 = 10^{-5} M_\odot$, $a_2 = 1$ au, $\mu = 1$, and $e_2 = 0.2$. Only the satellite's semimajor axis $a_1$ is varied. Results are shown for the following models: Heppenheimer (1978; red curve), Marchal (1990; blue curve), Andrade-Ines et al. (2016; green curve) and Giuppone et al. (2011; magenta curve) models. The numerically determined solution can be considered the most accurate one (black dots). The vertical dashed lines represent the nominal positions of the 7:1 (orange) and 8:1 (magenta) mean-motion resonances. The numerical solution is omitted for $\alpha < 0.025$, as this is a region where the external perturbations tend to zero and all models start to coincide.</figcaption></figure>

integration of $T = 12\pi/g_H$, that is six times the secular period of the system. Throughout this work we shall consider the solutions for $\epsilon_E$ and $g_E$ obtained from direct N-body simulations to be the most accurate ones. They will serve as calibration values for the results originating from analytic and semi-analytic methods. Though fast and accurate, the exclusively numerical methodology to calculate secular quantities has its limitations as well. For instance, if the orbits are chaotic or in a mean-motion resonance (MMR) the method will not converge (Andrade-Ines et al. 2016). Similarly, care has to be taken when comparing initial conditions used for the same orbit in both the complete problem and the models. A given set of osculating initial conditions of the complete problem is not the same as the averaged initial conditions of the secular models, especially for the case when there are strong short-period perturbations involved (Andrade-Ines et al. 2016). In this work, we consider as initial conditions of the secular problem $(k(0), h(0))$ the average of the same variables $(k(t), h(t))$ over one orbital period of the outer body calculated with the complete problem.

#### 3.1. Results of the Comparison

In order to check the fidelity of the methods presented in the previous sections, we calculated the stationary secular solutions for several systems using each of the previously mentioned models. Figure 3 contains a comparison among all models, where the fully numerically computed secular results serve as a reference. The model by Heppenheimer (1978) is referred to as Hepp, Marchal (1990)'s model as March, and the one by

<footnote n="8">⁸ Non-secular, i.e., not averaged over the fast angles.</footnote>

<footnote n="9">⁹ At least several secular periods.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Relative error of the secular frequency *g* calculated using *Hepp* (left column), *March* (middle column), and *Corr* (right column) with respect to the reference solution as a function of the semimajor axis ratio (α). The three rows represent different perturber eccentricity values (e₂) and the color code shows results for various mass ratios of the binary (μ). Similar to Figure 3, values for α < 0.025 were not calculated as the perturbations become negligible and all methods start to provide accurate approximations to the planetary orbit.</figcaption></figure>

Andrade-Ines et al. (2016) as A-I. In the same figure, we include the results given by the empirical correction derived by Giuppone et al. (2011; referred as Giup) designed for the γ Cephei system. The same abbreviations remain valid throughout the remainder of this work. One can see that, for small values of the semimajor axis ratio α, all of the models coincide with the reference solution. For large α values, however, none of the models provide accurate representations of the system dynamics. It becomes evident from Figure 3, though, that more complex models of higher order tend to work better for larger values of α. Although Giup is an empirical approximation to a second-order analytical model, it provides a poor fit to the forced eccentricity for the particular case presented in Figure 3. This is understandable, as the approximation was derived to study the γ Cephei system and shall be only applicable to systems sharing similar parameters. The gap in the curve of the reference solution at α ≈ 0.2 is due to the 7:1 mean-motion resonance (MMR) between the satellite and the perturber. As discussed in the previous section, crossing an MMR impedes the convergence of the iterative method used to construct the reference solution. Such gaps are going to appear whenever MMRs begin to be significant for the dynamical evolution.

**4. Introduction of Empirical Corrections**

A comparison between Equations (3) and (15) shows that by
limiting expansions to low eccentricities of the satellite we can

write the Hamiltonian in a quadratic form permitting a general
solution akin to Equation (8). Since *Hepp* and *March* models
are not identical, differences in the actual values of the secular
frequency and the forced eccentricity are to be expected. A
glance at Equations (16) and (17) reveals that the higher-order
terms can be interpreted as corrections to the first-order
solutions for the forced eccentricity *ϵ*<sub>H</sub> and the secular
frequency *g*<sub>H</sub>. The second-order solution can take a very
complex form, though, and worse, it may still not describe the
secular motion for all stable orbits in the relevant parameter-
space with the desired accuracy (Andrade-Ines et al. 2016).
Figure 3 contains some of the systems where the second-order
model does not match the reference solution at all, especially
for large semimajor axis ratios. Those systems may, perhaps,
be properly described by third- or higher-order solutions,
which, unfortunately, further increase in complexity.

In order to circumvent this issue, one can take an alternative
approach by introducing *empirical* corrections. Those are
coefficients added to a simple model, which are constructed
from the reference solution through a fitting procedure. This
was successfully done by Giuppone et al. (2011), where the
authors presented corrected values of *g* and *ϵ*. Those, however,
were only valid for the particular case of the γ-Cephei system.
In this work, we will present a general expression for these
quantities, valid for a much wider range of parameters. Let us
define the relative difference (error) between the reference
solution and Heppenheimer's model in the secular frequency
<figure><img src="image_6.png" /><figcaption>Figure 5. Same as Figure 4 but for the forced eccentricity $\epsilon$.</figcaption></figure>

$g_H$ and the forced eccentricity $\epsilon_H$, respectively, as

$$ \delta_g(\alpha, e_2, \mu) = 1 - \frac{g_E}{g_H}, \qquad (20) $$

$$ \delta_e(\alpha, e_2, \mu) = 1 - \frac{\epsilon_E}{\epsilon_H}, \qquad (21) $$

where $g_E$ and $\epsilon_E$ are obtained numerically. Those errors are evaluated for several mass ratios of the binary $\mu = (0.1, 0.2, 0.5, 1, 2, 5, 10)$, given a set of binary eccentricities $e_2 = (0.1, 0.2, 0.3, 0.4, 0.5, 0.6)$ and semimajor axis ratios $\alpha \in [0.01, 0.4]$. The sampling interval of the latter was varied between 0.01 to 0.05 depending on the system, such that there were at least 25 different values of $\alpha$ for each pair of parameters ($\mu, e_2$). For all simulations, we fixed the mass of the central body at $m_0 = 1M_\odot$ and the semimajor axis of the perturber at $a_2 = 1$ au. We then applied the method of least squares to minimize the error functions $\delta_g$ and $\delta_e$, assuming that the correction functions have a polynomial dependence on the parameters of the problem $\mu = m_2/m_0$, $\alpha = a_1/a_2$ and $e_2$, in the form of

$$ \delta_g = \sum_{i=1}^{N_g} A_i \alpha^{p_i} e_2^{q_i} \mu^{l_i}, \qquad (22) $$

$$ \delta_e = \sum_{i=1}^{N_e} A'_i \alpha^{p'_i} e_2^{q'_i} \mu^{l'_i}. \qquad (23) $$

Naturally, one can expect that the larger the number of terms in the fit, the smaller the total error. However, the higher the order of the fitting polynomials, the faster the fit deviates from the nominal solution outside of the fitting region. Furthermore, high-order polynomials quickly become unwieldy themselves, defeating the purpose of simplifying the secular model. Limiting our expressions to less than 20 terms, a series of different functions $\delta_e$ and $\delta_g$ were tested and we applied an iterative method to determine the optimum number of terms and the order ($p_i, q_i, l_i$) of the polynomial functions that resulted in the smallest error for the fitted quantities. The results are the following expressions:

$$ \begin{aligned} \delta_g = & \alpha^{3/2} [-4.6274\mu^{1/2} - 4.0190\mu + 0.25041\mu^2 \\ & - 3.41e_2^2\mu^{1/2} + 11.09e_2^2\mu - 0.9823e_2^2\mu^2 \\ & - 20.13e_2^4\mu^{1/2} - 85.49e_2^4\mu + 4.996e_2^4\mu^2] \\ & + \alpha^{9/2} [123.67\mu^{1/2} - 799.20\mu - 201.49\mu^2 \\ & + 180e_2^2\mu^{1/2} - 5555e_2^2\mu - 617.7e_2^2\mu^2 \\ & + 2.671 \times 10^4 e_2^4\mu^{1/2} - 1.0229 \times 10^5 e_2^4\mu \\ & - 23076e_2^4\mu^2], \end{aligned} \qquad (24) $$
$$
\begin{equation}
\begin{aligned}
\delta_e ={}& \alpha^{3/2}[29.494e_2\mu^{1/2} + 9.220e_2\mu - 99.85e_2^2\mu^{1/2} \\
& \qquad - 31.50e_2^2\mu + 124.60e_2^3\mu^{1/2} + 35.69e_2^3\mu] \\
& + \alpha^{9/2}[1073.0e_2\mu^{1/2} + 4280e_2\mu - 1609.8e_2\mu^2 \\
& \qquad - 4161e_2^2\mu^{1/2} - 2.978 \times 10^4e_2^2\mu + 6429e_2^2\mu^2 \\
& \qquad + 1.82 \times 10^3e_2^3\mu^{1/2} + 7.449 \times 10^4e_2^3\mu \\
& \qquad - 8681e_2^3\mu^2].
\end{aligned}
\tag{25}
\end{equation}
$$

The fit deduced uncertainties of the coefficients in the above equations are tabulated in Appendix A. Note that the above formulae do not represent a series expansion in the parameters $\alpha$, $\mu$, and $e_2$. They represent the best multivariate polynomial fit to numerical simulations. As such, the coefficients do not need to fulfill convergence criteria. In fact, including additional terms in the fit functions would drastically change the numerical coefficients presented. Given the above expressions for $\delta_g$ and $\delta_e$, the corrections to the secular frequency and the forced eccentricity read

$$g_C = g_H [1 - \delta_g], \quad (26)$$

$$\epsilon_C = \epsilon_H [1 - \delta_e]. \quad (27)$$

Applying this correction is straight forward. The secular orbit evolution can be obtained simply by replacing $g_H$ and $\epsilon_H$ by $g_C$ and $\epsilon_C$, respectively, in Equation (8). The resulting secular quantities are valid for small satellite eccentricities and for the range of parameters $0.1 \le e_2 \le 0.6$ and $0.1 \le m_2/m_0 \le 10$, for all stable and non-resonant $\alpha$.

### 4.1. Quality of the Correction

The corrections presented in Equations (24) and (25) are finite polynomial approximations and therefore are still bound to deviate from the reference solution. To assess the quality of the correction, we compare results obtained from the analytical models by Heppenheimer and Marchal with those achieved using Heppenheimer plus our correction and gauge them on the reference solution.<fnref n="10" /> In order to quantify the deviation of the respective approach from the reference solution, we define the standard relative error as

$$\Delta x_y = \frac{|x_y - x_{\text{ref}}|}{x_{\text{ref}}}, \quad (28)$$

where $x = (g, \epsilon)$ are the secular parameters, the index ref stands for the reference solution obtained from the N-body simulations, and $y = (\text{Hepp}, \text{March}, \text{Corr})$ represents the model choice. Figures 4 and 5 show the relative errors in the secular quantities as a function of the semimajor axis ratio $\alpha$.

Both Figures 4 and 5 underline that the results obtained from the corrected model (Corr) are excellent. The errors of the corrected secular frequencies as well as forced eccentricities are much smaller than those resulting from any other method. For the majority of the parameter-space ($\alpha, e_2, \mu$), the relative error is smaller than 5%. The points exhibiting the largest errors are all located close to the stability limits in very localized regions of the ($\alpha, e_2, \mu$) space. The non-optimum behavior there is caused by the presence of significant mean-motion resonances that the numerical method itself is not able to resolve (see Figure 3). The left panels in Figures 4 and 5 show that the error for Hepp increases with larger values of $\alpha, e_2$, and $\mu$, which is

<table><caption>**Table 1**
Initial Parameters of the Examples and Secular Parameters Calculated With Hepp and Corr Models for The Forced Eccentricity and Secular Frequency</caption><thead><tr><th></th><th>HD 41004 Bb<sup>a</sup> (a)</th><th>γ Cephei Ab<sup>b</sup> (b)</th><th>(c)</th><th>(d)</th></tr></thead><tbody><tr><td>m<sub>0</sub> (M<sub>☉</sub>)</td><td>0.42</td><td>1.4</td><td>1</td><td>1</td></tr><tr><td>m<sub>1</sub> (M<sub>☉</sub>)</td><td>1.743 × 10<sup>-2</sup></td><td>1.765 × 10<sup>-3</sup></td><td>10<sup>-4</sup></td><td>10<sup>-4</sup></td></tr><tr><td>m<sub>2</sub> (M<sub>☉</sub>)</td><td>0.7</td><td>0.41</td><td>1</td><td>10</td></tr><tr><td>a<sub>1</sub> (au)</td><td>0.0177</td><td>2.05</td><td>0.17</td><td>0.1</td></tr><tr><td>a<sub>2</sub> (au)</td><td>20.0</td><td>20.2</td><td>1</td><td>1</td></tr><tr><td>e<sub>1</sub></td><td>0.081</td><td>0.05</td><td>0.01</td><td>0.01</td></tr><tr><td>e<sub>2</sub></td><td>0.4</td><td>0.41</td><td>0.2</td><td>0.1</td></tr><tr><td>ε<sub>H</sub></td><td>5.27 × 10<sup>4</sup></td><td>0.063</td><td>0.044</td><td>0.0123</td></tr><tr><td>ε<sub>C</sub></td><td>5.27 × 10<sup>4</sup></td><td>0.057</td><td>0.030</td><td>0.0103</td></tr><tr><td>g<sub>H</sub> (rad yr<sup>-1</sup>)</td><td>1.95 × 10<sup>-6</sup></td><td>7.66 × 10<sup>-4</sup></td><td>0.351</td><td>1.46</td></tr><tr><td>g<sub>C</sub> (rad yr<sup>-1</sup>)</td><td>1.95 × 10<sup>-6</sup></td><td>9.01 × 10<sup>-4</sup></td><td>0.709</td><td>3.98</td></tr></tbody></table>

**Notes.** We assumed for simplicity $M_1 = M_2 = \bar{\omega}_1 = \bar{\omega}_2 = 0$ and we considered masses to be maximal ($m_1 \sin i = m_1$).

<sup>a</sup> Santos et al. (2002), Zucker et al. (2004), Roell et al. (2012).

<sup>b</sup> Neuhäuser et al. (2007), Endl et al. (2011), Reffert & Quirrenbach (2011).

to be expected. The panels in the middle column of both Figures 4 and 5 show that Marchal’s model leads to much smaller errors compared to Heppenheimer’s, especially for small values of $e_2, \mu$. Regarding the forced eccentricity, March presents a peculiar non-monotonous behavior: the error begins increasing, then decreasing, and finally, for some cases, crosses the origin while increasing again. While this suggests that March is applicable to a wider range of parameters, the same is not true for the secular frequency (Figure 4).

## 5. Applicability and Discussion of the Empirically Corrected Secular Solution

In order to demonstrate possible applications and limits of applicability of the corrected secular solution presented in this work, we have selected four different S-type three-body configurations. Two systems correspond to actually detected exoplanets orbiting a single star with a giant planet perturber and one component of a binary star. The other two systems are fictitious, selected to showcase the limits of the various secular theories. Their physical and orbital parameters are presented in Table 1. We also include in Table 1 the values obtained by Hepp and Corr for the forced eccentricity and the secular frequency for each of the systems as well as the stability limits. Figure 6 compares the orbit evolution of the satellite in the systems given in Table 1 as predicted by the four secular models (Hepp, March, A-I, Corr) to the results of a direct N-body integration including non-secular and short-period variations.

The HD 41004 Bb system, example (a), is characterized by a very hierarchical structure ($\alpha \ll 1$) with little short-period activity. All models were able to properly describe its secular dynamics. In contrast, the planet orbiting $\gamma$ Cephei A, presented as example (b), experiences much stronger perturbations, since the system hosts a secondary star relatively close to the primary. This causes Hepp to no longer provide accurate estimates for the secular frequency and forced eccentricity.

If one considers even stronger perturbations, such as in example (c), for instance, we see that both Hepp and March start to diverge from the N-body integration. Note that short-period terms that are not considered in secular dynamics start to play an important role even far from mean-motion resonances (Georgakarakos 2003, 2005). Finally, in example (d), only Corr is

<footnote n="10">¹⁰ We have omitted the A-I model in this comparison for computational reasons, as the evaluation takes somewhat longer than for the other methods.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Examples of the orbits obtained from the integration of the different models for the systems presented in Table 1 compared to the integration of the exact equations of motion.</figcaption></figure>

able to capture the secular dynamics of the system, even when short-period terms are dominating the direct *N*-body integration result. Due to the high-amplitude short-period oscillations present in both the eccentricity $e_1$ and the secular angle $\Delta\varpi$, the secular variable $k = e_1 \cos(\Delta\varpi)$ reaches negative values. In both examples (c) and (d), the free eccentricity $e_p$ is smaller than the forced one $\epsilon$, which leads to a secular oscillation of $\Delta\varpi$ around 0 (see Figure 2). However, in the complete problem, the secular angle alternates between circulation and oscillation around $\Delta\varpi = 0$. This particular characteristic may give the impression that none of the models are capable of reproducing the the *N*-body integration results for $e_1$. Nevertheless, we note that the secular variables (k, h) can still be accurately described by Corr. High-amplitude short-period oscillations are also an indication that higher-order theories should be applied when constructing analytical models to properly describe secular dynamics (Andrade-Ines et al. 2016).

## 6. Conclusions

Simple and accurate models for the long-term orbit evolution of small satellites under the influence of an external perturber are highly desirable, as they find multiple applications in various fields of astrophysics. Although several analytic models exist that describe moderately perturbed systems sufficiently well, they tend to become unwieldy and imprecise in the case of strong perturbations. In this work, we have mapped out an escape route from this dilemma by providing a simple means to calculate the secular dynamics of heavily perturbed three-body systems based on empirical corrections of the simple model by Heppenheimer

(1978). Those corrections were obtained through a fitting procedure bringing secular frequency (g) and forced eccentricity (e) predictions closer to more precise results obtained from direct *N*-body integrations. The resulting Equations (26) and (27) are extremely simple and straight forward to apply. We have shown that our empirical correction can describe all stable secular orbits of planets in hierarchical three-body configurations of “S-type” near perfectly for a wide range of initial conditions, i.e., all semimajor axis ratios $\alpha < 0.4$, perturber eccentricities $e_1 < 0.6$, and perturber mass ratios $0.1 \le \mu \le 10$. Applying our predictions to several exoplanetary systems, we confirm an excellent agreement between our simple secular model and direct *N*-body calculation results.

This research has in part been funded by CNPq Project 204873/2014-2 and by the Jet Propulsion Laboratory through the California Institute of Technology postdoctoral fellowship program, under a contract with the National Aeronautics and Space Administration. The authors would, furthermore, like to acknowledge the support of the IMCCE, Observatoire de Paris, France. The authors would like to thank the anonymous referee who helped to improve this manuscript.

## Appendix A

### Coefficients and Errors of the Correction

In Table 2, the coefficients used in the empirical correction factors for the planet's secular frequency (g, Equation (22)) and its forced eccentricity (e, Equation (23)) are presented. The fit-derived uncertainties for the parameters are also shown.
<table><caption>**Table 2**
Numerical Coefficients of the Empirical Corrections for the Secular Frequency $A_i$ (Equation (22)) and for the Forced Eccentricity $A'_i$ (Equation (23)) and Their Errors $\sigma A_i$ and $\sigma A'_i$, Respectively</caption><thead><tr><th>i</th><th>p<sub>i</sub></th><th>q<sub>i</sub></th><th>l<sub>i</sub></th><th>A<sub>i</sub></th><th>σA<sub>i</sub></th><th>p'<sub>i</sub></th><th>q'<sub>i</sub></th><th>l'<sub>i</sub></th><th>A'<sub>i</sub></th><th>σA'<sub>i</sub></th></tr></thead><tbody><tr><td>1</td><td>1.5</td><td>0</td><td>0.5</td><td>-4.6274</td><td>0.0058</td><td>1.5</td><td>1</td><td>0.5</td><td>29.494</td><td>0.056</td></tr><tr><td>2</td><td>1.5</td><td>0</td><td>1</td><td>-4.0190</td><td>0.0052</td><td>1.5</td><td>1</td><td>1</td><td>9.220</td><td>0.032</td></tr><tr><td>3</td><td>1.5</td><td>0</td><td>2</td><td>0.25041</td><td>0.00039</td><td>1.5</td><td>2</td><td>0.5</td><td>-99.85</td><td>0.30</td></tr><tr><td>4</td><td>1.5</td><td>2</td><td>0.5</td><td>-3.41</td><td>0.11</td><td>1.5</td><td>2</td><td>1</td><td>-31.50</td><td>0.17</td></tr><tr><td>5</td><td>1.5</td><td>2</td><td>1</td><td>11.09</td><td>0.10</td><td>1.5</td><td>3</td><td>0.5</td><td>124.60</td><td>0.36</td></tr><tr><td>6</td><td>1.5</td><td>2</td><td>2</td><td>-0.9823</td><td>0.0077</td><td>1.5</td><td>3</td><td>1</td><td>35.69</td><td>0.20</td></tr><tr><td>7</td><td>1.5</td><td>4</td><td>0.5</td><td>-20.13</td><td>0.31</td><td>4.5</td><td>1</td><td>0.5</td><td>1073.0</td><td>7.5</td></tr><tr><td>8</td><td>1.5</td><td>4</td><td>1</td><td>-85.49</td><td>0.29</td><td>4.5</td><td>1</td><td>1</td><td>4280</td><td>14</td></tr><tr><td>9</td><td>1.5</td><td>4</td><td>2</td><td>4.996</td><td>0.023</td><td>4.5</td><td>1</td><td>2</td><td>-1609.8</td><td>3.5</td></tr><tr><td>10</td><td>4.5</td><td>0</td><td>0.5</td><td>123.67</td><td>0.48</td><td>4.5</td><td>2</td><td>0.5</td><td>-4161</td><td>67</td></tr><tr><td>11</td><td>4.5</td><td>0</td><td>1</td><td>-799.20</td><td>0.88</td><td>4.5</td><td>2</td><td>1</td><td>-2.978 × 10<sup>4</sup></td><td>1.2 × 10<sup>2</sup></td></tr><tr><td>12</td><td>4.5</td><td>0</td><td>2</td><td>-201.49</td><td>0.22</td><td>4.5</td><td>2</td><td>2</td><td>6429</td><td>28</td></tr><tr><td>13</td><td>4.5</td><td>2</td><td>0.5</td><td>180</td><td>23</td><td>4.5</td><td>3</td><td>0.5</td><td>1.82 × 10<sup>3</sup></td><td>1.3 × 10<sup>2</sup></td></tr><tr><td>14</td><td>4.5</td><td>2</td><td>1</td><td>-5555</td><td>40</td><td>4.5</td><td>3</td><td>1</td><td>7.449 × 10<sup>4</sup></td><td>2.3 × 10<sup>2</sup></td></tr><tr><td>15</td><td>4.5</td><td>2</td><td>2</td><td>-617.7</td><td>8.8</td><td>4.5</td><td>3</td><td>2</td><td>-8681</td><td>54</td></tr><tr><td>16</td><td>4.5</td><td>4</td><td>0.5</td><td>2.671 × 10<sup>4</sup></td><td>15</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>17</td><td>4.5</td><td>4</td><td>1</td><td>-1.0229 × 10<sup>5</sup></td><td>2.6 × 10<sup>2</sup></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>18</td><td>4.5</td><td>4</td><td>2</td><td>-23076</td><td>54</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

## Appendix B
Notation

Table 3 contains a list of the variables used in this work together with a brief description.

<table><caption>**Table 3**
Description of the Notation in This Work</caption><thead><tr><th>Variable</th><th>Description</th></tr></thead><tbody><tr><td>m<sub>i</sub></td><td>Mass (central body: i = 0, planet: i = 1, secondary body: i = 2)</td></tr><tr><td>μ</td><td>Mass ratio m<sub>2</sub>/m<sub>0</sub></td></tr><tr><td>a<sub>i</sub></td><td>Semimajor axis (planet: i = 1, secondary body: i = 2)</td></tr><tr><td>α</td><td>Semimajor axis ratio a<sub>1</sub>/a<sub>2</sub></td></tr><tr><td>e<sub>i</sub></td><td>Eccentricity (planet: i = 1, secondary body: i = 2)</td></tr><tr><td>ω<sub>i</sub></td><td>Longitude of the pericenter (planet: i = 1, secondary body: i = 2)</td></tr><tr><td>g</td><td>Gravitational constant</td></tr><tr><td>P<sub>i</sub></td><td>Legendre polynomials</td></tr><tr><td>n<sub>i</sub></td><td>Mean motion (planet: i = 1, secondary body: i = 2)</td></tr><tr><td>k, h</td><td>Coordinates of the planet's eccentricity vector</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>k<sub>0</sub>, h<sub>0</sub></td><td>Initial values of the coordinates of the eccentricity vector</td></tr><tr><td>g<sub>y</sub></td><td>Secular frequency (Hepp: y = H, March: y = M, A-I: y = A, Corr: y = C)</td></tr><tr><td>ε<sub>y</sub></td><td>Forced eccentricity (Hepp: y = H, March: y = M, A-I: y = A, Corr: y = C)</td></tr><tr><td>e<sub>p</sub></td><td>Proper eccentricity</td></tr><tr><td>φ</td><td>Initial phase</td></tr><tr><td>δ<sub>x</sub></td><td>Correction term (secular frequency: x = g, forced eccentricity: x = ε)</td></tr><tr><td>A<sub>i</sub>, A'<sub>i</sub></td><td>Numerical correction coefficients</td></tr><tr><td>p<sub>i</sub>, q<sub>i</sub>, l<sub>i</sub>, p'<sub>i</sub>, q'<sub>i</sub>, l'<sub>i</sub></td><td>Exponents of α, ε<sub>2</sub> and μ in the correction</td></tr><tr><td>Δx<sub>y</sub></td><td>Error (secular frequency: x = g, forced eccentricity: x = ε) for each model (Hepp: y = H, March: y = M, A-I: y = A, Corr: y = C)</td></tr></tbody></table>

## References

Andrade-Ines, E., Beaugé, C., Michtchenko, T., & Robutel, P. 2016, *CeMDA*, 124, 405

Brouwer, D., & Clemence, G. M. 1961, Methods of Celestial Mechanics (New York: Academic Press)

Couetdic, J., Laskar, J., Correia, A. C. M., Mayor, M., & Udry, S. 2010, *A&A*, 519, A10

Dvorak, R. 1986, *A&A*, 167, 379

Eggl, S., & Dvorak, R. 2010, in Lecture Notes in Physics Vol. 790, Dynamics of Small Solar System Bodies and Exoplanets, ed. J. Souchay & R. Dvorak (Berlin: Springer), 431

Endl, M., Cochran, W. D., Hatzes, A. P., & Wittenmyer, R. A. 2011, in AIP Conf. Ser. 1331, Planetary Systems Beyond the Main Sequence, ed. S. Schuh, H. Drechsel, & U. Heber (Melville, NY: AIP), 88

Ferraz-Mello, S. 2007, in *Canonical Perturbation Theories: Degenerate Systems and Resonance*, Vol. 345, Degenerate Systems and Resonance (New York: Springer), 345

Georgakarakos, N. 2003, *MNRAS*, 345, 340

Georgakarakos, N. 2005, *MNRAS*, 362, 748

Giuppone, C. A., Leiva, A. M., Correa-Otto, J., & Beaugé, C. 2011, *A&A*, 530, A103

Heppenheimer, T. A. 1978, *A&A*, 65, 421

Hori, G. 1966, *PASJ*, 18, 287

Kaula, W. M. 1962, *AJ*, 67, 300

Laskar, J. 1990, *Icar*, 88, 266

Laskar, J., & Boué, G. 2010, *A&A*, 522, A60

Leiva, A. M., Correa-Otto, J. A., & Beaugé, C. 2013, *MNRAS*, 436, 3772

Marchal, C. 1990, The Three-body Problem (Amsterdam: Elsevier)

Michtchenko, T. A., Lazzaro, D., Ferraz-Mello, S., & Roig, F. 2002, *Icar*, 158, 343

Michtchenko, T. A., & Malhotra, R. 2004, *Icar*, 168, 237

Neuhäuser, R., Mugrauer, M., Fukagawa, M., Torres, G., & Schmidt, T. 2007, *A&A*, 462, 777

Noyelles, B., Lemaîre, A., & Vienne, A. 2008, *A&A*, 478, 959

Reffert, S., & Quirrenbach, A. 2011, *A&A*, 527, A140

Roell, T., Neuhäuser, R., Seifahrt, A., & Mugrauer, M. 2012, *A&A*, 542, A92

Santos, N. C., Mayor, M., Naef, D., et al. 2002, *A&A*, 392, 215

Sundman, K. F. 1907, Recherches Sur le Problème des Trois Corps (6th ed.; Ex Officina Typographica Societatis Litterariae Fennicae)

Thebault, P. 2011, *CeMDA*, 111, 29

Zucker, S., Mazeh, T., Santos, N. C., Udry, S., & Mayor, M. 2004, *A&A*, 426, 695