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

TL;DR

This paper introduces a semi-analytical correction to the classical secular orbit model, improving its accuracy for systems with strong external perturbers by using polynomial correction factors validated against numerical simulations.

Contribution

A simple, accurate correction model for secular orbit evolution that enhances Heppenheimer's classical solution using polynomial functions based on system parameters.

Findings

Corrected secular equations match numerical simulations more closely.

Model applicable across a wide range of parameters.

Provides limits of applicability for the correction.

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 (1978). A comparison between analytical predictions and numerical simulations allows us to determine corrective factors for the secular frequency and forced eccentricity in the co-planar 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 Heppenheimer's 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.

Tables

Table 1. Table 1: Initial parameters of the examples and secular parameters calculated with Hepp and Corr models for the forced eccentricity and secular frequency.

	HD 41004Bba (a)	$\gamma$ Cephei Abb (b)	(c)	(d)
$m_0\left(M_{\odot }\right)$	0.42	1.4	1	1
$m_1\left(M_{\odot }\right)$	$1.743\times {10}^{-2}$	$1.765\times {10}^{-3}$	${10}^{-4}$	${10}^{-4}$
$m_2\left(M_{\odot }\right)$	0.7	0.41	1	10
$a_1\left(au\right)$	0.0177	2.05	0.17	0.1
$a_2\left(au\right)$	20.0	20.2	1	1
$e_1$	0.081	0.05	0.01	0.01
$e_2$	0.4	0.41	0.2	0.1
${ϵ}_H$	$5.27\times {10}^4$	0.063	0.044	0.0123
${ϵ}_C$	$5.27\times {10}^4$	0.057	0.030	0.0103
$g_H\left(rad/yr\right)$	$1.95\times {10}^{-6}$	$7.66\times {10}^{-4}$	0.351	1.46
$g_C\left(rad/yr\right)$	$1.95\times {10}^{-6}$	$9.01\times {10}^{-4}$	0.709	3.98

Table 2. Table 2: Numerical coefficients of the empirical corrections for the secular frequency A i subscript 𝐴 𝑖 A_{i} (Equation ( 24 ))and for the forced eccentricity A i ′ superscript subscript 𝐴 𝑖 ′ A_{i}^{\prime} (Equation ( 25 )) and their errors σ ​ A i 𝜎 subscript 𝐴 𝑖 \sigma A_{i} and σ ​ A i ′ 𝜎 superscript subscript 𝐴 𝑖 ′ \sigma A_{i}^{\prime} , respectively.

$i$	$p_i$	$q_i$	$l_i$	$A_i$	$\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$	$-$	$-$	$-$	$-$	$-$

Table 3. Table 3: Description of the notation in this work.

Variable	Description
$m_i$	Mass (central body: $i=0$, planet: $i=1$, secondary body: $i=2$)
$\mu$	Mass ratio $m_2/m_0$
$a_i$	Semimajor axis (planet: $i=1$, secondary body: $i=2$)
$\alpha$	Semimajor axis ratio $a_1/a_2$
$e_i$	Eccentricity (planet: $i=1$, secondary body: $i=2$)
${\varpi }_i$	Longitude of the pericenter (planet: $i=1$, secondary body: $i=2$)
$𝒢$	Gravitational constant
$P_i$	Legendre polynomials
$n_i$	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
$g_y$	Secular frequency (Hepp: $y=H$, March: $y=M$, A-I: $y=A$, Corr: $y=C$)
${ϵ}_y$	Forced eccentricity (Hepp: $y=H$, March: $y=M$, A-I: $y=A$, Corr: $y=C$)
$e_p$	Proper eccentricity
$\varphi$	Initial phase
${\delta }_x$	Correction term (secular frequency: $x=g$, forced eccentricity: $x=ϵ$)
$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 $\alpha$, $e_2$ and $\mu$ in the correction
$\Delta x_y$	Error (secular frequency: $x=g$, forced eccentricity: $x=ϵ$)
	for each model (Hepp: $y=H$, March: $y=M$, A-I: $y=A$, Corr: $y=C$)

Equations

$$e_{1}(t)=\epsilon\big{[}4-4\cos{(gt)}\big{]}^{1/4},\qquad\varpi_{1}(t)=\frac{g}{2}\;t,$$

$${\mathcal{R}}_H=\frac{\mathcal{G}{m}_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],$$

$${\mathcal{R}}_H=n_1{a}_1^2{g}_H\left[\frac{e_1^2}{2}-{ϵ}_H{e}_1\cos \Delta \varpi \right],$$

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

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

$$\displaystyle\begin{array}[]{rcl}k&=e_{1}\cos\Delta\varpi,\\ h&=e_{1}\sin\Delta\varpi,\end{array}$$

$$\begin{array}[]{rcl}\vspace{0.2cm}\displaystyle\frac{dh}{dt}=&\displaystyle\frac{1}{n_{1}a_{1}^{2}}\displaystyle\frac{\partial{\cal R}}{\partial k}&=g_{H}\left(k-\epsilon_{H}\right),\\ \displaystyle\frac{dk}{dt}=&-\displaystyle\frac{1}{n_{1}a_{1}^{2}}\displaystyle\frac{\partial{\cal R}}{\partial h}&=-g_{H}h.\end{array}$$

$$\begin{array}[]{rl}k=&e_{p}\cos(g_{H}t+\phi)+\epsilon_{H},\\ h=&e_{p}\sin(g_{H}t+\phi),\end{array}$$

$$e_p=\sqrt{(k_0-{ϵ}_H{)}^2+h_0^2},$$

$$\cos \varphi =\frac{k_0-{ϵ}_H}{e_p},\sin \varphi =\frac{h_0}{e_p}.$$

$${\mathcal{R}}_M=Q_1+Q_2+Q_3$$

$$Q_1=\frac{1}{4}\mathcal{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}},$$

$$\begin{array}[]{rl}\vspace{0.4cm}Q_{2}=&-\displaystyle\frac{15}{16}{\cal G}\displaystyle\frac{m_{0}m_{2}(m_{0}-m_{1})}{(m_{0}+m_{1})^{2}}\displaystyle\frac{a_{1}^{3}}{a_{2}^{4}}\\ \vspace{0.4cm}\hfil&\times\ e_{1}\left(1+\displaystyle\frac{3}{4}e_{1}^{2}\right)\displaystyle\frac{e_{2}}{(1-e_{2}^{2})^{5/2}}\cos\Delta\varpi\end{array}$$

$$\begin{array}[]{rl}\vspace{0.4cm}Q_{3}=&\displaystyle\frac{15}{64}{\cal G}\displaystyle\frac{m_{0}m_{2}^{2}}{(m_{0}+m_{1})^{3/2}(m_{0}+m_{1}+m_{2})^{1/2}}\\ \vspace{0.4cm}\hfil&\times\ e_{1}^{2}(1-e_{1}^{2})^{1/2}\displaystyle\frac{\left[5(3+2e_{2}^{2})+3e_{2}^{2}\cos(2\Delta\varpi)\right]}{(1-e_{2}^{2})^{3}}\end{array}$$

$${\mathcal{R}}_M=n_1{a}_1^2{g}_M\left[\frac{e_1^2}{2}-{ϵ}_M{e}_1\cos \Delta \varpi \right],$$

$$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],$$

$${ϵ}_M={ϵ}_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}.$$

$$\frac{\partial{\cal H}_{A}}{\partial e_{1}}\Big{|}_{e_{1}=\epsilon_{A},\Delta\varpi=0}=0.$$

$${ϵ}_E=⟨e_1(t)\cos \Delta \varpi (t){⟩}_T,$$

$${\delta }_g(\alpha ,e_2,\mu )=1-\frac{g_E}{g_H},$$

$${\delta }_{ϵ}(\alpha ,e_2,\mu )=1-\frac{{ϵ}_E}{{ϵ}_H},$$

$${\delta }_g=\sum _{i=1}^{N_g}{A}_i{\alpha }^{p_i}{e}_2^{q_i}{\mu }^{l_i},$$

$${\delta }_{ϵ}=\sum _{i=1}^{N_e}{A}_i^{\prime }{\alpha }^{p_i^{\prime }}{e}_2^{q_i^{\prime }}{\mu }^{l_i^{\prime }}.$$

$$\begin{array}[]{rl}\vspace{0.2cm}\delta_{g}=\alpha^{3/2}&[-4.6274\mu^{1/2}-4.0190\mu+0.25041\mu^{2}\\ \vspace{0.2cm}\hfil&-3.41e_{2}^{2}\mu^{1/2}+11.09e_{2}^{2}\mu-0.9823e_{2}^{2}\mu^{2}\\ \vspace{0.2cm}\hfil&-20.13e_{2}^{4}\mu^{1/2}-85.49e_{2}^{4}\mu+4.996e_{2}^{4}\mu^{2}]\\ \vspace{0.2cm}+\alpha^{9/2}&[123.67\mu^{1/2}-799.20\mu-201.49\mu^{2}\\ \vspace{0.2cm}\hfil&+180e_{2}^{2}\mu^{1/2}-5555e_{2}^{2}\mu-617.7e_{2}^{2}\mu^{2}\\ \vspace{0.2cm}\hfil&+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{array}$$

$$\begin{array}[]{rl}\vspace{0.2cm}\delta_{\epsilon}=\alpha^{3/2}&[29.494e_{2}\mu^{1/2}+9.220e_{2}\mu-99.85e_{2}^{2}\mu^{1/2}\\ \vspace{0.2cm}\hfil&-31.50e_{2}^{2}\mu+124.60e_{2}^{3}\mu^{1/2}+35.69e_{2}^{3}\mu]\\ \vspace{0.2cm}+\alpha^{9/2}&[1073.0e_{2}\mu^{1/2}+4280e_{2}\mu-1609.8e_{2}\mu^{2}\\ \vspace{0.2cm}\hfil&-4161e_{2}^{2}\mu^{1/2}-2.978\times 10^{4}e_{2}^{2}\mu+6429e_{2}^{2}\mu^{2}\\ \vspace{0.2cm}\hfil&+1.82\times 10^{3}e_{2}^{3}\mu^{1/2}+7.449\times 10^{4}e_{2}^{3}\mu\\ &-8681e_{2}^{3}\mu^{2}].\\ \end{array}$$

$$g_C=g_H[1-{\delta }_g],$$

$${ϵ}_C={ϵ}_H[1-{\delta }_{ϵ}].$$

$$\Delta x_y=\frac{|x_y-x_{ref}|}{x_{ref}},$$

Full text

celestial mechanics – planets and satellites: dynamical evolution and stability – planets and satellites: formation – (stars:) binaries: general

\AuthorCallLimit

=2 \fullcollaborationNameThe Friends of AASTeX Collaboration

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, 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 masses111That is, neglecting the dynamical effects due to the shapes of the bodies.. 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

[TABLE]

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, Equations (1) completely describe the secular evolution of the system222This 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.. The non-linear 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). Yet, for large scale planet formation simulations in particular, the time requirements to individually propagate each particle through numerical integrations are still prohibitive (Leiva et al., 2013; Thebault, 2011). 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 Equations (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.

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 $\vec{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 $|{\vec{r}_{1}}|<|\vec{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, the long-term evolution of the planets’s 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

[TABLE]

where ${\cal 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

[TABLE]

where

[TABLE]

is the secular frequency and

[TABLE]

is the forced eccentricity, with $n_{1}=\sqrt{{\cal 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

[TABLE]

the Lagrange-Laplace planetary equations up to order $\mathcal{O}(e_{1}^{2})$ that determine the satellite’s orbit evolution read (Brouwer & Clemence, 1961):

[TABLE]

The system of Equations (9) admits the analytical solution

[TABLE]

where $e_{p}$ and $\phi$ are constants of integration, obtained from the initial conditions $e_{1}(0)$ and $\varpi_{1}(0)$ as

[TABLE]

[TABLE]

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 Equations (10) 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 Equations (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 same secular frequency, no matter the initial values of $\varpi_{1}$ and $e_{1}$.333This statement holds for models up to order ${\mathcal{O}}(e_{1}^{2})$. As the planet’s 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 Equations (10) reduce to Equations (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-order444in $\mu$ and $\alpha$, see also Georgakarakos (2003, 2005) 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

[TABLE]

where

[TABLE]

[TABLE]

originate in the Legendre polynomial expansion up to $P_{3}$ and

[TABLE]

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 (13) becomes

[TABLE]

where

[TABLE]

[TABLE]

The disturbing function ${\cal R}_{M}$ in Equation (17) has a form similar to ${\cal 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 Equations (10) 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 (18) and (19) relate Marchal’s secular quantities with those in Heppenheimer’s model. Only if $\mu\ll 1$, (i.e. $m_{2}\ll m_{0}$), 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-order555in $\alpha^{2}\cdot\mu$ 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 (${\cal 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

[TABLE]

Equation (20) 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 Equations (10), 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 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 complete666non-secular, i.e. not averaged over the fast angles 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}\textrm{e}^{{\textrm{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

[TABLE]

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 long777at least several secular periods 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 integration of $T=12\pi/g_{H}$, that is 6 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, specially 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 \big{(}k(0),h(0)\big{)} the average of the same variables \big{(}k(t),h(t)\big{)} 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 $\gamma$ 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 $\alpha$, all the models coincide with the reference solution. For large $\alpha$ values, however, none of the models provides 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 $\alpha$. 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 $\gamma$ 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 a 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 (17) 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 Equations (10). 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 (18) and (19) 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 $\epsilon$. Those, however, were only valid for the particular case of the $\gamma$-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

[TABLE]

[TABLE]

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\textrm{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_{2}$, in the form of

[TABLE]

[TABLE]

Naturally, one can expect 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}$ was 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:

[TABLE]

[TABLE]

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 fulfil 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:

[TABLE]

[TABLE]

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 Equations (10). The resulting secular quantities are valid for small satellite eccentricities and for the range of parameters $0.1\leq e_{2}\leq 0.6$ and $0.1\leq m_{2}/m_{0}\leq 10$, for all stable and non-resonant $\alpha$.

4.1 Quality of the Correction

The corrections presented in Equations (26) and (27) 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.888We have omitted the A-I model in this comparison for computational reasons, as the evaluation takes somewhat longer than for the other methods. In order to quantify the deviation of the respective approach from the reference solution we define the standard relative error as

[TABLE]

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)$ represent 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 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 begin increasing, then decrease and finally, for some cases, cross 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 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 seconday 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_{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 is 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, 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 ($\epsilon$) predictions closer to more precise results obtained from direct N-body integrations. The resulting Equations (28) and (29) 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\leq\mu\leq 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 (24)) and its forced eccentricity ($\epsilon$, Equation (25)) are presented. The fit-derived uncertainties for the parameters are also shown.

Appendix B Notation

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