Secular Dynamics in Hierarchical Three-Body Systems

# Secular Dynamics in Hierarchical Three-Body Systems

## Abstract

The secular approximation for the evolution of hierarchical triple configurations has proven to be very useful in many astrophysical contexts, from planetary to triple-star systems. In this approximation the orbits may change shape and orientation, on time scales longer than the orbital time scales, but the semimajor axes are constant. For example, for highly inclined triple systems, the Kozai-Lidov mechanism can produce large-amplitude oscillations of the eccentricities and inclinations. Here we revisit the secular dynamics of hierarchical triple systems. We derive the secular evolution equations to octupole order in Hamiltonian perturbation theory. Our derivation corrects an error in some previous treatments of the problem that implicitly assumed a conservation of the z-component of the angular momentum of the inner orbit (i.e., parallel to the total angular momentum of the system). Already to quadrupole order, our results show new behaviors including the possibility for a system to oscillate from prograde to retrograde orbits. At the octupole order, for an eccentric outer orbit, the inner orbit can reach extremely high eccentricities and undergo chaotic flips in its orientation. We discuss applications to a variety of astrophysical systems, from stellar triples to merging compact binaries and planetary systems. Our results agree with those of previous studies done to quadrupole order only in the limit in which one of the inner two bodies is a massless test particle and the outer orbit is circular; our results agree with previous studies at octupole order for the eccentricity evolution, but not for the inclination evolution.

## 1 Introduction

Triple star systems are believed to be very common (e.g., Tokovinin, 1997; Eggleton et al., 2007). From dynamical stability arguments these must be hierarchical triples, in which the (inner) binary is orbited by a third body on a much wider orbit. Probably more than 50% of bright stars are at least double (Tokovinin, 1997; Eggleton et al., 2007). Given the selection effects against finding faint and distant companions we can be reasonably confident that the proportion is actually substantially greater. Tokovinin (1997) showed that of binary stars with period  d in which the primary is a dwarf () have at least one additional companion. He found that the fraction of triples and higher multiples among binaries with period (d) is . Moreover, Pribulla and Rucinski (2006) have surveyed a sample of contact binaries, and noted that among 151 contact binaries brighter than 10 mag., 42 are at least triple.

Many close stellar binaries with two compact objects are likely produced through triple evolution. Secular effects (i.e., coherent interactions on timescales long compared to the orbital period), and specifically Kozai-Lidov cycling (Kozai, 1962; Lidov, 1962, see below), have been proposed as an important element in the evolution of triple stars (e.g. Harrington, 1969; Mazeh and Shaham, 1979; Söderhjelm, 1982; Kiseleva et al., 1998; Fabrycky and Tremaine, 2007; Perets and Fabrycky, 2009; Thompson, 2011; Shappee and Thompson, 2012). In addition, Kozai-Lidov cycling has been suggested to play an important role in both the growth of black holes at the centers of dense star clusters and the formation of short-period binary black holes (Wen, 2003; Miller and Hamilton, 2002; Blaes et al., 2002). Recently, Ivanova et al. (2010) showed that the most important formation mechanism for black hole XRBs in globular clusters may be triple-induced mass transfer in a black hole-white dwarf binary.

Secular perturbations in triple systems also play an important role in planetary system dynamics. Kozai (1962) studied the effects of Jupiter’s gravitational perturbation on an inclined asteroid in our own solar system. In the assumed hierarchical configuration, treating the asteroid as a test particle, Kozai (1962) found that its inclination and eccentricity fluctuate on timescales much larger than its orbital period. Jupiter, assumed to be in a circular orbit, carries most of the angular momentum of the system. Due to Jupiter’s circular orbit and the negligible mass of the asteroid, the system’s potential is axisymmetric and thus the component of the inner orbit’s angular momentum along the total angular momentum is conserved during the evolution. Kozai (1979) also showed the importance of secular interactions for the dynamics of comets (see also Quinn et al., 1990; Bailey et al., 1992; Thomas and Morbidelli, 1996). The evolution of the orbits of binary minor planets is dominated by the secular gravitational perturbation from the sun (Perets and Naoz, 2009); properly accounting for the resulting secular effects—including Kozai cycling—accurately reproduces the binary minor planet orbital distribution seen today (Naoz et al., 2010; Grundy et al., 2011). In addition Kinoshita and Nakai (1991), Vashkov’yak (1999), Carruba et al. (2002), Nesvorný et al. (2003), Ćuk and Burns (2004) and Kinoshita and Nakai (2007) suggested that secular interactions may explain the significant inclinations of gas giant satellites and Jovian irregular satellites.

Similar analyses have been applied to the orbits of extrasolar planets (e.g., Innanen et al., 1997; Wu and Murray, 2003; Fabrycky and Tremaine, 2007; Wu et al., 2007; Naoz et al., 2011; Veras and Ford, 2010; Correia et al., 2011). Naoz et al. (2011) considered the secular evolution of a triple system consisting of an inner binary containing a star and a Jupiter-like planet at several AU, orbited by a distant Jupiter-like planet or brown-dwarf companion. Perturbations from the outer body can drive Kozai-like cycles in the inner binary, which, when planet-star tidal effects are incorporated, can lead to the capture of the inner planet onto a close, highly-inclined or even retrograde orbit, similar to the orbits of the observed retrograde “hot Jupiters.” Many other studies of exoplanet dynamics have considered similar systems, but with a very distant stellar binary companion acting as perturber. In such systems, the outer star completely dominates the orbital angular momentum, and the problem reduces to test-particle evolution (see Lithwick and Naoz, 2011; Katz et al., 2011; Naoz et al., 2012a). If the lowest level of approximation is applicable (e.g., the outer perturber is on a circular orbit), the -component of the inner orbit’s angular momentum is conserved (e.g., Lidov and Ziglin, 1974).

In early studies of high-inclination secular perturbations (Kozai, 1962; Lidov, 1962), the outer orbit was circular and again dominated the orbital angular momentum of the system. In this situation, the component of the inner orbit’s angular momentum along the z-axis is conserved. In many later studies the assumption that the -component of the inner orbit’s angular momentum is constant was built into the equations (e.g. Eggleton et al., 1998; Mikkola and Tanikawa, 1998; Zdziarski et al., 2007). In fact these studies are only valid in the limit of a test particle forced by a perturber on a circular orbit. To leading order in the ratio of semimajor axes, the double averaged potential of the outer orbit is axisymmetric (even for an eccentric outer perturber), thus if taken to the test particle limit, this results in a conservation of the -component of the inner orbit’s angular momentum. We refer to this limit as the “standard” treatment of Kozai oscillations, i.e. quadrupole-level approximation in the test particle limit (test particle quadrupole, hereafter TPQ).

In this paper we show that a common mistake in the Hamiltoniano treatment of these secular systems can lead to the erroneous conclusion that the -component of the inner orbit’s angular momentum is constant outside the TPQ limit; in fact, the -component of the inner orbit’s angular momentum is only conserved by the evolution in the test-particle limit and to quadrupole order. To demonstrate the error we focus on the quadrupole (non-test-particle) approximation in the main body of the paper, but we include the full octupole–order equations of motion in an appendix.

In what follows we show the applications of these two effects (i.e., correcting the error and including the full octupole–order equations of motion) by considering different astrophysical systems. Note that the applications illustrated in the text are inspired by real systems; however, we caution that we consider here only Newtonian point mass dynamics, while in reality other effects such as tides and general relativity can greatly effect the evolution. For example, general relativity may alter the evolution of the system, which can give rise to a resonant behavior of the inner orbitÕs eccentricity (e.g., Ford et al., 2000a; Naoz et al., 2012b). Furthermore, tidal forces can suppress the eccentricity growth of the inner orbit, and thus significantly modify the evolutionary track of the system (e.g., Mazeh and Shaham, 1979; Söderhjelm, 1984; Kiseleva et al., 1998). In particular tides, in some cases, can considerably suppress the chaotic behavior that arises in the presence of the the octupoleÐlevel of approximation (e.g., Naoz et al., 2011, 2012a). Therefore, while the examples presented in this paper are inspired by real astronomical systems, the true evolutionary behavior will be modified from what we show once the eccentricity becomes too high.

This paper is organized as follows. We first present the general framework (§2); we then derive the complete formalism for the quadrupole-level approximation and the equations of motion (§3), we also develop the octupole-level approximation equations of motion in §4. We discuss a few of the most important implications of the correct formalism in §5. We also compare our results with those of previous studies (§5) and offer some conclusions in §6.

## 2 Hamiltonian Perturbation Theory for Hierarchical Triple Systems

Many gravitational triple systems are in a hierarchical configuration—two objects orbit each other in a relatively tight inner binary while the third object is on a much wider orbit. If the third object is sufficiently distant, an analytic, perturbative approach can be used to calculate the evolution of the system. In the usual secular approximation (e.g., Marchal, 1990), the two orbits torque each other and exchange angular momentum, but not energy. Therefore the orbits can change shape and orientation (on timescales much longer than their orbital periods), but not semimajor axes (SMA).

We first define our basic notations. The system consists of a close binary (bodies of masses and ) and a third body (mass ). It is convenient to describe the orbits using Jacobi coordinates (Murray and Dermott, 2000, p. 441-443). Let be the relative position vector from to and the position vector of relative to the center of mass of the inner binary (see fig. 1). Using this coordinate system the dominant motion of the triple can be divided into two separate Keplerian orbits: the relative orbit of bodies 1 and 2, and the orbit of body 3 around the center of mass of bodies 1 and 2. The Hamiltonian for the system can be decomposed accordingly into two Keplerian Hamiltonians plus a coupling term that describes the (weak) interaction between the two orbits. Let the SMAs of the inner and outer orbits be and , respectively. Then the coupling term in the complete Hamiltonian can be written as a power series in the ratio of the semi-major axes (e.g., Harrington, 1968). In a hierarchical system, by definition, this parameter is small.

The complete Hamiltonian expanded in orders of is (e.g., Harrington, 1968),

 H = k2m1m22a1+k2m3(m1+m2)2a2 +k2a2∞∑j=2αjMj(r1a1)j(a2r2)j+1Pj(cosΦ) ,

where is the gravitational constant, are Legendre polynomials, is the angle between and (see Figure 1) and

 Mj=m1m2m3mj−11−(−m2)j−1(m1+m2)j . (2)

Note that we have followed the convention of Harrington (1969) and chosen our Hamiltonian to be the negative of the total energy, so that for bound systems.

We adopt the canonical variables known as Delaunay’s elements, which provide a particularly convenient dynamical description of our three-body system (e.g. Valtonen and Karttunen, 2006). The coordinates are chosen to be the mean anomalies, and , the longitudes of ascending nodes, and , and the arguments of periastron, and , where subscripts denote the inner and outer orbits, respectively. Their conjugate momenta are

 L1 = m1m2m1+m2√k2(m1+m2)a1 , (3) L2 = m3(m1+m2)m1+m2+m3√k2(m1+m2+m3)a2 ,
 G1=L1√1−e21 ,G2=L2√1−e22 , (4)

and

 H1=G1cosi1 ,H2=G2cosi2 , (5)

where () is the inner (outer) orbit eccentricity. Note that and are also the magnitudes of the angular momentum vectors ( and ), and and are the -components of these vectors. Figure 2 shows the resulting configuration of theses vectors. The following geometric relations between the momenta follow from the law of cosines:

 cositot = G2tot−G21−G222G1G2 , (6) H1 = G2tot+G21−G222Gtot , (7) H2 = G2tot+G22−G212Gtot , (8)

where is the (conserved) total angular momentum, and the angle between and defines the mutual inclination . From eqs. (7) and (8) we find that the inclinations and are determined by the orbital angular momenta:

 cosi1 = G2tot+G21−G222GtotG1 , (9) cosi2 = G2tot+G22−G212GtotG2 . (10)

In addition to these geometrical relations we also have that

 H1+H2=Gtot=const . (11)

The canonical relations give the equations of motion:

 dLjdt=∂H∂lj ,dljdt=−∂H∂Lj , (12) dGjdt=∂H∂gj ,dgjdt=−∂H∂Gj , (13) dHjdt=∂H∂hj ,dhjdt=−∂H∂Hj , (14)

where . Note that these canonical relations have the opposite sign relative to the usual relations (e.g., Goldstein, 1950) because of the sign convention we have chosen for our Hamiltonian. Finally we write the Hamiltonian through second order in as (e.g., Kozai, 1962)

 H = β12L21+β22L22+ 4β3(L41L62)(r1a1)2(a2r2)3(3cos2Φ+1) ,

where the mass parameters are

 β1 = k2m1m2L21a1 , (16) β2 = k2(m1+m2)m3L22a2 (17) and β3 = k416(m1+m2)7m73(m1m2)3(m1+m2+m3)3 . (18)

## 3 Secular Evolution to the Quadrupole Order

In this section, we derive the secular quadrupole–level Hamiltonian. in Appendix A we develop the complete quadrupole-level secular approximation and in particular in Appendix A.3 we present the quadrupole–level equations of motion. The main difference between the derivation shown here (see also Appendix A) and those of previous studies lies in the “elimination of nodes” (e.g., Kozai, 1962; Jefferys and Moser, 1966). This is related to the transition to a coordinate system with the total angular momentum along the z-axis, which is known as the invariable plane (e.g., Murray and Dermott, 2000). In this coordinate system (see Figure 2), the longitudes of the ascending nodes differ by , i.e.,

 h1−h2=π . (19)

Conservation of the total angular momentum implies that this relation holds at all times. Many previous works have exploited it to explicitly simplify the Hamiltonian by setting before deriving the equations of motion. After the substitution, the Hamiltonian is independent of the longitudes of ascending nodes ( and ), and this can lead to the incorrect conclusion that when the canonical equations of motion are derived. Some previous studies incorrectly concluded that the -components of the orbital angular momenta are always constant (see also Appendix C). The substitution is incorrect at the Hamiltonian level because it unduly restricts variations in the trajectory of the system to those where . After deriving the equations of motion, however, we can exploit the relation , which comes from the conservation of angular momentum. This considerably simplifies the evolution equations. We show (Appendices A.3 and B) that one can still use the Hamiltonian with the nodes eliminated found in previous studies (e.g., Kozai, 1962; Harrington, 1969) as long as the evolution equations for the inclinations are derived from the total angular momentum conservation, instead of using the canonical relations. Of course, the correct evolution equations can also be calculated from the correct Hamiltonian (without the nodes eliminated), which we derive in this section.

We note that there are some other derivations of the secular evolution equations that avoid the elimination of the nodes (Farago and Laskar, 2010; Laskar and Boué, 2010; Mardling, 2010; Katz and Dong, 2011), and thus do not suffer from this error1.

The secular Hamiltonian is given by the average over the rapidly-varying and in equation (2) (see Appendix A for more details)

 H2 = C28{[1+3cos(2i2)]([2+3e21][1+3cos(2i1)] + 30e21cos(2g1)sin2(i1))+3cos(2Δh)[10e21cos(2g1) × {3+cos(2i1)}+4(2+3e21)sin(i1)2]sin2(i2) + 12{2+3e21−5e21cos(2g1)}cos(Δh)sin(2i1)sin(2i2) + 120e21sin(i1)sin(2i2)sin(2g1)sin(Δh) − 120e21cos(i1)sin2(i2)sin(2g1)sin(2Δh)} ,

where

 C2=k416(m1+m2)7(m1+m2+m3)3m73(m1m2)3L41L32G32 . (21)

Making the usual (incorrect) substitution (i.e. eliminating the nodes), we get the quadrupole-level Hamiltonian that has appeared in many previous works (see, e.g. Ford et al., 2000a):

 H2(Δh→π) = C2{(2+3e21)(3cos2itot−1) + 15e21sin2itotcos(2g1)} ,

where we have set . Because this Hamiltonian is missing the longitudes of ascending nodes ( and ), many previous studies concluded that the -components (i.e. vertical components) of the angular momenta of the inner and outer orbits (i.e., and ) are constants.

We derive the quadrupole–level equations of motions in Appendix A.3. In particular, we give the equations of motion of the z-component of the angular momentum of the inner and outer orbits derived from the Hamiltonian in Eq. (3). As we show in the subsequent sections the evolution of produces a qualitatively different evolutionary route for many astrophysical systems considered in previous works.

In Appendix A.4 we show that the quadrupole approximation leads to well-defined minimum and maximum eccentricity and inclination. The eccentricity of the inner orbit and the inner (and mutual) inclination oscillate. In the test-particle limit, our formalism gives the critical initial mutual inclination angles for large oscillations of with nearly-zero initial inner eccentricity, in agreement with Kozai (1962).

It is easy to show that and are constant only in the TPQ limit without using the explicit equations of motion in Appendix A. Because the Hamiltonian in Eq. (3) is independent of , at the quadrupole level. Combining this with the geometric relation in Eq. (7), , and the constantcy of the total angular momentum, , we have that

 ˙H1=˙G1G1Gtot. (23)

In the TPQ limit, , so , and the -component of each orbit’s angular momentum is conserved. Outside this limit, when is not negligable, and cannot be constant. Note that the TPQ limit, where , is equivalent to the limit where appearing in many previous works.

## 4 Octupole-Level Evolution

In Appendix B, we derive the secular evolution equations to octupole order. Many previous octupole–order derivations provided correct secular evolution equations for at least some of the elements, in spite of using the elimination of nodes substitution at the Hamiltonian level (e.g. Harrington, 1968, 1969; Sidlichovsky, 1983; Krymolowski and Mazeh, 1999; Ford et al., 2000a; Blaes et al., 2002; Lee and Peale, 2003; Thompson, 2011). This is because the evolution equations for , , and can be found correctly from a Hamiltonian that has had and eliminated by the relation ; the partial derivatives with respect to the other coordinates and momenta are not affected by the substitution. The correct evolution of and can then be derived, not from the canonical relations, but from total angular momentum conservation. We discuss in more details the comparison between this work and previous analyses in §5.

The octupole-level terms in the Hamiltonian can become important when the eccentricity of the outer orbit is non-zero, and if is large enough. We quantify this by considering the ratio between the octupole to quadrupole-level coefficients, which is

 C3C2=154(m1−m2m1+m2)(a1a2)11−e22 , (24)

where is the octupole-level coefficient [eq. (71)] and is the quadrupole-level coefficient [eq. (21)]. We define

 ϵM=(m1−m2m1+m2)(a1a2)e21−e22 , (25)

which gives the relative significance of the octupole-level term in the Hamiltonian. This parameter has three important parts; first the eccentricity of the outer orbit (), second, the mass difference of the inner binary ( and ) and the SMA ratio2. In the test particle limit (i.e., ) is reduced to the octupole coefficient introduced in Lithwick and Naoz (2011) and Katz et al. (2011),

 ϵ=(a1a2)e21−e22 . (26)

We call the octupole-level behavior of a system for which is not satisfied the “eccentric Kozai-Lidov” (EKL) mechanism.

The octupole terms vanish when . Therefore if one artificially held , in the test-particle limit the inner body’s orbit would be given by the equations derived by Kozai (1962), i.e. by the test particle quadrupole equations. However, at octupole order the value of evolves in time if the inner body is massive. Furthermore, even if the inner body is massless, if the outer body has then the inner body’s behavior will also be different than in Kozai’s treatment. For example, Lithwick and Naoz (2011) and Katz et al. (2011) find that the inner orbit can flip orientation (see below) even in the test-particle, octupole limit. The octupole-level effects can change qualitatively the evolution of a system. Compared to the quadrupole-level behavior, the eccentricity of the inner orbit in the EKL mechanism can reach a much higher value. In some cases these excursions to very high eccentricities are accompanied by a “flip” of the orbit with respect to the total angular momentum, i.e. starting with the inner orbit can eventually reach (see Figures 69 for examples). Chaotic behavior is also possible at the octupole level (Lithwick and Naoz, 2011), but not at the quadrupole–level.

Given the large, qualitative changes in behavior moving from quadrupole to octupole order in the Hamiltonian, is it possible that similar changes in the secular evolution may occur at even higher orders? The answer to this question probably lies in the elimination of as an integral of motion at octupole order, leaving only four integrals of motion: the energy of the system, and the three components of the total angular momentum. There are no more integrals of motion to be eliminated, and thus one might expect no more dramatic changes in the evolution when moving to even higher orders. It is possible to see this quantitatively for specific initial conditions through comparisons with direct -body integrations. We compared our octupole equations with direct -body integrations, using the Mercury software package (Chambers and Migliorini, 1997). We used both Burlisch-Stoer and symplectic integrators (Wisdom and Holman, 1991) and found consistent results between the two. We present the results of a typical integration compared to the integration of the octupole-level secular equations in Figure 3. The initial conditions (see caption) for this system are those of Naoz et al. (2011), Figure 1. We find good agreement between the direct integration and the secular evolution at octupole order. Both show a beat-like pattern of eccentricity oscillations, suggesting an interference between the quadrupole and octupole terms, and both methods show similar flips of the inner orbit.

## 5 Implications and Comparison with Previous Studies

The Kozai (1962) and Lidov (1962) equations of motions are correct to quadrupole order and for a test particle, but differ from the correct evolution equations for non-test-particle inner orbits and/or at octupole order. In this Section we show how these differences give rise to qualitatively different evolutionary behaviors than those assumed in some previous works.

### 5.1 Massive Inner Object at the Quadrupole Level

The danger with working in the wrong limit is apparent if we consider an inner object that is more massive then the outer object. While the TPQ formalism incorrectly assumes that the orbit of the outer body is fixed in the invariable plane, and therefore the inner body’s vertical angular momentum is constant, the quadrupole-level equations presented in Appendix A.3 do not.

We compare the two formalisms in Figure 4. We consider the triple system PSR B162026 located near the core of the globular cluster M4. The inner binary contains a millisecond radio pulsar of and a companion of (McKenna and Lyne, 1988). Following Ford et al. (2000b), we adopt parameters for the outer perturber of and . Note that Ford et al. (2000b) found , but it is interesting to show that even for an axisymmetric outer potential the evolution of the system is qualitatively different then the TPQ approximation (see the caption for a full description of the initial conditions). Note that the actual measured inner binary eccentricity is , however in order to illustrate the difference we adopt a higher value ( ). For these initial conditions , so a careful analysis would require incorporating the octupole–order terms in the motion; nevertheless, we consider the evolution of the system to quadrupole order for comparison with the TPQ formalism. We have verified, however, that the neglected octupole–order effects do not qualitatively change the behavior of the system. This is because the outer companion mass is low, and hence the inner orbit does not exhibit large amplitude oscillations3.

For the comparison, we do not compare the (constant) from the TPQ formalism to the (varying) of the correct formalism. Instead, we compare the (varying) from the correct formalism (solid red line) with (dashed blue line), which is the vertical angular momentum that would be inferred in our formalism if the outer orbit were instantaneously in the invariable plane, as assumed in the TPQ formalism.

In Figure 4, the mutual inclination oscillates between to , and thus crosses . These oscillations are mostly due to the oscillations of the outer orbit’s inclination, while does not change by more than in each cycle. Clearly, the outer orbit does not lie in the fixed invariable plane. Figure 4, bottom panel, shows , which, in the TPQ limit, is the vertical angular momentum of the inner body.

We can evaluate analytically the error introduced by the application of the TPQ formalism to this situation. We compare the vertical angular momentum () as calculated here to . The relative error between the formalisms is . In Figure 5 we show the ratio between the inner orbit’s vertical angular momentum in the TPQ limit (i.e., ) and equation (52) as a function of the total angular momentum ratio, , for various inclinations. Note that this error can be calculated without evolving the system by using angular momentum conservation, Eq. (6). The TPQ limit is only valid when .

### 5.2 Octupole–Level Planetary Dynamics

Recent measurements of the sky-projected angle between the orbits of several hot Jupiters and the spins of their host stars have shown that roughly one in four is retrograde (Gaudi and Winn, 2007; Triaud et al., 2010; Albrecht et al., 2012). If these planets migrated in from much larger distances through their interaction with the protoplanetary disk (Lin and Papaloizou, 1986; Masset and Papaloizou, 2003), their orbits should have low eccentricities and inclinations4. Disk migration scenarios therefore have difficulty accounting for the observed retrograde hot Jupiter orbits. An alternative migration scenario that can account for the retrograde orbits is the secular interaction between a planet and a binary stellar companion (Wu and Murray, 2003; Fabrycky and Tremaine, 2007; Wu et al., 2007; Takeda et al., 2008; Correia et al., 2011). For an extremely distant and massive companion () the quadrupole test-particle approximation applies, and is nearly constant (where the planet is the massless body). Although this forbids orbits that are truly retrograde (with respect to the total angular momentum of the system), if the inner orbit begins highly inclined relative to the outer star’s orbit and aligned with the spin of the inner star, then the star-planet spin-orbit angle can change by more than during the secular evolution of the system, producing apparently retrograde orbits (Fabrycky and Tremaine, 2007; Correia et al., 2011). Nonetheless, a difficulty with this “stellar Kozai” mechanism is that even with the most optimistic assumptions it can only produce of hot Jupiters (Wu et al., 2007).

Wu and Murray (2003), Wu et al. (2007), Fabrycky and Tremaine (2007) and Correia et al. (2011) studied the evolution of a Jupiter-mass planet in stellar binaries in the TPQ formalism. For example, the case of HD 80606b (Wu and Murray (2003, Fig. 1); Fabrycky and Tremaine (2007, Fig. 1) and Correia et al. (2011, also Fig. 1)) was considered with an outer stellar companion at  AU. However, if the companion is assumed to be eccentric is not negligible, and the system is more appropriately described with the test particle octupole–level approximation (e.g., Lithwick and Naoz, 2011; Katz et al., 2011). Furthermore, the statistical distribution for closer stellar binaries in Wu et al. (2007) and Fabrycky and Tremaine (2007) is only valid in the approximation where the outer orbit’s eccentricity is zero. In fact, for the systems considered in those studies is not negligible and the octupole–level approximation results in dramatically different behavior as was shown in Naoz et al. (2012a). The same dramatic difference in behavior also exists in the analysis of triple stars (e.g., Fabrycky and Tremaine, 2007; Perets and Fabrycky, 2009), see §5.4.

A dramatic difference between the octupole and quadrupole–level of approximation is that the former often generates extremely high eccentricities. In real systems, such high eccentricities can be suppressed by tides or GR (e.g., Söderhjelm, 1984; Eggleton et al., 1998; Kiseleva et al., 1998; Borkovits et al., 2004). Flips can also be prevented because they typically occur shortly after extreme eccentricities (see Teyssandier et al. in prep.). In our previous studies that include tides, planetary perturbers typically allow flips to happen, while stellar perturbers mostly suppress them (Naoz et al., 2011, 2012a) But in both cases, tides quantitatively affect the evolution.

Naoz et al. (2011) considered planet-planet secular interactions with tidal interactions as a possible source of retrograde hot Jupiters. In this situation is not small, requiring computation of the octupole-level secular dynamics. In Figures 6 and 7 we show the evolution of a representative configuration (see the caption for a full description of the initial conditions). For this configuration, . Flips of the inner orbit are associated with evolution to very high eccentricity (see Figures 6 and 7).

### 5.3 Octupole–Level Solar System Dynamics

Kozai (1962) studied the dynamical evolution of an asteroid due to Jupiter’s secular perturbations. He assumed that Jupiter’s eccentricity is strictly zero. However, Jupiter’s eccentricity is , and thus studying the evolution of a test particle in the asteroid belt ( AU) places the evolution in a regime where the eccentric Kozai-Lidov effect could be significant, with (Lithwick and Naoz, 2011; Katz et al., 2011).

We considered the evolution of asteroid at  AU (assumed to be a test particle) due to Jupiter at  AU with eccentricity of (see the caption for a full description of the initial conditions). The asteroid is a test particle and therefore . In Figure 8 we compare the evolution of an asteroid using the TPQ limit (e.g., Kozai, 1962; Thomas and Morbidelli, 1996; Kinoshita and Nakai, 2007) and the octupole-level evolution discussed here. For this value of , the eccentric Kozai-Lidov effect significantly alters the evolution of the asteroid, even driving it to such high inclination that the orbit becomes retrograde. Though we deal only with point masses in this work, note that the eccentricity is so high that the inner orbit’s pericenter lies well within the sun.

The value of here is mainly due to the relative high in the problem (an issue raised in the original work on this problem (Kozai, 1962)). The system is very packed which raises questions with regards to the validity of the hierarchical approximation. Even in the EKL formalism, such high eccentricities occur that the asteriod collides with the sun and the apo-center of the asteroid approaches about  AU from Jupiter’s orbit. To determine the importance of these effects, we ran an -body simulation using the Mercury software package (Chambers and Migliorini, 1997). We used both Bulirsch-Stoer and symplectic integrators (Wisdom and Holman, 1991). The results are depicted at Figure 8, which show that the TPQ limit is indeed inadequate for the system. In addition the octupole–level approximation has some deviations from the direct -body integration, particularly in the high eccentricity regime. Note that the evolution of the asteroid in the direct integration resulted in a collision with the Sun5. In reality, it is likely that a planetary encounter would remove the asteroid from the solar system before this point. In contrast to the EKL mechanism, assuming zero eccentricity for Jupiter results in consistent results between the secular evolution and the direct integration (Thomas and Morbidelli, 1996).

As shown in Figure 8, taking into account Jupiter’s eccentricity (), produces a dramatically different evolutionary behavior, including retrograde orbits for the asteroid. Thomas and Morbidelli (1996) applied the TPQ formalism to the asteroid-Jupiter setting (see for example their Figure 2 for  AU). Kinoshita and Nakai (2007) developed an analytical solution for the TPQ limit (see also Kinoshita and Nakai, 1991, 1999).

The TPQ formalism has also been applied to the study of the outer solar system. Kinoshita and Nakai (2007) applied their analytical solution to Neptune’s outer satellite Laomedeia. This system has and thus the TPQ limit there is justified. In addition, Perets and Naoz (2009) have studied the evolution of binary minor planets using the TPQ approximation. In this problem and thus the TPQ approximation is valid.

Lidov and Ziglin (1976, sections 3–4) also solved analytically the quadrapole–level approximation but, unlike Kinoshita and Nakai (2007), they did not restrict themselves to the TPQ limit, and used the total angular momentum conservation law in order to calculate the inclinations. Thus, their formalism is equivalent to ours at quadrupole–order. Later, Mazeh and Shaham (1979) also derived evolution equations outside the TPQ limit (their eqs. A1-A8), allowing for small eccentricities and inclinations of the outer body.

### 5.4 Octupole–Level Perturbations in Triple Stars

The evolution of triple stars has been studied by many authors using the standard (TPQ) formalism (e.g., Mazeh and Shaham, 1979; Eggleton et al., 1998; Kiseleva et al., 1998; Mikkola and Tanikawa, 1998; Eggleton and Kiseleva-Eggleton, 2001; Fabrycky and Tremaine, 2007; Perets and Fabrycky, 2009). In some cases the corrected formalism derived here can give rise to qualitatively different results. We show that some of the previous studies should be repeated in order to account for the correct dynamical evolution, and give one example where the eccentric Kozi-Lidov mechanism dramatically changes the evolution.

Fabrycky and Tremaine (2007) studied the distribution of triple star properties using Monte Carlo simulations. We choose a particular system from their triple-star suite of simulations to illustrate how the dynamics including the octupole order can be qualitatively different from what would be seen at quadrupole order (see the caption for a full description of the initial conditions). For this system (and ). The evolution of the system is shown in Figure 9. At octupole order, the inclination of the inner orbit oscillates between about and , often becoming retrograde (relative to the total angular momentum), while the quadrupole–order behavior is very different and the inner orbit remains always prograde. The octupole–order treatment also gives rise to much higher eccentricities (Krymolowski and Mazeh, 1999; Ford et al., 2000a). In Figure 10 we compare the octupole–level evolution (of the same system) with direct 3–body integration.

The evolution shown in Figure 9 is for point-mass stars; in reality, these high-eccentricity excursions would actually drive the inner binary to its Roche limit, leading to mass transfer. For these high eccentricities tides will play an important role and thus in reality flips in similar systems may be suppressed. Similarly the high eccentricities often excited through the eccentric Kozai mechanism can also lead to compact object binary merger.

The possibility of forming blue stragglers through secular interactions in triple star systems has been suggested by Perets and Fabrycky (2009) and Geller et al. (2011). As shown in Krymolowski and Mazeh (1999); Ford et al. (2000a) and in the example above the minimum pericenter distance of the inner binary can differ significantly between the TPQ and EKL formalisms. This suggests that using the correct EKL formalism could significantly increase the computed likelihood of such a formation mechanism for blue stragglers.

For many years CH Cygni was considered to be an interesting triple candidate because it exhibits two clear distinguishable periods (e.g. Donnison and Mikulskis, 1995; Skopal et al., 1998; Mikkola and Tanikawa, 1998; Hinkle et al., 1993). However, a triple system model based on the TPQ Kozai mechanism (Mikkola and Tanikawa, 1998) did not reproduce the observed masses of the system (Hinkle et al., 1993, 2009). Applying the corrected formalism in this paper to the system parameters derived in Mikkola and Tanikawa (1998) gives a very different evolution than in the TPQ formalism6. Therefore, it seems likely that an analysis based on the formalism discussed in this paper would give a significantly different fit. In Figure 11 we illustrate the differences between the TPQ, correct quadrupole, and octupole evolution of the system. The best-fit parameters of the system are taken from Mikkola and Tanikawa (1998) where (see the caption for a full description of the initial conditions, where we allowed for a freedom in our choice of and since the best fit was found using the TPQ limit, at which is fixed). Note that the choice of the inner eccentricity does not strongly influence the evolution while the choice of the outer orbit’s eccentricity does. Most importantly, the rather large for this system implies that the system is not stable, i.e., the averaging over the orbits is not justified. From direct integration we found that the system undergoes strong encounters and the inner binary collides in this example.

It is also interesting to investigate a system for which the eccentric Kozai mechanism is suppressed due to comparable masses for the inner orbit, and low eccentricity of the outer orbit (i.e., ). Kiseleva et al. (1998) and Eggleton and Kiseleva-Eggleton (2001) studied the Algol triple system (Lestrade et al., 1993) using the TPQ equations. The TPQ equations were also used in the paper that introduced the influential KCTF mechanism (Mazeh and Shaham, 1979; Eggleton et al., 1998). Note that tides dominate the evolution of the Algol system today (e.g.,  Söderhjelm, 2006). Figure 12 compares the evolution computed in the (incorrect) TPQ formalism, the correct quadrupole formalism, and the octupole-level EKL formalism applied to an Algol–like system. The correct quadrupole formalism decreases the minimum value of by almost a factor of 2 relative to the TPQ formalism. The reduced pericenter distance would strongly increase the effects of tidal friction (not included here), which may lead to rapid circularization of the inner orbit. The octupole-level computation decreases the minimum pericenter distance by a further 40%.

Note that the masses and orbital parameters used in Kiseleva et al. (1998) and Eggleton and Kiseleva-Eggleton (2001) are out of date. New observations (e.g., Baron et al., 2012) find the secondary mass to be smaller then the primary and the mutual inclination to be closer to . In Figure 13 we show the octupole–level evolution of the system considering the new parameters. In the absence of any additional physical mechanism, such as general relativity, tides, mass transfer, etc., the EKL mechanism could play a very important role in the dynamical evolution of the system.

We note that the inner binary in the Algol system is dominated by tidal effects (Söderhjelm, 1975; Kiseleva et al., 1998; Eggleton and Kiseleva-Eggleton, 2001) and figures 12 and 13 do not represent the system today but an Algol–like analogy. We use the Algol parameters here only to show hypothetical outcomes of the correct dynamical evolution. It would be interesting to study stellar evolution including tides in the context of the EKL mechanism, for a system such as Algol.

### 5.5 The Danger of the Quadrupole–Level of Approximation

The octupole-level Hamiltonian and equations of motion were previously derived by Harrington (1968, 1969); Sidlichovsky (1983); Marchal (1990); Krymolowski and Mazeh (1999); Ford et al. (2000a); Blaes et al. (2002) and Lee and Peale (2003). Most of the equations of motion can be derived correctly when applying the elimination of the nodes—only the and equations are affected. These authors calculated the time evolution of the inclinations (i.e.  and ) from the total (conserved) angular momentum, and thus avoided the problem that arises when eliminating the nodes from the Hamiltonian. In appendix B we show the complete set of equations of motion for the octupole-level approximation, derived from a correct Hamiltonian, including the nodal terms.

As displayed here the octupole-level approximation gives rise to a qualitatively different evolutionary behavior for cases where [see eq. (25)] is not negligible. We note that many previous studies applied the quadrupole-level approximation, which may lead to significantly different results (e.g., Mazeh and Shaham, 1979; Quinn et al., 1990; Bailey et al., 1992; Innanen et al., 1997; Eggleton et al., 1998; Mikkola and Tanikawa, 1998; Eggleton and Kiseleva-Eggleton, 2001; Valtonen and Karttunen, 2006; Fabrycky and Tremaine, 2007; Wu et al., 2007; Zdziarski et al., 2007; Takeda et al., 2008; Perets and Fabrycky, 2009). Neglecting the octupole-level approximation can cause changes in the dynamics varying from a few percent to completely different qualitative behavior.

Some other derivations of octupole–order equations of motion dealt with the secular dynamics in a general way, without using Hamiltonian perturbation theory or elimination of the nodes (Farago and Laskar, 2010; Laskar and Boué, 2010; Mardling, 2010; Katz and Dong, 2011). In these works there were no references to the discrepancy between these derivations and the previous studies. Also, note that the results of Holman et al. (1997) are based on a direct N-body integration, and thus are not subject to the errors mentioned above.

## 6 Conclusions

We have shown that the “standard” TPQ Kozai formalism (Kozai, 1962; Lidov, 1962) has been applied in inappropriate situations. A common error in the implementation of the relevant Hamiltonian mechanics (premature elimination of the nodes) leads to the (incorrect) conclusion that the conservation of the -component of each orbit’s angular momentum from the TPQ dynamics generalizes beyond the TPQ approximation. Correcting the formalism we find that the -components of both the inner and outer orbits’ angular momenta in general change with time at both the quadrupole and octupole level. The conservation of the inner orbit’s -component of the angular momentum (the famous ) only holds in the quadrupole-level test particle approximation. We have explained in details the source of the error in previous derivations (Appendix C).

We have re-derived the secular evolution equations for triple systems using Hamiltonian perturbation theory to the octupole-level of approximation (Section 2 and Appendix A, 4 and Appendix B). We have also shown that one can use the simplified Hamiltonian found in the literature (e.g., Ford et al., 2000a) as long as the equations of motion for the inclinations are calculated from the total angular momentum.

The correction shown here has important implications to the evolution of triple systems. We discussed a few interesting implications in Section 5. We showed that already at the quadrupole-level approximation the explicit assumption that the vertical angular momentum is constant can lead to erroneous results, see for example Figure 4. In this Figure we showed that far from the test particle limit in the quadrupole-level one can already find a significant difference in the evolutionary behavior. The correct results agree with the test particle limit only when (see Figure 5). We show in Appendix A.4 that at the quadrupole level of approximation, the inner eccentricity and the mutual inclination have a well defined maximum and minimum irrespective of the mass of the inner bodies. In the test particle limit these values converge to the well-known critical inclinations () for large oscillatory amplitudes.

The most notable outcome of the results presented here happens in the octupole-level of approximation (which we call the EKL formalism), when the inner orbit flips from prograde to retrograde with respect to the total angular momentum. Just before the flip the inner orbit has an excursion of extremely high eccentricity. In the presence of tidal forces (not included in this study) the outcome of a system can be different than the one assumed while using the TPQ formalism. Krymolowski and Mazeh (1999), Ford et al. (2000a), Blaes et al. (2002), Lee and Peale (2003) and Laskar and Boué (2010) present the correct octupole equations of motion. Had these authors integrated their equations for systems such as those presented in this paper, they could already have discovered the possibility of flipping the inner orbit.

## Acknowledgments

We thank Boaz Katz, Rosemary Mardling and Eugene Chiang for useful discussions. We thank Staffan Sderhjelm, our referee for very useful comments the improved the manuscript in great deal. We also thank Keren Sharon and Paul Kiel for comments on the manuscript. S.N. supported by NASA through a Einstein Postdoctoral Fellowship awarded by the Chandra X-ray Center, which is op- erated by the Smithsonian Astrophysical Observatory for NASA under contract PF2-130096. Y.L. acknowledges support from NSF grant AST-1109776. Simulations for this project were performed on the HPC cluster fugu funded by an NSF MRI award.

## Appendix A The Quadrupole level of Approximation

We develop the complete quadrupole-level secular approximation in this section. As mentioned, the main difference between the derivation shown here and those of previous studies lies in the “elimination of nodes” (e.g., Kozai, 1962; Jefferys and Moser, 1966), which relates to the transition the invariable plane (e.g., Murray and Dermott, 2000) coordinate system, where the total angular momentum lies along the -axis.

### a.1 Transformation to the Invariable Plane

We choose to work in a coordinate system where the total initial angular momentum of the system lies along the axis (see Figure 2),; the - plane in this coordinate system is known as the invariable plane (e.g., Murray and Dermott, 2000), and therefore we call this coordinate system the invariable coordinate system. We begin by expressing the vectors and each in a coordinate system where the periapse of the orbit is aligned with the x-axis and the orbit lies in the x-y plane, called the “orbital coordinate system,” and then rotating each vector to the invariable coordinate system. The rotation that takes the position vector in the orbital coordinate system to the position in the invariable coordinate system is given by (see Murray and Dermott, 2000, chapter 2.8, and Figure 2.14 for more details)

 r1,inv=Rz(h1)Rx(i1)Rz(g1)r1,orb , (27)

where the subscript “inv” and “orb” refer to the invariable and orbital coordinate systems, respectively. The rotation matrices and as a function of rotation angle, , are

 Rz(θ)=⎛⎜⎝cosθ−sinθ0sinθcosθ0001⎞⎟⎠ (28)

and

 Rx(θ)=⎛⎜⎝1000cosθ−sinθ0sinθcosθ⎞⎟⎠ . (29)

Thus, the angle between and is given by:

 cosΦ=^rT2,orbR−1z(g2)R−1x(i2)R−1z(h2)Rz(h1)Rx(i1)Rz(g1)^r1,orb, (30)

where are unit vectors that point along . In the orbital coordinate system, we have

 ^r1,2,orb=⎛⎜⎝cos(f1,2)sin(f1,2)0⎞⎟⎠, (31)

where () is the true anomaly for the inner (outer) orbit. Note that , so the Hamiltonian will depend on the difference in the longitudes of the ascending nodes; in a similar manner, the Hamiltonian depends on and only through expressions of the form and . Replacing in the Hamiltonian, eq. (2), we can now integrate over the the mean anomaly angles using the Kepler relations between the mean and true anomalies:

 dli=1√1−e2i(riai)2dfi , (32)

where for the outer orbit one should simply replace the subscript “1” with “2”.

### a.2 Transformation to Eliminate Mean Motions

Because we are interested in the long-term dynamics of the triple system, we now describe the transformation that eliminates the short-period terms in the Hamiltonian that depend of and . The technique we will use is known as the Von Zeipel transformation (for more details, see Brouwer, 1959).

Write the triple-system Hamiltonian in eq. (2) as

 H=HK1+HK2+H2, (33)

where and are the Kepler Hamiltonians that describe the inner and outer elliptical orbits in the triple system and describes the quadrupole interaction between the orbits. Note that is , and is the only term in that depends on or . We seek a canonical transformation that can eliminate the and terms from . Such a transformation must be close to the identity, since ; let the generating function be

 S(L∗j,G∗j,H∗j,lj,gj,hj) = 2∑j=1[L∗jlj+G∗jgj+H∗jhj] + α2S2(L∗j,G∗j,H∗j,lj,gj,hj) ,

where we indicate the new momenta with a superscript asterix, and is the non-identity piece of the transformation that we will use to eliminate . The relationship between the new and old canonical variables is

 pi=∂S∂qi=p∗i+α2∂S2∂qi (35)

and

 q∗i=∂S∂p∗i=qi+α2∂S2∂p∗i, (36)

where the momenta , and the coordinates . Because our generating function is time-independent, the new and old Hamiltonians agree when evaluated at the corresponding points in phase space:

 H(qi,pi)=H∗(q∗i,p∗i) (37)

when the phase space coordinates satisfy equations (35) and (36). Inserting these relations into the un-transformed Hamiltonian, and expanding to lowest order in , we have

 H(q∗i,p∗i)+α2∂H∂pi∂S2∂qi−α2∂H∂qi∂S2∂p∗i=H∗(q∗i,p∗i). (38)

Equating terms order-by-order in gives

 HK1(q∗i,p∗i)=H∗K1(q∗i,p∗i), (39)
 HK2(q∗i,p∗i)=H∗K2(q∗i,p∗i), (40)

and

 H2(q∗i,p∗i)+α22∑i=1∂H∂pi∂S2∂qi−α22∑i=1∂H∂qi∂S2∂p∗i=H∗2(q∗i,p∗i). (41)

Since the last two terms on the left-hand side of this latter equation are already , only the and parts of contribute. These Kepler Hamiltonians only depend on and , so there are only two non-zero partials of at order :

 H2(q∗i,p∗i)+α2∂HK1∂L1∂S2∂l1+α2∂HK2∂L2∂S2∂l2=H∗2(q∗i,p∗i). (42)

We must use the terms that depend on to cancel any terms in that depend on and . Note that is periodic in and with period (see equations (30) and (31)), so we can write

 H2(q∗i,p∗i)=α2h0+α2∞∑k1,k2=1hk1k2e−ik1l∗1−ik2