Orbital flips in hierarchical triple systems: relativistic effects and third-body effects to hexadecapole order

TL;DR

This paper extends the analytical understanding of hierarchical triple systems by including hexadecapole and relativistic effects, showing that orbital flips and high eccentricities are robust features influenced by higher-order perturbations.

Contribution

The study derives secular evolution equations for hierarchical triples up to hexadecapole order, incorporating relativistic effects, and demonstrates the persistence of orbital flips under these higher-order perturbations.

Findings

Orbital flips are generally robust against relativistic and hexadecapole effects.

Hexadecapole contributions can induce flips in equal-mass inner binaries.

Relativistic precession does not suppress orbital flips in most cases.

Abstract

We analyze the secular evolution of hierarchical triple systems in the post-Newtonian approximation to general relativity. We expand the Newtonian three-body equations of motion in powers of the ratio $a/A$, where $a$ and $A$ are the semimajor axis of the inner binary's orbit and of the orbit of the third body relative to the center of mass of the inner binary, respectively. The leading order "quadrupole" terms, of order $(a/A)^3$ relative to the $1/a^2$ acceleration within the inner binary, are responsible for the well-known Kozai-Lidov oscillations of orbital inclination and eccentricity. The octupole terms, of order $(a/A)^4$ have been shown to allow the inner orbit to "flip" from prograde relative to the outer orbit to retrograde and back, and to permit excursions to very large eccentricities. We carry the expansion of the equations of motion to hexadecapole order, corresponding to contributions of order $(a/A)^5$. We also include the leading orbital effects of post-Newtonian theory, namely the pericenter precessions of the inner and outer orbits. Using the Lagrange planetary equations for the orbit elements of both binaries, we average over orbital timescales, obtain the equations for the secular evolution of the elements through hexadecapole order, and employ them to analyze cases of astrophysical interest. We find that, for the most part, the orbital flips found at octupole order are robust against both relativistic and hexadecapole perturbations. We show that, for equal-mass inner binaries, where the octupole terms vanish, the hexadecapole contributions can alone generate orbital flips and excursions to very large eccentricities.

Figures

Tables

Table 1. Table 1: Physical parameters and initial conditions for selected case studies

System	$m_1$	$m_2$	$m_3$	$a$ (a.u.)	$A(a.u.)$	$e$	$E$	$z$	$\omega$	${\omega }_3$
Hot Jupiters	$M_J$	$M_{\odot }$	$40M_J$	6	100	0.001	0.6	65	45	0
Coplanar Flips	$M_J$	$M_{\odot }$	$0.03M_{\odot }$	4	50	0.8	0.6	5	0	0
Asteroid-Jupiter	$0$	$M_{\odot }$	$M_J$	2	5	0.2	0.05	65	0	0
Triple star	$0.25M_{\odot }$	$M_{\odot }$	$0.6M_{\odot }$	60	800	0.01	0.6	98	0	0
CH Cygni	$0.5M_{\odot }$	$3.51M_{\odot }$	$0.909M_{\odot }$	0.05	0.21	0.32	0.6	72	145	0

Equations

$${\mathbf{a}}_a=-\frac{G{m}_b{\mathbf{x}}_{ab}}{r_{ab}^3}-\frac{G{m}_c{\mathbf{x}}_{ac}}{r_{ac}^3},$$

$${\mathbf{X}}_c$$

$${\mathbf{x}}_c$$

$$\frac{x_{a3}^j}{r_{a3}^3}=\frac{x_{c3}^j}{r_{c3}^3}-\sum _{\ell =1}^{\infty }\frac{1}{\ell !}{x}_{ac}^L{\partial }_c^{jL}\left(\frac{1}{r_{c3}}\right),$$

$$\Delta$$

$$\eta$$

$$a^j$$

$$\times n^L{N}^{⟨jL⟩}\left[{\alpha }_2^{\ell }-(-{\alpha }_1{)}^{\ell }\right],$$

$$A^j$$

$$\times n^{L+1}{N}^{⟨j(L+1)⟩}\left[{\alpha }_2^{\ell }-(-{\alpha }_1{)}^{\ell }\right],$$

$$\mathcal{E}$$

$$=\frac{1}{2}\left(m\eta v^2+M{\eta }_3{V}^2\right)-\frac{G{m}^2\eta }{r}-\frac{G{M}^2{\eta }_3}{R}$$

$$-\eta {\eta }_3\frac{G{M}^{2}r}{R^2}\sum _{\ell =1}^{\infty }\frac{(2\ell +1)!!}{(\ell +1)!}{\left(\frac{r}{R}\right)}^{\ell }$$

$$\times n^{L+1}{N}^{⟨L+1⟩}\left[{\alpha }_2^{\ell }-(-{\alpha }_1{)}^{\ell }\right],$$

$$\mathbf{J}$$

$$=m\eta \left(\mathbf{x}\times \mathbf{v}\right)+M{\eta }_3\left(\mathbf{X}\times \mathbf{V}\right),$$

$$a^j$$

$$+\frac{3}{2}\frac{G{m}_3{r}^2}{R^4}\Delta \left(5N^j{N}_n^2-2n^j{N}_n-N^j\right)$$

$$+\frac{1}{2}\frac{G{m}_3{r}^3}{R^5}(1-3\eta )(35N^j{N}_n^3-15n^j{N}_n^2$$

$$-15N^j{N}_n+3n^j),$$

$$A^j$$

$$-\frac{1}{2}\frac{GM{r}^3}{R^5}\eta \Delta (35N^j{N}_n^3-15n^j{N}_n^2$$

$$-15N^j{N}_n+3n^j)$$

$$-\frac{5}{8}\frac{GM{r}^4}{R^6}\eta (1-3\eta )(63N^j{N}_n^4-28n^j{N}_n^3$$

$$-42N^j{N}_n^2+12n^j{N}_n+3N^j),$$

$$\mathcal{E}$$

$$\displaystyle\quad-\frac{1}{2}\frac{GM^{2}\eta\eta_{3}r^{2}}{R^{3}}\biggl{[}(3N_{n}^{2}-1)+\Delta\left(\frac{r}{R}\right)N_{n}\left(5N_{n}^{2}-3\right)$$

$$\displaystyle\qquad+\frac{1}{4}(1-3\eta)\left(\frac{r}{R}\right)^{2}(35N_{n}^{4}-30N_{n}^{2}+3)\biggr{]}\,.$$

$$r$$

$$\mathbf{x}$$

$$\mathbf{n}$$

$$\mathbf{\lambda }$$

$$\mathbf{h}$$

$$R$$

$$\mathbf{X}$$

$$\mathbf{N}$$

$$\mathbf{\Lambda }$$

$$\mathbf{H}$$

$$\mathbf{J}$$

$$+J_3\left[\sin {\iota }_3\left(\sin {\Omega }_3{\mathbf{e}}_X-\cos {\Omega }_3{\mathbf{e}}_Y\right)+\cos {\iota }_3{\mathbf{e}}_Z\right],$$

$${\Omega }_3=\Omega +\pi .$$

$$\beta$$

$$z$$

$$\cot \iota =\frac{\beta +\cos z}{\sin z},\cot {\iota }_3=\frac{{\beta }^{-1}+\cos z}{\sin z}.$$

$$\frac{dp}{dt}$$

Full text

I Introduction and summary

The hierarchical three-body problem, in which a close binary system is in orbit with a distant third body, is as old as Newton’s gravity, but continues to yield surprises. The first surprise came when Newton himself failed to account for the advance of the lunar perigee caused by the perturbing effect of the distant Sun (although a correct calculation does exist in his unpublished papers). Clairaut published a correct solution in 1749. Another notable surprise was LeVerrier’s failure in 1859 to account for the advance of Mercury’s perihelion by including perturbations of the Sun-Mercury binary system due to the distant planets. The solution to this surprise was famously provided by Einstein.

A contemporary surprise was the discovery in the 1960s of the Kozai-Lidov mechanism, in which, over long timescales, there is an interchange between the eccentricity of the two-body inner orbit and its inclination relative to the plane of the third body. This remarkable effect was discovered independently by Lidov Lidov (1962), who was studying orbits of artifical satellites, and Kozai Kozai (1962) who was studying asteroid orbits. For many years, interest in the Kozai-Lidov effect was largely confined to solar-system research, until the discoveries of unusual exoplanet and multiple star systems brought the phenomenon into the astrophysical realm. Because the mechanism could generate orbits with high eccentricity, it even became of interest for general relativistic astrophysics because of the possible enhancement of relativistic effects such as the pericenter advance and the emission of gravitational radiation.

The Kozai-Lidov mechanism is obtained by expanding the perturbing acceleration in the inner binary’s motion caused by the third body in powers of $\epsilon\equiv a/A$, the ratio of the two semimajor axes, and keeping the leading term, which is proportional to $\epsilon^{3}$ (conventionally called the “quadrupole order” term). A similar expansion is performed on the acceleration of the third body. The equations of motion are averaged over time to suppress periodic effects and to reveal the long-timescale, secular changes in the orbits. One immediate consequence is that the two semimajor axes, $a$ and $A$ do not experience secular changes. In the limit where one of the inner bodies is much less massive than the other, the component of the angular momentum of the inner orbit that is perpendicular to the plane of the outer orbit turns out to be constant. Since this component is proportional to $\cos z\sqrt{1-e^{2}}$, where $z$ is the inclination angle between the normals to the two orbital planes and $e$ is the eccentricity of the inner orbit, we see that, as $z$ decreases, $e$ increases, and vice versa. The variables $e$ and $z$ oscillate between well-defined maxima and minima, depending on the initial conditions. In addition, if $z$ is initially less than $90^{\rm o}$, so that the inner orbit is prograde relative to the outer orbit, the orbit stays prograde. If the inner orbit is initially retrograde ($z>90^{\rm o}$), it stays retrograde. The inner orbit cannot “flip” relative to the outer orbit.

The next surprise came in 2011. It was known by then that in about 25 percent of exoplanet systems with “hot Jupiters”, that is Jovian-mass planets close to the host star, the planet was in a retrograde orbit relative to the spin of the star. If, as in the solar system, the star and the other planets rotated in the same direction, how did these Jupiters end up in retrograde orbits? Naoz et al. Naoz et al. (2011) pointed out that, if one included the terms in the perturbing acceleration at the next order in $\epsilon$, namely $\epsilon^{4}$, (called “octupole order” terms) then orbital flips could occur. In addition, unlike the modest variations in eccentricity allowed by the Kozai-Lidov mechanism, excursions to eccentricities very close to unity could occur. (These behaviors had actually been noticed almost a decade earlier Krymolowski and Mazeh (1999); Ford et al. (2000); Blaes et al. (2002), but at the time there was no obvious astrophysical application.) As a result a “run of the mill” Jupiter, perturbed by a more distant planet, could be flipped to a retrograde orbit and also brought very close to the star, where tidal and other dissipative processes could circularize the orbit, thus producing a retrograde “hot Jupiter”.

In follow-up papers, Naoz and collaborators Naoz et al. (2013a); Li et al. (2014) studied other situations in which orbital flips could occur. Naoz et al. Naoz et al. (2013b) studied the effects of post-Newtonian general relativistic (GR) corrections, including gravitational radiation reaction, on the generation of orbital flips and extreme eccentricities, while Liu et al. Liu et al. (2015) studied the impact of short-range forces induced by tidal, rotational and GR effects on these extreme phenomena. Lithwick and Naoz Lithwick and Naoz (2011) and Katz et al. Katz et al. (2011) studied the case where one of the inner bodies is a “test” particle. Naoz Naoz (2016) provides a thorough review of these effects in hierarchical triple systems and discusses their astrophysical implications.

Given the complexity of the hierarchical three-body problem, it is natural to ask, are there more surprises? To that end, we have gone to the next order in the expansion of the perturbing acceleration, to order $\epsilon^{5}$, called “hexadecapole order”. Other authors have addressed this level of approximation in a range of contexts, mainly using the canonical approach of Delaunay variables. Laskar and Boué Laskar and Boué (2010) obtained the disturbing function in the Hamiltonian formally to all orders and explicitly to very high orders in $\epsilon$; they did not derive the explicit equations of motion at hexadecapole order. Hamers derived the secular equations through hexadecapole order (unpublished), and Hamers and Portegies Zwart Hamers and Portegies Zwart (2016) expanded the Hamiltonian for an $N$-body system in a sequence of hierarchical orbits to hexadecapole and dotriocontupole ($\epsilon^{6}$) orders. Antognini Antognini (2015) derived (though did not display) the secular equations through hexadecapole order in both the Delaunay approach and in a method using eccentricity and angular momentum vectors, and made the code publicly available. Carvalho et al. Carvalho et al. (2016) derived the contributions to the disturbing function at hexadecapole and dotriocontupole orders, but under the assumption that the orbital plane of the third body is fixed.

We use the approach of “osculating orbit elements” whereby each of the orbits is characterized by its semimajor axis and eccentricity, its inclination and angle of ascending node relative to a reference coordinate system, and its angle of pericenter measured from the ascending node. The equations of motion for the two orbits can then be rewritten as the “Lagrange planetary” equations for the orbit elements, which take the generic form $dY^{\alpha}/dt=Q^{\alpha}(X^{\beta},Z^{\gamma},t)$, where $X^{\beta}$ and $Z^{\gamma}$ denote orbit elements of the inner and outer binary, respectively, and $Y^{\alpha}$ denotes an orbit element of either binary. We then carry out the conventional average over an orbit of both the inner binary and the outer binary, arriving at equations for the secular changes in the orbit elements. To quadrupole and octupole orders, our equations for the secular evolution of the elements agree completely with those derived using the Delaunay variables approach, and presented in Sec. 2.2 of Ford et al. Ford et al. (2000), or in Appendix A and B of Naoz et al. Naoz et al. (2013a).

We incorporate the effects of general relativity (GR) by adding to the secular equations the relativistic pericenter advances of both the inner and outer orbits at the first post-Newtonian order (we do not include additional GR terms studied in Naoz et al. (2013b)). We then apply these hexadecapole-order equations including GR to a number of case studies presented in the literature, particularly those where orbital flips and large eccentricity excursions occur at octupole order Lithwick and Naoz (2011); Naoz et al. (2013a); Li et al. (2014). We also explore the special case where the masses comprising the inner binary are equal. In this case the octupole terms vanish identically. Nevertheless, we find a “sweet spot” in the space of orbits where the hexadecapole terms alone can generate orbital flips and large eccentricity excursions.

The remainder of this paper presents details. In Sec. II we present the detailed derivation of the secular Lagrange planetary equations through hexadecapole order. In Sec. III we present five case studies analyzed using these higher-order equations. Section IV considers the equal-mass case. Section V presents concluding remarks. In an Appendix, we present a dictionary for converting between the language of osculating orbit elements and the Delaunay variables approach used in Ford et al. (2000); Naoz et al. (2013a). Hereafter, we refer to the two papers by Naoz et al. Naoz et al. (2011, 2013a) collectively as NFLRT.

II Secular evolution of hierarchical triple systems

II.1 Equations of motion and conserved quantities

We begin with the Newtonian equations of motion for a three-body system,

[TABLE]

where $a=(1,\,2,\,3)$, $b\neq c\neq a$, $G$ is Newton’s constant, $\bm{x}_{ab}\equiv\bm{x}_{a}-\bm{x}_{b}$, and $r_{ab}\equiv|\bm{x}_{ab}|$. Bodies 1 and 2 will be taken to be the “inner” binary, with body 3 taken to be the “outer” perturbing body. We define the centers of mass of the system and of the inner binary to be

[TABLE]

where $M\equiv m_{1}+m_{2}+m_{3}$ and $m\equiv m_{1}+m_{2}$. A “hierarchical” triple system is one in which the orbital separation of the inner binary is small compared to that of the outer binary, so we expand the equations of motion in powers of that small ratio. This can be carried out by writing $\bm{x}_{a3}=\bm{x}_{ac}+\bm{x}_{c3}$, where $a=(1,\,2)$, with $|\bm{x}_{ac}|\ll|\bm{x}_{c3}|$, and using the Taylor expansion

[TABLE]

where the superscript $L$ is a multi-index, with the interpretation $z^{L}\equiv z^{i_{1}}z^{i_{2}}\dots z^{i_{\ell}}$; similarly, $\partial_{c}^{jL}$ is a multi-partial derivative with respect to $\bm{x}_{c}$, and a contraction over the $\ell$ repeated indices is assumed. We now define $\bm{x}\equiv\bm{x}_{1}-\bm{x}_{2}$, $r\equiv|\bm{x}|$, $\bm{n}\equiv\bm{x}/r$, $\bm{X}\equiv\bm{x}_{3c}$, $R\equiv|\bm{X}|$, $\bm{N}\equiv\bm{X}/R$; note that $\bm{X}$ is chosen to point from the inner binary to the third body. We also define the dimensionless mass coefficients $\alpha_{i}\equiv m_{i}/m$, ($i=1,\,2$), with $\alpha_{1}+\alpha_{2}=1$. We define the dimensionless mass difference and dimensionless reduced mass

[TABLE]

and choose body 1 to be less than or equal in mass to body 2, so that $\Delta=\sqrt{1-4\eta}\geq 0$; recall that $\eta$ ranges between zero and $1/4$. We note that $\bm{x}_{1c}=\alpha_{2}\bm{x}$ and $\bm{x}_{2c}=-\alpha_{1}\bm{x}$. With this convention, in the limit $m_{1}\to 0$, $\eta\to 0$, $\Delta\to 1$, and the relative orbit $\bm{x}$ and the actual orbit ${\bm{x}}_{1c}$ of body 1 coincide. We also recall that $\partial_{c}^{jL}r_{c3}=(2\ell+1)!!N^{\langle jL\rangle}/R^{\ell+2}$, where the superscript $\langle\dots\rangle$ denotes a symmetric tracefree product (for a review see Poisson and Will (2014)).

We can then express the equation of motion of the inner binary and of the third body relative to the inner center of mass in the general form

[TABLE]

where $\bm{a}\equiv d^{2}\bm{x}/dt^{2}$ and $\bm{A}\equiv d^{2}\bm{X}/dt^{2}$. Note that the perturbing terms in the equation for $A^{j}$ depend on the inner binary’s reduced mass parameter $\eta$; this is to be expected, since in the limit in which body 1 is a test body, $\eta\to 0$, and the third body moves on an unperturbed Keplerian orbit around the massive body 2.

The equations admit conserved total energy and angular momentum, given by

[TABLE]

where $\eta_{3}\equiv mm_{3}/M^{2}$, and we have chosen the coordinates so that $\bm{X}_{c}=0$.

Beginning with $\ell=1$, the terms in the expansions over $\ell$ are conventionally denoted quadrupole, octupole, hexadecapole, dotriocontupole, etc. We expand the equations through hexadecapole order, leading to the final forms

[TABLE]

where $N_{n}\equiv\bm{N}\cdot\bm{n}$. The octupole perturbations depend on the dimensionless mass difference $\Delta$, while the hexadecapole perturbations depend on the factor $(1-3\eta)$. To hexadecapole order, the energy is given by

[TABLE]

II.2 Osculating orbit elements and the Lagrange planetary equations

We define the osculating orbit elements of the inner and outer orbits in the standard manner: for the inner orbit, we have the orbit elements $p$, $e$, $\omega$, $\Omega$ and $\iota$, with the definitions

[TABLE]

where (${\bm{e}}_{X},\,{\bm{e}}_{Y},\,{\bm{e}}_{Z}$) define a reference basis, to be specified below. From the given definitions, it is evident that ${\bm{v}}=\dot{r}{\bm{n}}+(h/r){\bm{\lambda}}$ and $\dot{r}=(he/p)\sin f$.

For the outer orbit, we have the elements $P$, $E$, $\omega_{3}$, $\Omega_{3}$, and $\iota_{3}$, with the definitions

[TABLE]

In a similar manner, ${\bm{V}}=\dot{R}{\bm{N}}+(H/R){\bm{\Lambda}}$ and $\dot{R}=(HE/P)\sin F$. The semimajor axes of the two orbits are defined by $a\equiv p/(1-e^{2})$ and $A\equiv P/(1-E^{2})$.

The total angular momentum is strictly conserved if the system is isolated and we ignore gravitational radiation reaction, therefore it is natural to choose the $Z$-axis to lie along $\bm{J}$, i.e., $\bm{e}_{Z}=\bm{J}/J$ (see Fig. 1). In general,

[TABLE]

where $J_{b}=m\eta\sqrt{Gmp}$ and $J_{3}=M\eta_{3}\sqrt{GMP}$. Thus, to align $\bm{J}$ with the $Z$-axis, we must impose $J_{b}\sin\iota\sin\Omega=-J_{3}\sin\iota_{3}\sin\Omega_{3}$ and $J_{b}\sin\iota\cos\Omega=-J_{3}\sin\iota_{3}\cos\Omega_{3}$; this implies that $\tan\Omega=\tan\Omega_{3}$ and ${\rm sgn}(\sin\Omega)=-{\rm sgn}(\sin\Omega_{3})$. Together, these imply that

[TABLE]

Another way of stating this result is that the components of the angular momenta of the two orbits in the $X-Y$ plane must be equal and opposite, and thus that the orbital planes must intersect the $X-Y$ plane along a common line, and the lines of ascending nodes must be parallel and in opposite directions. We then have that $J_{b}\sin\iota=J_{3}\sin\iota_{3}$. Defining

[TABLE]

it is straightforward to obtain the relations

[TABLE]

It will turn out that only the relative inclination $z$ between the two orbits is dynamically relevant; given an evolution for $z$ and $\beta$, the individual orbital inclinations can be recovered algebraically from Eqs. (14).

From Eqs. (7), we define the perturbing accelerations $\delta{\bm{a}}\equiv\bm{a}+Gm{\bm{n}}/r^{2}$ and $\delta{\bm{A}}\equiv\bm{A}+GM{\bm{N}}/R^{2}$. Then, for the inner binary, we define the radial $\cal R$, cross-track $\cal S$ and out-of-plane $\cal W$ components of the perturbing acceleration $\delta{\bm{a}}$, defined respectively by ${\cal R}\equiv{\bm{n}}\cdot\delta{\bm{a}}$, ${\cal S}\equiv{\bm{\lambda}}\cdot\delta{\bm{a}}$ and ${\cal W}\equiv\bm{\hat{h}}\cdot\delta{\bm{a}}$, and we write down the “Lagrange planetary equations” for the evolution of the orbit elements,

[TABLE]

The auxiliary variable $\varpi$ is defined such that the change in pericenter angle is given by $\dot{\omega}=\dot{\varpi}-\dot{\Omega}\cos\iota$.

For the outer binary, the analogous components of the perturbing acceleration $\delta{\bm{A}}$, are defined by ${\cal R}_{3}\equiv{\bm{N}}\cdot\delta{\bm{A}}$, ${\cal S}_{3}\equiv{\bm{\Lambda}}\cdot\delta{\bm{A}}$ and ${\cal W}_{3}\equiv\bm{\hat{H}}\cdot\delta{\bm{A}}$. The planetary equations for the outer binary take the form of Eqs. (15), with suitable replacements of all the relevant variables, $p\to P$, $e\to E$, $m\to M$, $f\to F$, and so on, and with $\dot{\omega}_{3}=\dot{\varpi}_{3}-\dot{\Omega}_{3}\cos\iota_{3}$.

Combining these equations and inserting the perturbing accelerations, it is straightforward to verify directly that

[TABLE]

reflecting the conservation of the three components of the total angular momentum.

II.3 Secular evolution of the orbit elements

We now use first-order perturbation theory to obtain the secular evolutions of the orbital elements. This is done by substituting constant values of the orbit elements into the right-hand sides of the planetary equations, and averaging over an orbit of both the inner and outer binaries. This is justified by the fact that the leading order perturbation is at $\epsilon^{3}$ and we are going out to order $\epsilon^{5}$. Were we including terms of dotriocontupole order ($\epsilon^{6}$) in the equations of motion, for example, it would have been necessary to invoke second-order perturbation theory for the quadrupole terms, in which the full (secular plus periodic) solutions at quadrupole order are substituted back into the Lagrange planetary equations and the orbital average carried out again.

Each planetary equation can be written in the generic form

[TABLE]

where the $X^{\beta}$ and $Z^{\beta}$ are orbit elements associated with the inner and outer binaries, respectively, and where the last step recognizes that every term on the right-hand side can be factorized into a product of terms depending only on one or the other orbital elements and on either $f$ or $F$. Then the average of $dY^{\alpha}/dt$ is approximated as product of averages of $A^{\alpha}$ and $B^{\alpha}$, in other words

[TABLE]

where the two orbital periods are given by $P_{\rm inner}=2\pi\sqrt{a^{3}/Gm}$ and $P_{\rm outer}=2\pi\sqrt{A^{3}/GM}$, with the assumption that $P_{\rm inner}\ll P_{\rm outer}$. In integrating over an orbit of the inner binary, it is useful to convert the angular variable from the true anomaly $f$, to the eccentric anomaly $u$, using the relations

[TABLE]

along with $r=a(1-e\cos u)$ and $dt=\sqrt{a^{3}/Gm}(1-e\cos u)du$. For the outer binary, we use the fact that $dt=\sqrt{P^{3}/GM}(1+E\cos F)^{-2}dF$. Thus the orbit averages may be written

[TABLE]

After carrying out the orbital averages, we convert from time $t$ to a dimensionless time scaled by the inner orbital period, namely

[TABLE]

With this scaling, the entire secular dynamics depends on the three dimensionless parameters:

[TABLE]

In terms of these parameters, the quantity $\beta=J_{b}/J_{3}$ is given by

[TABLE]

Through hexadecapole order, we recover the well-known result that $p$, $e$, $P$, and $E$ evolve in such a way that the semimajor axes $a$ and $A$ are constant, in other words

[TABLE]

The secular evolution of the remaining orbit elements is given as follows:

Quadrupole order

[TABLE]

At quadrupole order, we recover many of the features of the well-known Kozai-Lidov behavior in hierarchical triple systems, such as the oscillation of $e$ and $\iota$ as the pericenter angle $\omega$ advances. When $\eta=0$, the outer orbit is a Keplerian ellipse with constant elements, and the quantity $[a(1-e^{2})]^{1/2}\cos z$ is constant under the secular evolution of $e$ and $\iota$; this is proportional to the component of the inner orbit’s angular momentum orthogonal to the plane of the third body. For general $\eta$, there is a fixed point of the orbit elements of the inner orbit ($\dot{e}=\dot{\iota}=\dot{\omega}$), when $\omega=\pi/2$ or $3\pi/2$, and when $e$ and $z$ satisfy the constraint

[TABLE]

For the outer orbit, the fixed point implies $\dot{\iota}_{3}=\dot{E}=0$, but $\dot{\omega}_{3}\neq 0$, in general.

Octupole order

[TABLE]

It is straightforward to verify that these results are completely equivalent to those of Ford et al. Ford et al. (2000) and NFLRT Naoz et al. (2013a). In Appendix A we provide a dictionary that translates between our osculating orbits language and the Delaunay variables language used in Ford et al. (2000); Naoz et al. (2013a).

Hexadecapole order

[TABLE]

At all three orders, these equations satisfy the three constraints (16) related to the conservation of total angular momentum.

Substituting the definitions (9) and (10) into the expression (8) for the conserved energy and averaging over time, we obtain the expression

[TABLE]

The quadrupole and octupole contributions agree with the corresponding contributions to the “energy function” $F$, displayed in Eqs. (8) - (11) of Lithwick and Naoz (2011); in that calculation, $\eta$ was chosen to vanish, the constant pericenter of the outer orbit was chosen to lie along the $X$-axis, and thus $\omega_{3}=\pi-\Omega$.

III Case studies of the effects of hexadecapole contributions

We now turn to the numerical analysis of the secular evolution of the orbital elements for cases of astrophysical interest. The two semimajor axes $a$ and $A$ are constants of the motion. The precession of the nodal angle $d\Omega/d\tau$ is of no internal dynamical interest; it represents an irrelevant rotation of the entire system about the conserved total angular momentum vector. None of the evolution equations depends on $\Omega$. The equation for $d\Omega/d\tau$ is useful only for constructing the equations of evolution for $\omega$ and $\omega_{3}$ using the relations $\dot{\omega}=\dot{\varpi}-\dot{\Omega}\cos\iota$ and $\dot{\omega}_{3}=\dot{\varpi}_{3}-\dot{\Omega}\cos\iota_{3}$. The individual inclinations $\iota$ and $\iota_{3}$ can be directly linked to the relative inclination angle $z$ via Eq. (14), and only $z$ appears in the equations. Thus the dynamical system reduces to five evolution equations for the five variables $e$, $E$, $z$, $\omega$ and $\omega_{3}$, depending only on the three dimensionless parameters $\alpha=m_{3}/m$, $\eta=m_{1}m_{2}/m^{2}$ and $\epsilon=a/A$. The only place where the actual mass or distance scale enters is in the conversion from the dimensionless time $\tau$ to real time $t$ via the scaling $t=P_{\rm inner}\tau=2\pi\tau(a^{3}/Gm)^{1/2}$.

The foregoing remarks apply only in Newtonian gravity. In the real world, general relativity should be included, and indeed it is well known that the simplest quadrupole-order Kozai-Lidov oscillations can be strongly suppressed if the rate of relativistic advance of the pericenter of the inner binary is large enough Holman et al. (1997). Including the leading contribution of general relativity forces us to introduce an additional dimensionless parameter $\delta$ to the problem, given by

[TABLE]

where $c$ is the speed of light. The dominant effect is to add to the pericenter advances of the two orbits the terms

[TABLE]

Additional relativistic effects, such as those studied in Naoz et al. (2013b), will be the subject of future work.

With three fundamental parameters (four if we include general relativity) and 5 dynamical variables, a complete exploration of the full parameter space is beyond the scope of this paper. Instead we will analyze the effects of the hexadecapole contribution on a selection of case studies that have appeared in the literature. Most of these have been presented by Naoz and collaborators Naoz et al. (2011, 2013a), who first pointed out examples where orbital flips and excursions to very large eccentricities induced by octupole-order terms were astrophysically interesting. We will find that, in almost all cases, the hexadecapole and GR contributions make only small quantitative differences, but do not impact the orbital-flip or large-eccentricity phenomena.

Table 1 lists the specific parameters and initial conditions for the cases studied in this section.

III.1 Hot Jupiters

In their seminal discussion of the possibility of hot Jupiters in retrograde orbits, NFLRT considered an inner binary of a Jupiter orbiting a solar-mass star with $a=6$ a.u., perturbed by a brown-dwarf star with a mass of $40M_{J}$ and $A=100$ a.u. In this case, with $M_{\odot}=1047M_{J}$, the parameters (including the GR parameter) take the values

[TABLE]

The initial conditions chosen by NFLRT were

[TABLE]

We evolve the secular planetary equations for $1.7\times 10^{6}$ orbits of the inner binary (corresponding to about $2.5\times 10^{7}$ years) for four cases, octupole order, with and without GR precessions and hexadecapole order, with and without GR precessions. The four cases yield very similar results and so we show only two of the cases. Figure 2 shows the inclination angle $z$ and $\log(1-e)$ against time. Plotted in blue is the Newtonian evolution at octupole order, matching very well the results of Naoz et al. (2011, 2013a). Initially the system undergoes Kozai-Lidov type oscillations in $z$ but with the maximum value of $z$ rising steadily; when $z$ reaches $90^{\rm o}$, the orbit becomes retrograde and the oscillations “flip”. Later the orbit flips back to prograde, and so on. Plotted in red is the full evolution including hexadecapole terms and the GR pericenter precessions. The pattern of flips and the excursions to large eccentricity are essentially the same as in the octupole case; only the timescale has been shortened slightly, in agreement with the $N$-body integrations carried out by Naoz et al. and shown in their Fig. 3 Naoz et al. (2013a). In this case, the hexadecapole and relativistic terms do not change the behavior to any significant degree.

We remark that Carvalho et al. Carvalho et al. (2016) found that hexadecapole contributions, derived assuming a fixed third-body orbit, appeared to produce somewhat anomalous flip behavior (their Fig. 8), only to be restored to behavior consistent with direct numerical integrations by the dotriocontupole terms, derived under the same assumption (their Fig. 9). In our approach, both orbits are perturbed consistently, and the hexadecapole order results are fully compatible with the $N$-body integrations of Naoz et al. (2013a).

III.2 Orbital flips from nearly coplanar orbits

Li et al. Li et al. (2014) discovered the possibility of generating orbital flips and large eccentricities from initially nearly coplanar orbits using the octupole-order equations. The inner system was again a Jupiter-Sun binary with $a=4$ a.u., perturbed by a brown dwarf, with $m_{3}=0.03\,M_{\odot}$ and $A=50$ a.u. The parameters then have the values

[TABLE]

and the initial conditions are

[TABLE]

We evolve the equations for $2.5\times 10^{5}$ inner orbits ($2\times 10^{6}$ years), for three cases: octupole and hexadecapole orders without GR precessions, and hexadecapole order with GR precessions. The results are shown in Fig. 3. At octupole order without GR (upper panel, plotted in black), the system oscillates about small values of $z$ for a while, then migrates quickly to a retrograde orbit, oscillates about the new values for a while, then migrates back. During the transition the eccentricity reaches extreme values close to unity (lower panel). Including the hexadecapole terms shortens the timescale slightly (plotted in blue), but otherwise preserves the basic behavior. These curves are in excellent agreement with results obtained by Li Liprivate and by Hamers Hamersprivate using both $N$-body codes and orbit element codes to the same multipolar order. However, including the GR precessions with the hexadecapole terms causes the first flip to abort (plotted in red); subsequent flips are then out of phase with those where GR is not included. It is evident that the transition from prograde to retrograde orbits is very sensitive to the phases of the two pericenter angles, $\omega$ and $\omega_{3}$ as the inclination angle $z$ approaches $90^{\rm o}$. The cumulative precessions in these angles induced by general relativity can turn a transition to retrograde into a bounce back to prograde, and vice versa.

We now vary the semimajor axis $a$ of the inner orbit in order to assess the effects of GR. Holding the other parameters and initial conditions fixed, we obtain the curves shown in Fig. 4. Here the time scales as $\tau(a/4)^{3}$, where $a$ is in astronomical units; this timescale is chosen so that similar numbers of Kozai cycles can appear on one plot. For $a=5$ (blue), the pattern of flips is very similar to that obtained without GR, shown in blue in Fig. 3. For $a=4$ (red), the curve is the same as that shown in red in Fig. 3. For $a=3$ (green) the migration to large inclinations is suppressed by the more rapid GR precessions, although migrations to large eccentricities still occur. Finally, for $a=2$, (black) the GR precessions permit only small Kozai-like oscillations about the initial values of $z$ and $e$.

For the nominal value $a=4$ a.u., we also show the sensitivity of orbital flips to the pericenter angles. Figure 5 shows evolutions for four initial pericenter angles of body 3: $0^{\rm o}$ (red, same as in Fig. 3), $45^{\rm o}$ (blue), $90^{\rm o}$ (green) and $180^{\rm o}$ (black). Notice that the initial values $\omega=0^{\rm o}$, $\omega_{3}=0^{\rm o}$ correspond to orbits with the initial pericenters pointing in opposite directions (Fig. 1), while the values $\omega=0^{\rm o}$, $\omega_{3}=180^{\rm o}$ correspond to initial orbits with the pericenters pointing in the same direction. This dependence is in agreement with results from $N$-body integrations by Li Liprivate .

III.3 An asteroid Jupiter system

NFLRT showed that octupole perturbations could induce orbital flips in a Sun-asteroid-Jupiter system. In this case, $a=2$ a.u. and $A=5$ a.u., and the parameters are

[TABLE]

The initial conditions are

[TABLE]

We evolve the planetary equations for $10^{6}$ orbital periods, corresponding to about 2.8 million years, with and without hexadecapole terms. We include the GR precessions, but they turn out to have negligible effect in this example. Figure 6 shows the resulting evolutions of $z$ and $e$. Including the hexadecapole terms stretches the timescale somewhat, in agreement with the full $N$-body numerical evolutions carried out by NFLRT (see Fig. 8 of NFLRT). As in the previous example, the initial choice $\omega_{3}=180^{\rm o}$ leads to no orbital flips.

III.4 A triple-star hierarchical system

Analyzing a set of hierarchical triple-star systems studied by Fabrycky and Tremaine Fabrycky and Tremaine (2007), NFLRT again found orbital flip behavior (Fig. 9 of Naoz et al. (2013a)). The system studied consists of an inner binary with $m_{1}=0.25M_{\odot}$, $m_{2}=M_{\odot}$, $a=60$ a.u., and an outer star, with $m_{3}=0.6M_{\odot}$, $a=800$ a.u. In this example, the parameters are

[TABLE]

The initial conditions are

[TABLE]

In this example, the initial inner orbit is already retrograde. We evolve for $5\times 10^{4}$ orbits, corresponding to $2\times 10^{7}$ years. The results are shown in Fig. 7. The octupole-order curves (blue) agree well with the curves displayed in Fig. 9 of Naoz et al. (2013a), while the hexadecapole contributions (red) preserve the flips with minor changes. However, if we include the hexadecapole orders and decrease $a$, making the inner binary more relativistic, while holding the other parameters and initial conditions fixed, then the flips to prograde become progressively more sporadic, finally disappearing completely when $a=18$ a.u.

III.5 The CH Cygni system

Using the best fit parameters for the triple system CH Cygni from Mikkola and Tanikawa Mikkola and Tanikawa (1998), NFLRT showed that including the octupole-order contributions changed the evolution from conventional Kozai oscillations to orbital flips and excursions to large eccentricity. The parameters are

[TABLE]

and the initial conditions are

[TABLE]

We evolve for 4000 orbits, corresponding to about 22 years, with results shown in Fig. 8. The octupole-order results (blue) closely match those of Naoz et al. (2013a), Fig. 11, showing both orbital flips and large eccentricity excursions. But in this case, with hexadecapole contributions (red), the flips are suppressed and the eccentricity excursions are reduced. GR precessions were included, but make no discernable difference in this example.

IV Equal-mass inner binaries

When the bodies making up the inner binary have equal masses, the octupole terms vanish, leaving only the quadrupole and hexadecapole contributions. It is therefore interesting to explore whether the hexadecapole terms alone can generate orbital flips and large eccentricities. Since precise equality of masses is rare, this special case might not be of generic astrophysical interest, although it might be relevant for inner binaries consisting of neutron stars, whose masses tend to cluster around $1.4M_{\odot}$.

In the equal-mass case, $\eta=1/4$ and $\Delta=0$, and thus the free parameters reduce to three: $\alpha$, $\epsilon$ and the GR parameter $\delta$. Because the hexadecapole terms are smaller than the quadrupole terms by a factor $\epsilon^{2}$, then if $\epsilon$ is too small, hexadecapole effects are too small to be of any consequence. One “sweet spot”, where orbital flips can be induced by hexadecapole terms alone occurs around the values $\epsilon\sim 0.1$ and $\alpha\sim 10$. Note that the combination $\alpha\epsilon^{3}$, which controls the leading quadrupole effects, is still small; this constraint must hold so that the problem remains within the perturbative regime.

The first example is displayed in Fig. 9. The chosen parameters are:

[TABLE]

The initial conditions are

[TABLE]

A specific system with these parameters consists of two $1.4\,M_{\odot}$ neutron stars orbiting a $30\,M_{\odot}$ star or black hole, with $a=7$ a.u. and $A=55$ a.u. Scaling all masses and semimajor axes by a common factor $\zeta$ yields identical evolutions, since the three parameters of Eq. (42) are unchanged. Only the timescale set by the inner orbital period changes, scaling by $\zeta$. Evolving the system for 500 orbits of the inner binary, we find that the evolution for $\omega_{3}=0^{\rm o}$ (initial pericenters in opposite directions along the line of nodes) is identical to that for $\omega_{3}=180^{\rm o}$ (initial pericenters in the same direction), resulting in an orbital flip (red curves in Fig. 9), while the evolution for $\omega_{3}=90^{\rm o}$ does not show flips. This is in contrast to the cases where octupole terms dominate, where $\omega_{3}=0^{\rm o}$ leads to flips while $\omega_{3}=180^{\rm o}$ does not. This makes sense because, as can be seen from Eqs. (27), the octupole terms change sign under the transformation $\omega_{3}\to\omega_{3}+\pi$, whereas the hexadecapole terms in Eqs. (28) are invariant under that transformation. On the other hand, many pieces of the hexadecapole terms change sign under the transformation $\omega_{3}\to\omega_{3}+\pi/2$, and as a consequence, the initial angle $\omega_{3}=90^{\rm o}$ yields no flips.

In the foregoing example, the evolutions are the same whether the GR precessions are included or not. We can investigate when GR effects become important by “dialing up” the GR parameter $\delta$ while holding $\alpha$ and $\epsilon$ fixed. This is equivalent either to reducing $a$ and $A$ by the same factor, holding the masses fixed, or to increasing all the masses by the same factor, holding $a$ and $A$ fixed. We find that orbital flips are preserved until $\delta$ is about $820$ times larger than the value shown in Eq. (42).

Another example generates orbital flips from nearly coplanar orbits, an analogue of the case discussed in Sec. III.2. The results are shown in Fig. 10. In this case the parameters are

[TABLE]

A sample system is again two $1.4\,M_{\odot}$ neutron stars with $a=7$ a.u., but now orbiting a $50\,M_{\odot}$ star or black hole at $A=80$ a.u. The initial conditions are

[TABLE]

The quadrupole-order evolution, shown in blue, displays the standard Kozai-Lidov cycles, whereas the hexadecapole-order evolution shows orbital flips and excursions to extreme eccentricities, well beyond (in the sense of $\log(1-e)$) the initial relatively large initial value of $e=0.99$. Here again, GR precessions play a negligible role, suppressing the flips only when the GR parameter $\delta$ is dialed up by a factor of about $80$.

As a final example, we display in Fig. 11 the effect of slightly unequal masses on the generation of orbital flips via octupole-order terms. We again consider an inner binary of total mass $2.8\,M_{\odot}$, with $a=4$ a.u., orbiting a star or black hole of mass $50\,M_{\odot}$ at $A=50$ a.u. The initial conditions are

[TABLE]

The equal-mass case shows no orbital flips in this case (blue curves in Fig. 11), basically because $\epsilon=0.08$ is smaller than in the previous cases, and the hexadecopole terms alone are not large enough to do the job. As we change the two inner masses holding the total mass fixed, the octupole terms kick in, but are initially too small to generate flips, until we reach $m_{1}=1.33\,M_{\odot}$, $m_{2}=1.47\,M_{\odot}$, whereupon orbital flips and large eccentricities are generated (red curves).

V Concluding remarks

We have extended the study of Kozai-Lidov type hierarchical triple systems to hexadecapole order, or to order $(a/A)^{5}$, and examined a number of astrophysically interesting cases to elucidate the effect of the higher order terms on extreme behavior such as orbital flips and excursions to large eccentricity. Given the complexity of the three-body problem, even in the hierarchical regime, it may come as no surprise that we find a complicated range of behaviors. In most cases, the hexadecapole terms have only small quantitative effects on the long-term evolution of the system.

In addition, in the astrophysical systems examined in Sec. III, the parameter $\epsilon$ ranged from $0.06$ to $0.24$; at the upper end of this range, the systems are not very hierarchical. Given the inherently chaotic nature of the three-body problem, it pays to be cautious in ascribing a specific phenomenon (such as orbital flips) solely to the presence of a higher-order term, as opposed to a possible slight change in initial conditions.

For equal-mass systems (and possibly for a range of nearly equal-mass systems), where the octupole terms vanish or are suppressed, we found a region of parameter space where orbital flips and excursions to large eccentricity could be generate by the hexadecapole terms.

We have derived and presented the equations in as clear a fashion as possible, to make it easy for other researchers to use them to explore the full parameter space of hierarchical triple behavior. For example, all the examples discussed in this paper are characterized by $\beta\ll 1$, whereby the system’s angular momentum resides primarily in the outer orbit. The other limit, $\beta\gg 1$ may yield interesting behavior when hexadecapole terms are included (see Antonini and Perets (2012); Naoz et al. (2017) for a studies at octupole order). Finally one should look at the interplay between these Newtonian $N$-body effects and GR effects beyond the basic pericenter precessions, including higher PN contributions, frame dragging effects, gravitational-radiation reaction damping, and effects arising from “cross-terms” between GR and quadrupole contributions Will (2014).

Acknowledgements.

This work was supported in part by the National Science Foundation, Grant Nos.  PHY 16-00188. We are particularly grateful to Smadar Naoz, Adrian Hamers, Gongjie Li, Jean Teyssandier, Fred Rasio and Todd Thompson for valuable comments on an earlier draft of this paper.

Appendix A A Delaunay/Osculating elements dictionary

Here we provide a dictionary that may be useful in translating between the language of osculating orbit elements used in this paper, and the language of Delaunay variables used in conventional treatments of many-body dynamics, and in particular in NFLRT Naoz et al. (2011, 2013a).

NFLRT used the subscript 2 to denote the orbit elements of the outer body, whereas we use the subscript 3; they use $k^{2}$ to denote the Newtonian constant $G$. There are six Delaunay “coordinates”: the two mean anomalies $\ell_{1}$ and $\ell_{2}$, which correspond roughly to our true anomalies $f$ and $F$, the longitudes of the ascending nodes $h_{1}$, and $h_{2}$, which correspond to $\Omega$ and $\Omega_{3}$ and the arguments of pericenter $g_{1}$ and $g_{2}$, which correspond to $\omega$ and $\omega_{3}$. The “conjugate momenta” to those variables are (Eqs. (3) – (8) of Naoz et al. (2013a)):

[TABLE]

Since $G_{1}=J_{b}=m\eta\sqrt{Gmp}$ and $G_{2}=J_{3}=M\eta_{3}\sqrt{GMP}$, it is straightforward to read off the correspondences $(a_{1},\,a_{2})\rightleftharpoons(a,\,A)$, $(e_{1},\,e_{2})\rightleftharpoons(e,\,E)$, $(\iota_{1},\,\iota_{2})\rightleftharpoons(\iota,\iota_{3})$, with $z=\iota_{\rm tot}=\iota+\iota_{3}$. Note that

[TABLE]

The parameters $C_{2}$ and $C_{3}$ of Naoz et al. (2013a), Eqs. (21) and (B1), are given by

[TABLE]

The ratio

[TABLE]

consistent with Eq. (24) of Naoz et al. (2013a). The amplitudes of the perturbing effects on the elements of each orbit are controlled in Naoz et al. (2013a) by the ratios

[TABLE]

These amplitudes correspond to those displayed in Eqs. (25) and (27). Finally, in making comparisons with Naoz et al. (2013a), it is useful to note that

[TABLE]

With these translations, it can be shown that at quadrupole order, our Eqs. (25) are identical to Eqs. (A26) - (A35), and that at octupole order, our Eqs. (27) are identical to Eqs. (B6) - (B17) of Naoz et al. (2013a).