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

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.³ Assuming the planet started out on a circular orbit in the same orbital plane as the massive bodies, its eccentricity (e₁) and longitude of pericenter (${\varpi }_{1}$) evolve as

$$\begin{eqnarray}&&{e}_{1}(t)=\epsilon {[4-4\cos ({gt})]}^{1/4},\qquad {\varpi }_{1}(t)=\displaystyle \frac{g}{2}\ t,\end{eqnarray} \tag{ 1 }$$

where and g are the so-called forced eccentricity and secular frequency, respectively. Both 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.⁴ 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.

Let us consider a system consisting of a main body with mass m₀, a planet with mass m₁, and a secondary body with mass m₂ in a reference frame centered on m₀. Let ${{\boldsymbol{r}}}_{i}$ be the position vector of the ith body with respect to m₀ (see Figure 1). We assume, furthermore, that all bodies move in the same plane. If ${m}_{1}\ll {m}_{0}\sim {m}_{2}$ and $| {{\boldsymbol{r}}}_{1}| \lt | {{\boldsymbol{r}}}_{2}|$ for all time and neglecting gravitational effects due to m₁, we arrive at the coplanar restricted three-body problem of S-type (Dvorak 1986). The orbit of m₂ around m₀ 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.

Zoom In Zoom Out Reset image size

2.1. Heppenheimer (1978)

Assuming there are no significant mean-motion resonances or high-amplitude short-period oscillations, the dynamical evolution of the planets's orbit (m₁) will be dominated by the secular interaction with the perturber (m₂). 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₃ (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:

$$\begin{eqnarray}{{ \mathcal R }}_{H} & = & \displaystyle \frac{{ \mathcal G }{m}_{2}}{{(1-{e}_{2}^{2})}^{3/2}}\displaystyle \frac{{a}_{1}^{2}}{{a}_{2}^{3}}\left[\displaystyle \frac{1}{4}+\displaystyle \frac{3}{8}{e}_{1}^{2}\right.\\ & & \left.-\displaystyle \frac{15}{16}\displaystyle \frac{{a}_{1}}{{a}_{2}}\displaystyle \frac{{e}_{1}{e}_{2}}{(1-{e}_{2}^{2})}\cos ({\rm{\Delta }}\varpi )\right],\end{eqnarray} \tag{ 2 }$$

where ${ \mathcal G }$ is the gravitational constant, ${\rm{\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 ith 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,

$$\begin{eqnarray}&&{{ \mathcal R }}_{H}={n}_{1}{a}_{1}^{2}{g}_{H}\left[\displaystyle \frac{{e}_{1}^{2}}{2}-{\epsilon }_{H}{e}_{1}\cos {\rm{\Delta }}\varpi \right],\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{g}_{H}=\displaystyle \frac{3}{4}\mu {\alpha }^{3}\displaystyle \frac{{n}_{1}}{{(1-{e}_{2}^{2})}^{3/2}}\end{eqnarray} \tag{ 4 }$$

is the secular frequency and

$$\begin{eqnarray}&&{\epsilon }_{H}=\displaystyle \frac{5}{4}\alpha \displaystyle \frac{{e}_{2}}{(1-{e}_{2}^{2})}\end{eqnarray} \tag{ 5 }$$

is the forced eccentricity, with ${n}_{1}=\sqrt{{ \mathcal G }{m}_{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{eqnarray}k & = & {e}_{1}\cos {\rm{\Delta }}\varpi ,\\ h & = & {e}_{1}\sin {\rm{\Delta }}\varpi ,\end{eqnarray} \tag{ 6 }$$

the Lagrange–Laplace planetary equations up to order ${ \mathcal O }({e}_{1}^{2})$, which determine the satellite's orbit evolution, read (Brouwer & Clemence 1961):

$$\begin{eqnarray}\displaystyle \frac{{dh}}{{dt}} & = & \displaystyle \frac{1}{{n}_{1}{a}_{1}^{2}}\displaystyle \frac{\partial { \mathcal R }}{\partial k}={g}_{H}(k-{\epsilon }_{H}),\\ \displaystyle \frac{{dk}}{{dt}} & = & -\displaystyle \frac{1}{{n}_{1}{a}_{1}^{2}}\displaystyle \frac{\partial { \mathcal R }}{\partial h}=-{g}_{H}h.\end{eqnarray} \tag{ 7 }$$

The system of Equation (7) admits the analytical solution

$$\begin{eqnarray}k & = & {e}_{p}\cos ({g}_{H}t+\phi )+{\epsilon }_{H},\\ h & = & {e}_{p}\sin ({g}_{H}t+\phi ),\end{eqnarray} \tag{ 8 }$$

where e_p and ϕ are constants of integration, obtained from the initial conditions ${e}_{1}(0)$ and ${\varpi }_{1}(0)$ as

$$\begin{eqnarray}&&{e}_{p}=\sqrt{{({k}_{0}-{\epsilon }_{H})}^{2}+{h}_{0}^{2}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\cos \phi =\displaystyle \frac{{k}_{0}-{\epsilon }_{H}}{{e}_{p}},\ \sin \phi =\displaystyle \frac{{h}_{0}}{{e}_{p}}.\end{eqnarray} \tag{ 10 }$$

Here, ${k}_{0}={e}_{1}(0)\ \cos {\rm{\Delta }}\varpi (0)$ and ${h}_{0}={e}_{1}(0)\sin {\rm{\Delta }}\varpi (0)$, with ${\rm{\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₁ and ${\rm{\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₁ and ${\rm{\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₁, $\alpha ={a}_{1}/{a}_{2}$, $\mu ={m}_{2}/{m}_{0}$, and e₂ 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 same secular frequency, no matter the initial values of ${\varpi }_{1}$ and e₁.⁵ 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}\gt 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.

Zoom In Zoom Out Reset image size

2.2. Marchal (1990)

Using the Von Zeipel averaging method on the three-body disturbing function, Marchal (1990) obtained a second-order⁶ 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

$$\begin{eqnarray}&&{{ \mathcal R }}_{M}={Q}_{1}+{Q}_{2}+{Q}_{3},\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{Q}_{1}=\displaystyle \frac{1}{4}{ \mathcal G }\displaystyle \frac{{m}_{0}{m}_{2}}{{m}_{0}+{m}_{1}}\displaystyle \frac{{a}_{1}^{2}}{{a}_{2}^{3}}\left(1+\displaystyle \frac{3}{2}{e}_{1}^{2}\right)\displaystyle \frac{1}{{(1-{e}_{2}^{2})}^{3/2}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}{Q}_{2} & = & -\displaystyle \frac{15}{16}{ \mathcal G }\displaystyle \frac{{m}_{0}{m}_{2}({m}_{0}-{m}_{1})}{{({m}_{0}+{m}_{1})}^{2}}\displaystyle \frac{{a}_{1}^{3}}{{a}_{2}^{4}}\\ & & \times \ {e}_{1}\left(1+\displaystyle \frac{3}{4}{e}_{1}^{2}\right)\displaystyle \frac{{e}_{2}}{{(1-{e}_{2}^{2})}^{5/2}}\cos {\rm{\Delta }}\varpi \end{eqnarray} \tag{ 13 }$$

originate in the Legendre polynomial expansion up to P₃ and

$$\begin{eqnarray}{Q}_{3} & = & \displaystyle \frac{15}{64}{ \mathcal G }\displaystyle \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}\displaystyle \frac{[5(3+2{e}_{2}^{2})+3{e}_{2}^{2}\cos (2{\rm{\Delta }}\varpi )]}{{(1-{e}_{2}^{2})}^{3}}\end{eqnarray} \tag{ 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 ${ \mathcal O }({e}_{1}^{3})$, Equation (11) becomes

$$\begin{eqnarray}&&{{ \mathcal R }}_{M}={n}_{1}{a}_{1}^{2}{g}_{M}\left[\displaystyle \frac{{e}_{1}^{2}}{2}-{\epsilon }_{M}{e}_{1}\cos {\rm{\Delta }}\varpi \right],\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{g}_{M}={g}_{H}\left[1+\displaystyle \frac{25}{8}\displaystyle \frac{\mu }{\sqrt{1+\mu }}{\alpha }^{3/2}\displaystyle \frac{3+2{e}_{2}^{2}}{{(1-{e}_{2}^{2})}^{3/2}}\right],\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\epsilon }_{M}={\epsilon }_{H}{\left[1+\displaystyle \frac{25}{8}\displaystyle \frac{\mu }{\sqrt{1+\mu }}{\alpha }^{3/2}\displaystyle \frac{3+2{e}_{2}^{2}}{{(1-{e}_{2}^{2})}^{3/2}}\right]}^{-1}.\end{eqnarray} \tag{ 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⁷ 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

$$\begin{eqnarray}&&{\displaystyle \frac{\partial {{ \mathcal H }}_{A}}{\partial {e}_{1}}\right|}_{{e}_{1}={\epsilon }_{A},{\rm{\Delta }}\varpi =0}=0.\end{eqnarray} \tag{ 18 }$$

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.

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 (α, μ, and e₂), it can be calculated in advance. Returning to Figure 2 we see that the secular evolution of the planet's orbit in the 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 ϕ, 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 ${\rm{\Delta }}\varpi (t)=0$. In other words, the planet's orbit (a₁, e₁, ${\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⁸ planar problem, however, we have three degrees of freedom, associated with the following three frequencies: the mean motions of the planet (n₁) and the perturber (n₂) and the secular frequency g. In this case, the stationary secular solution corresponds to a quasi-periodic orbit in osculating elements with frequencies n₁ and n₂ 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{\rm{\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

$$\begin{eqnarray}&&{\epsilon }_{E}=\langle {e}_{1}(t)\cos {\rm{\Delta }}\varpi (t){\rangle }_{T},\end{eqnarray} \tag{ 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⁹ 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₁ and n₂. 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 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 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 $\alpha \approx 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.

Zoom In Zoom Out Reset image size

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 ${\epsilon }_{H}$ and the secular frequency g_H. 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 g_H and the forced eccentricity ${\epsilon }_{H}$, respectively, as

$$\begin{eqnarray}&&{\delta }_{g}(\alpha ,{e}_{2},\mu )=1-\displaystyle \frac{{g}_{E}}{{g}_{H}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\delta }_{\epsilon }(\alpha ,{e}_{2},\mu )=1-\displaystyle \frac{{\epsilon }_{E}}{{\epsilon }_{H}},\end{eqnarray} \tag{ 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 α for each pair of parameters $(\mu ,{e}_{2})$. For all simulations, we fixed the mass of the central body at ${m}_{0}=1{M}_{\odot }$ and the semimajor axis of the perturber at ${a}_{2}=1\,\mathrm{au}$. We then applied the method of least squares to minimize the error functions ${\delta }_{g}$ and ${\delta }_{\epsilon }$, 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₂, in the form of

$$\begin{eqnarray}&&{\delta }_{g}=\displaystyle \sum _{i=1}^{{N}_{g}}\,{A}_{i}{\alpha }^{{p}_{i}}{e}_{2}^{{q}_{i}}{\mu }^{{l}_{i}},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\delta }_{\epsilon }=\displaystyle \sum _{i=1}^{{N}_{e}}\,{A}_{i}^{\prime }{\alpha }^{{p}_{i}^{\prime }}{e}_{2}^{{q}_{i}^{\prime }}{\mu }^{{l}_{i}^{\prime }}.\end{eqnarray} \tag{ 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 }_{\epsilon }$ 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{eqnarray}{\delta }_{g} & = & {\alpha }^{3/2}[-4.6274{\mu }^{1/2}-4.0190\mu +0.25041{\mu }^{2}\\ & & -3.41{e}_{2}^{2}{\mu }^{1/2}+11.09{e}_{2}^{2}\mu -0.9823{e}_{2}^{2}{\mu }^{2}\\ & & -20.13{e}_{2}^{4}{\mu }^{1/2}-85.49{e}_{2}^{4}\mu +4.996{e}_{2}^{4}{\mu }^{2}]\\ & & +{\alpha }^{9/2}[123.67{\mu }^{1/2}-799.20\mu -201.49{\mu }^{2}\\ & & +180{e}_{2}^{2}{\mu }^{1/2}-5555{e}_{2}^{2}\mu -617.7{e}_{2}^{2}{\mu }^{2}\\ & & +2.671\times {10}^{4}{e}_{2}^{4}{\mu }^{1/2}-1.0229\times {10}^{5}{e}_{2}^{4}\mu \\ & & -23076{e}_{2}^{4}{\mu }^{2}],\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}{\delta }_{\epsilon } & = & {\alpha }^{3/2}[29.494{e}_{2}{\mu }^{1/2}+9.220{e}_{2}\mu -99.85{e}_{2}^{2}{\mu }^{1/2}\\ & & -31.50{e}_{2}^{2}\mu +124.60{e}_{2}^{3}{\mu }^{1/2}+35.69{e}_{2}^{3}\mu ]\\ & & +{\alpha }^{9/2}[1073.0{e}_{2}{\mu }^{1/2}+4280{e}_{2}\mu -1609.8{e}_{2}{\mu }^{2}\\ & & -4161{e}_{2}^{2}{\mu }^{1/2}-2.978\times {10}^{4}{e}_{2}^{2}\mu +6429{e}_{2}^{2}{\mu }^{2}\\ & & +1.82\times {10}^{3}{e}_{2}^{3}{\mu }^{1/2}+7.449\times {10}^{4}{e}_{2}^{3}\mu \\ & & -8681{e}_{2}^{3}{\mu }^{2}].\end{eqnarray} \tag{ 25 }$$

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 α, μ, and e₂. 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 }_{\epsilon }$, the corrections to the secular frequency and the forced eccentricity read

$$\begin{eqnarray}&&{g}_{C}={g}_{H}[1-{\delta }_{g}],\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\epsilon }_{C}={\epsilon }_{H}[1-{\delta }_{\epsilon }].\end{eqnarray} \tag{ 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\leqslant {e}_{2}\leqslant 0.6$ and $0.1\leqslant {m}_{2}/{m}_{0}\leqslant 10$, for all stable and non-resonant α.

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.¹⁰ In order to quantify the deviation of the respective approach from the reference solution, we define the standard relative error as

$$\begin{eqnarray}&&{\rm{\Delta }}{x}_{y}=\displaystyle \frac{| {x}_{y}-{x}_{\mathrm{ref}}| }{{x}_{\mathrm{ref}}},\end{eqnarray} \tag{ 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=({Hepp},{March},{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 α.

Zoom In Zoom Out Reset image size

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 α, e₂, and μ, which is 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₂, μ. Regarding the forced eccentricity, March presents a peculiar non-monotonous behavior: the error begin 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).

Zoom In Zoom Out Reset image size

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 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.

Zoom In Zoom Out Reset image size

Table 1.  Initial Parameters of the Examples and Secular Parameters Calculated With Hepp and Corr Models for The Forced Eccentricity and Secular Frequency

	HD 41004 Bba (a)	γ Cephei Abb (b)	(c)	(d)
${m}_{0}({M}_{\odot })$	0.42	1.4	1	1
${m}_{1}({M}_{\odot })$	$1.743\times {10}^{-2}$	$1.765\times {10}^{-3}$	10−4	10−4
${m}_{2}({M}_{\odot })$	0.7	0.41	1	10
${a}_{1}({au})$	0.0177	2.05	0.17	0.1
${a}_{2}({au})$	20.0	20.2	1	1
e1	0.081	0.05	0.01	0.01
e2	0.4	0.41	0.2	0.1
${\epsilon }_{H}$	$5.27\times {10}^{4}$	0.063	0.044	0.0123
${\epsilon }_{C}$	$5.27\times {10}^{4}$	0.057	0.030	0.0103
${g}_{H}(\mathrm{rad}\,{\mathrm{yr}}^{-1})$	$1.95\times {10}^{-6}$	$7.66\times {10}^{-4}$	0.351	1.46
${g}_{C}(\mathrm{rad}\,{\mathrm{yr}}^{-1})$	$1.95\times {10}^{-6}$	$9.01\times {10}^{-4}$	0.709	3.98

Notes. We assumed for simplicity ${M}_{1}={M}_{2}={\varpi }_{1}={\varpi }_{2}=0$ and we considered masses to be maximal (${m}_{1}\sin i={m}_{1}$).

^aSantos et al. (2002), Zucker et al. (2004), Roell et al. (2012). ^bNeuhäuser et al. (2007), Endl et al. (2011), Reffert & Quirrenbach (2011).

Download table as:  ASCIITypeset image

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 γ 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 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₁ and the secular angle ${\rm{\Delta }}\varpi$, the secular variable $k={e}_{1}\cos ({\rm{\Delta }}\varpi )$ reaches negative values. In both examples (c) and (d), the free eccentricity e_p is smaller than the forced one , which leads to a secular oscillation of ${\rm{\Delta }}\varpi$ around 0 (see Figure 2). However, in the complete problem, the secular angle alternates between circulation and oscillation around ${\rm{\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₁. 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).

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 () 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 \lt 0.4$, perturber eccentricities ${e}_{1}\lt 0.6$, and perturber mass ratios $0.1\leqslant \mu \leqslant 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.

In Table 2, the coefficients used in the empirical correction factors for the planet's secular frequency (g, Equation (22)) and its forced eccentricity (, Equation (23)) are presented. The fit-derived uncertainties for the parameters are also shown.

Table 2.  Numerical Coefficients of the Empirical Corrections for the Secular Frequency A_i (Equation (22)) and for the Forced Eccentricity ${A}_{i}^{\prime }$ (Equation (23)) and Their Errors $\sigma {A}_{i}$ and $\sigma {A}_{i}^{\prime }$, Respectively

i	pi	qi	li	Ai	$\sigma {A}_{i}$	${p}_{i}^{\prime }$	${q}_{i}^{\prime }$	${l}_{i}^{\prime }$	${A}_{i}^{\prime }$	$\sigma {A}_{i}^{\prime }$
1	1.5	0	0.5	−4.6274	0.0058	1.5	1	0.5	29.494	0.056
2	1.5	0	1	−4.0190	0.0052	1.5	1	1	9.220	0.032
3	1.5	0	2	0.25041	0.00039	1.5	2	0.5	−99.85	0.30
4	1.5	2	0.5	−3.41	0.11	1.5	2	1	−31.50	0.17
5	1.5	2	1	11.09	0.10	1.5	3	0.5	124.60	0.36
6	1.5	2	2	−0.9823	0.0077	1.5	3	1	35.69	0.20
7	1.5	4	0.5	−20.13	0.31	4.5	1	0.5	1073.0	7.5
8	1.5	4	1	−85.49	0.29	4.5	1	1	4280	14
9	1.5	4	2	4.996	0.023	4.5	1	2	−1609.8	3.5
10	4.5	0	0.5	123.67	0.48	4.5	2	0.5	−4161	67
11	4.5	0	1	−799.20	0.88	4.5	2	1	$-2.978\times {10}^{4}$	$1.2\times {10}^{2}$
12	4.5	0	2	−201.49	0.22	4.5	2	2	6429	28
13	4.5	2	0.5	180	23	4.5	3	0.5	$1.82\times {10}^{3}$	$1.3\times {10}^{2}$
14	4.5	2	1	−5555	40	4.5	3	1	$7.449\times {10}^{4}$	$2.3\times {10}^{2}$
15	4.5	2	2	−617.7	8.8	4.5	3	2	−8681	54
16	4.5	4	0.5	$2.671\times {10}^{4}$	15	⋯	⋯	⋯	⋯	⋯
17	4.5	4	1	$-1.0229\times {10}^{5}$	$2.6\times {10}^{2}$	⋯	⋯	⋯	⋯	⋯
18	4.5	4	2	−23076	54	⋯	⋯	⋯	⋯	⋯

Download table as:  ASCIITypeset image

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

Table 3.  Description of the Notation in This Work

Variable	Description
mi	Mass (central body: i = 0, planet: i = 1, secondary body: i = 2)
μ	Mass ratio ${m}_{2}/{m}_{0}$
ai	Semimajor axis (planet: i = 1, secondary body: i = 2)
α	Semimajor axis ratio ${a}_{1}/{a}_{2}$
ei	Eccentricity (planet: i = 1, secondary body: i = 2)
${\varpi }_{i}$	Longitude of the pericenter (planet: i = 1, secondary body: i = 2)
${ \mathcal G }$	Gravitational constant
Pi	Legendre polynomials
ni	Mean motion (planet: i = 1, secondary body: i = 2)
$k,\ h$	Coordinates of the planet's eccentricity vector
t	Time
${k}_{0},\ {h}_{0}$	Initial values of the coordinates of the eccentricity vector
gy	Secular frequency (Hepp: y = H, March: y = M, A-I: y = A, Corr: y = C)
${\epsilon }_{y}$	Forced eccentricity (Hepp: y = H, March: y = M, A-I: y = A, Corr: y = C)
ep	Proper eccentricity
ϕ	Initial phase
${\delta }_{x}$	Correction term (secular frequency: x = g, forced eccentricity: $x=\epsilon$)
${A}_{i},\ {A}_{i}^{\prime }$	Numerical correction coefficients
${p}_{i},{q}_{i},{l}_{i},\ {p}_{i}^{\prime },{q}_{i}^{\prime },{l}_{i}^{\prime }$	Exponents of α, e2 and μ in the correction
${\rm{\Delta }}{x}_{y}$	Error (secular frequency: x = g, forced eccentricity: $x=\epsilon$)
	for each model (Hepp: y = H, March: y = M, A-I: y = A, Corr: y = C)

Download table as:  ASCIITypeset image