The Eccentric Kozai-Lidov Effect and Its Applications

The Eccentric Kozai-Lidov Effect and Its Applications

Smadar Naoz Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA;
snaoz@astro.ucla.edu
Abstract

The hierarchical triple body approximation has useful applications to a variety of systems from planetary and stellar scales to supermassive black holes. In this approximation, the energy of each orbit is separately conserved and therefore the two semi-major axes are constants. On timescales much larger than the orbital periods, the orbits exchange angular momentum which leads to eccentricity and orientation (i.e., inclination) oscillations. The orbits’ eccentricity can reach extreme values leading to a nearly radial motion, which can further evolve into short orbit periods and merging binaries. Furthermore, the orbits’ mutual inclination may change dramatically from pure prograde to pure retrograde leading to misalignment and a wide range of inclinations. This dynamical behavior is coined as the “eccentric Kozai-Lidov” mechanism. The behavior of such a system is exciting, rich and chaotic in nature. Furthermore, these dynamics are accessible from a large part of the triple body parameter space and can be applied to diverse range of astrophysical settings and used to gain insights to many puzzles.

1 Introduction

Triple systems are common in the Universe. They are found in many different astrophysical settings covering a large range of mass and physical scales, such as triple stars (e.g., Tokovinin, 1997; Eggleton et al., 2007; Tokovinin, 2014b, a), and accreting compact binaries with a companion (such as companions to X-ray binaries e.g., Grindlay et al., 1988; Prodan and Murray, 2012). In addition, it seems that supermassive black hole binaries and higher multiples are common and thus, any star in their vicinity forms a triple system (e.g., Valtonen, 1996; Di Matteo et al., 2005; Khan et al., 2012; Kulkarni and Loeb, 2012). Furthermore, considering the solar system, binaries composed of near earth objects, asteroids or dwarf planets (for which a substantial fraction seems to reside in a binary configuration, e.g., Polishook and Brosch, 2006; Nesvorný et al., 2011; Margot et al., 2015) naturally form a triple system with our Sun. Lastly, it was shown that Hot Jupiters are likely to have a far away companion, forming a triple system of star-Hot Jupiter binary with a distant perturber (e.g., Knutson et al., 2014; Ngo et al., 2015; Wang et al., 2015). Stability requirements yield that most of these systems will be hierarchical in scale, with a tight inner binary orbited by a tertiary on a wider orbit, forming the outer binary. Therefore, in most cases the dynamical behavior of these systems takes place on timescales much longer than the orbital periods.

The study of secular perturbations (i.e., long term phase average evolution over timescales longer than the orbital periods) in triple systems can be dated back to Lagrange, Laplace and Poincare. Many years later, the study of secular hierarchical triple system was addressed by Lidov (1961, where the English translation version was published only in 1962). He studied the orbital evolution of artificial satellites due to gravitational perturbations from an axisymmetric outer potential. Short time after that, Kozai (1962) studied the effects of Jupiter’s gravitational perturbations on an inclined asteroid in our own solar system. In these settings a relatively tight inner binary composed of a primary and a secondary (in these initial studies assumed to be a test particle), is orbited by a far away companion. We denote the inner (outer) orbit semi-major axis as (). In this setting the secular approximation can be utilized. This implies that the energy of each orbit is conserved separately (as well as the energy of the entire system), thus and are constants during the evolution. The dynamical behavior is a result of angular momentum exchange between the two orbits. Kozai (1962), for example, expanded the 3-body Hamiltonian in semi-major axis ratio (since the outer orbit is far away, is a small parameter). He then averaged over the orbits and lastly truncated the expansion to the lowest order, called the quadrupole, which is proportional to . Both Kozai (1962) and Lidov (1962) found that the inner test particle inclination and eccentricity oscillate on timescales much larger than its orbital period. In these studies the outer perturber was assumed to carry most of the angular momentum, and thus under the assumption of an axisymmetric outer potential the inner and outer orbits z-component of the angular momenta (along the total angular momentum) are conserved. This led to large variations between the eccentricity and inclination of the test particle orbit.

While the Kozai-Lidov111Although Lidov has published his work first, we are using here the alphabetical order for the name of the mechanism. mechanism seemed interesting it was largely ignored for many years. However, about 15-20 years ago, probably correlating with the detection of the eccentric planet 6 Cyg B, (Cochran et al., 1996), or the close to perpendicular stellar Algol system (Eggleton et al., 1998; Baron et al., 2012), the Kozai-Lidov mechanism received its deserved attention. However, while the mechanism seemed very promising in addressing these astrophysical phenomena, it was limited to a narrow parts of the parameter space (favoring close to perpendicular initial orientation between the two orbits, e.g., Marchal, 1990; Morbidelli, 2002; Valtonen and Karttunen, 2006; Fabrycky and Tremaine, 2007) and produced only moderate eccentricity excitations. Most of the studies that investigated different astrophysical applications of the Kozai-Lidov mechanism used the Kozai (1962) and Lidov (1962) test particle, axisymmetric outer orbit quadrupole-level approximation.

This approximation has an analytical solution which describes (for initially highly inclined orbits , see below) the large amplitude oscillations between the inner orbit’s eccentricity and inclination with respect to the outer orbit (e.g., Kinoshita and Nakai, 1999; Morbidelli, 2002). These oscillations have a well defined maximum and minimum eccentricity and inclination and limits the motion to either prograde () or retrograde () with respect to the outer orbit. The axisymmetric outer orbit quadrupole-level approximation is applicable for an ample amount of systems. For example, this approximation has appropriately described the motion of Earth’s artificial satellites under the influence of gravitational perturbations from the moon (e..g Lidov, 1962). Other astrophysical systems for which this approximation is applicable include (but are not limited to) the effects of the Sun’s gravitational perturbation on planetary satellites, since in this case indeed the satellite mass is negligible compared to the other masses in the system, and the planet’s orbit around the Sun is circular. Indeed it was shown that the axisymmetric outer orbit quadrupole-level approximation can successfully be used to study the inclination distribution of the Jovian irregular satellites (e.g., Carruba et al., 2002; Nesvorný et al., 2003) or in general the survival of planetary outer satellites (e.g., Kinoshita and Nakai, 1991), as well as the dynamical evolution of the orbit of a Kuiper Belt object satellite due to perturbation form the sun (e.g., Perets and Naoz, 2009; Naoz et al., 2010). Indeed this approximation is useful and can be applied in the limit of a circular outer orbit and a test particle inner object.

Recently, Naoz et al. (2011, 2013a) showed that relaxing either one of these assumptions leads to qualitative different dynamical evolution. Considering systems beyond the test particle approximation, or a circular orbit, requires the next level of approximation, called the octupole–level of approximation (e.g. Harrington, 1968, 1969; Ford et al., 2000b; Blaes et al., 2002). This level of approximation is proportional to . In the octupole–level of approximation, the inner orbit eccentricity can reach extremely high values, and does not have a well defined value, as the system is chaotic in general (Ford et al., 2000b; Naoz et al., 2013a; Li et al., 2014b, a; Teyssandier et al., 2013). In addition, the inner orbit inclination can flip its orientation from prograde, with respect to the total angular momentum, to retrograde (Naoz et al., 2011). We refer to this process as the eccentric Kozai–Lidov (EKL) mechanism. We note that here we follow the literature coined acronym “EKL” as oppose to the more chronologically accurate acronym “ELK.”

As will be discussed below the EKL mechanism taps into larger parts of the parameter space (i.e., beyond the range), and results in a richer and far more exciting dynamical evolution. As a consequence this mechanism is applicable to a wide range of systems that allow for eccentric orbits, or three massive bodies, from exoplanetary orbits over stellar interactions to black hole dynamics. The prospect of forming eccentric or short period planets through three body interactions was the source of many studies (e.g., Innanen et al., 1997; Wu and Murray, 2003; Fabrycky and Tremaine, 2007; Wu et al., 2007; Veras and Ford, 2010; Correia et al., 2011; Batygin et al., 2011; Naoz et al., 2011, 2012; Petrovich, 2015b, a). It also promoted many interesting application for stellar dynamics from stellar mergers (e.g., Perets and Fabrycky, 2009; Naoz and Fabrycky, 2014; Witzel et al., 2014; Stephan et al., 2016) to compact binary merger which may prompt supernova explosions for double white dwarf merger (e.g., Thompson, 2011; Katz and Dong, 2012), or gravitational wave emission for neutron star or black hole binary merger (e.g., Blaes et al., 2002; Seto, 2013).

2 The hierarchical three body secular approximation

In the three-body approximation, dynamical stability requires that either the system has circular, concentric, coplanar orbits or a hierarchical configuration, in which the inner binary is orbited by a third body on a much wider orbit, the outer binary (Figure 1). In this case the secular approximation (i.e., phase averaged, long term evolution) can be applied, where the interactions between two non-resonant orbits is equivalent to treating the two orbits as massive wires (e.g., Marchal, 1990). Here the line-density is inversely proportional to orbital velocity and 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 semi-major axes of the orbits. The gravitational potential is then expanded in semi-major axis ratio of , where () is the semi-major axis of the inner (outer) body (Kozai, 1962; Lidov, 1962). This ratio is a small parameter due to the hierarchical configuration.

Figure 1: Schematic description of the coordinate system and the angles used (not to scale). Left: The three bodies and the relative vectors. Here ’c.m.’ denotes the center of mass of the inner binary, containing objects of masses and . The separation vector points from to ; points from ‘c.m.’ to . The angle between the vectors and is . Right: Geometry of the angular momentum vectors and the definition of the relevant inclination angles. We show the total angular momentum vector (), the angular momentum vector of the inner orbit () with inclination with respect to and the angular momentum vector of the outer orbit () with inclination with respect to . The angle between and defines the mutual inclination . The invariable plane is perpendicular to , in other words the z axis is parallel to .

The hierarchical three body system consists of a tight binary ( and ) and a third body (). We define to 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 reduced into two separate Keplerian orbits: the first describing the relative tight orbit of bodies 1 and 2, and the second describes the wide 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 semi-major axes (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),

(0)

and in terms of the semi-major axises and :

(0)

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

(0)

The right term is often called the perturbing function as it describes the gravitational perturbations between the two orbits. The left two terms in Equation (2) are simply the energy of the inner and outer Kepler orbits. Note that the sign convention for this Hamiltonian is positive .

The frame of reference chosen throughout this review is the invariable plane for which the z axis is set along the total angular momentum, which is conserved during the secular evolution of the system (see figure 1), (e.g., Lidov and Ziglin, 1974). Another description used in the literature is the vectorial form (e.g. Katz et al., 2011; Boué and Fabrycky, 2014a), which has been proven to be useful to address different astrophysical setting. Considering the invariable plane it is convenient to adopt the canonical variables known as Delaunay’s elements, (e.g. Valtonen and Karttunen, 2006). These describe for each orbit three angles and three conjugate momenta.

The first set of angles are the mean anomalies, and (also often denote in the literature as and ), which describes the position of the object in their orbit. Their conjugate momenta are:

(0)

where subscripts denote the inner and outer orbits, respectively. The second set of angles are the arguments of periastron, and ( and ), which describes the position of the eccentricity vector (in the plane of the ellipse). Their conjugate momenta are the magnitude of the angular momenta vector of each orbit and (often used as and ):

(0)

where () is the inner (outer) orbit eccentricity. The last set of angles are the longitudes of ascending nodes, and ( and ). Their conjugate momenta are

(0)

(often denote as and ). Note that and are the magnitudes of the angular momentum vectors ( and ), and and are the -components of these vectors, (recall that the -axis is chosen to be along the total angular momentum ). In Figure 1 we show the configuration of the angular momentum vectors of the inner and outer orbit ( and , respectively) and and are the -components of these vectors, where the -axis is chosen to be along the total angular momentum . This conservation of the total angular momentum yields a simple relation between the component of the angular momenta and the total angular momentum magnitude:

(0)

The equations of motion are given by the canonical relations (for these equations we will use the notation):

(0)
(0)
(0)

where . Note that these canonical relations have the opposite sign relative to the usual relations (e.g., Goldstein, 1950) because of the sign convention typically chosen for this Hamiltonian.

As apparent from the Hamiltonian Eq. (2), if the semi-major axis ratio is indeed a small parameter then for the zeroth approximation each orbit can be described as a Keplerian orbit, for which its energy is conserved. Thus, we can average over the short timescale and focus on the long-term dynamics of the triple system. This process is known as the secular approximation, where the energy (semi-major axis) is conserved and the orbits exchange angular momentum. The short timescales terms in the Hamiltonian depend on and , and eliminating them needs is done via a canonical transformation. The technique used is known as the Von Zeipel transformation (Brouwer, 1959). In this canonical transformation, a time independent generating function is defined to be periodic in and , which allows the elimination of the short-period terms in the Hamiltonian and the details of this procedure are described in Naoz et al. (2013a) Appendix A2. Eliminating these angles from the Hamiltonian yields that their conjugate momenta and are conserved [see E1. (2)], thus yielding and , as expected. In the most general case of this three body secular approximation there are only two parameters which are conserved, i.e., the energy of the system (which also means that the energy of the inner and the outer orbits are conserved separately), and the total angular momentum .

The time evolution for the eccentricity and inclination of the system can be easily achieved from Equations (2)-(2)

(0)

and

(0)

where for the inner and outer orbit. See full set of equations of motions in Equations (8)-(8).

The lowest order of approximation, which is proportional to is called quadrupole level, and we find that an artifact of the averaging process results in conservation of the outer orbit angular momentum , in other words the system is symmetric for the rotation of the outer orbit. This was coined as the “happy coincidence” by Lidov and Ziglin (1976). Its significant consequence is that the this approximation should be only used for an axisymmetric outer potential such as circular outer orbits (Naoz et al., 2013a).

The next level of approximation, the octupole, is proportional to (see below) and thus the TPQ approximation can be successfully applied when this parameter is small for low inclinations (see below for numerical studies). However, close to perpendicular systems are extremely sensitive to this parameter.

A popular procedure which was done in earlier studies (e.g., Kozai, 1962) used the “elimination of nodes” (e.g., Jefferys and Moser, 1966). This describes the a simplification of the Hamiltonian by setting

(0)

This relation holds in the invariable plane when the total angular momentum is conserved, such as in our case. Some studies that exploited explicitly this relation in the Hamiltonian incorrectly concluded [using Equation (2)] that the -components of the orbital angular momenta are always constant. As showed in Naoz et al. (2011, 2013a), this leads to qualitatively different evolution for the triple body system. We can still use the Hamiltonian with the nodes eliminated as long as the equations of motions for the inclinations are derived from the total angular momentum conservation, instead of using the canonical relations (Naoz et al., 2013a).

2.1 Physical picture

Considering the quadrupole level of approximation (which is valid for axisymmetric outer orbit potential) for an inner test particle (either or goes to zero) the conserved quantities are the energy and the z-component of the angular momentum. In other words the Hamiltonian does not depend on longitude of acceding nodes () and thus the z-component of the inner orbit angular momentum, , is conserved and the system is integrable. In this case the equal precession rate of the inner orbit’s longitude of ascending nodes () and the longitude of the periapsis () means that an eccentric inner orbit feels an accumulating effect on the orbit. The the resonant angle , will librate around or which cause large amplitude eccentricity oscillations of the inner orbit.

In that case (circular outer orbit, in the test particle approximation) the conservation of the z component of the angular momentum  yields oscillations between the eccentricity and inclination. The inner orbit will be more eccentric for smaller inclinations and less eccentric for larger inclinations.

2.2 Circular outer body

In this case the gravitational potential set by the outer orbit is axisymmetric, and thus the quadrupole level of approximation describes the behavior of the hierarchical system well. We will consider two possibility, in the first one of the members of the inner orbit is a test particle, (i.e., either or are zero). In the second we will allow for all three masses to be non-negligible.

2.2.1 Axisymmetric Potential and Inner Test Particle - TPQ

Following Lithwick and Naoz (2011) we call this case the test particle approximation quadrupole (TPQ). Without loss of generality we take , the Hamiltonian of this system is very simple and can be written as:

(0)

where

(0)

where (e.g., Yokoyama et al., 2003; Lithwick and Naoz, 2011)222Note that unlike the Hamiltonian that will be presented in the next section [Equation (2.2.2)] this Hamiltonian only describes the test particle approximation. .

Figure 2: Cross section trajectory of the TPQ in the (top panels) and (bottom panels) planes. We define . The dashed horizontal lines in the top panels shows the critical inclination for which . The separatrix is associated with for and for , as depicted in the Figure. Left panels shows the case for and and (librating) and and (circulating). Right panels shows the case for and (librating) and and (circulating). Figure adopted from Lithwick and Naoz (2011).

At this physical setting the octupole level of approximation is zero and the inner orbit’s angular momentum along the z axis is conserved (), where is the specific z component of the angular momentum. Since both and are conserved, a new constant of motion can be defined. It is convenient (for reasons that will be identified in Section 2.3.1) to define the following constant (Katz et al., 2011):

(0)

which is a simple function of the initial conditions. The system is integrable and has a well defined maximum and minimum eccentricity and inclination. To find the extreme points we set in the time evolution equation [see Equation (8), quadrupole part] and find that the values of the argument of periapsis that satisfy this condition are , where . Thus, the resonant angle has two classes of trajectories, librating and circulating. On circulating trajectories, at , the eccentricity is smallest and the inclination is largest, and visa versa for . In Figure 2, librating trajectories (or libration modes) are associated with bound oscillations of and circulating trajectories (or circulation modes) are not constrained to a specific regime. The separatrix is the trajectory which separates the two modes of behavior, as we elaborate below.

The conservation of implies:

(0)

where and are the initial values. Note that in this case (TPQ) . Since the energy is also conserved, plugging in for the circulating mode we find

(0)

and for in equation (2.2.1) we find:

(0)

where represents the initial conditions plugged in equation (2.2.1). From equations (2.2.1) and (2.2.1) one can easily find the minimum eccentricity and maximum inclination, and from equations (2.2.1) and (2.2.1) the maximum eccentricity and the minimum inclination. A special and useful case is found by setting initially and , for this case the maximum eccentricity is

(0)

Solving the equations for instead we can find

(0)

which gives and , known as Kozai angles. These angles represent the regime where large eccentricity and inclination oscillations are expected to take place. The value marks the seperatrix depicted in Figure 2.

2.2.2 Axisymmetric Potential Beyond the Test Particle Approximation

In this case we still keep the outer orbit circular, thus the quadrupole level of approximation still valid, but we will relax the test particle approximation. The quadrupole level hamiltonian can be written as:

(0)

where

(0)

Note that in this form of Hamiltonian the nodes ( and ) have been eliminated, allowing for a cleaner format, however this does not mean that the z-component of the inner and outer angular momenta are constant of motion (as explained in Naoz et al., 2011, 2013a).

Relaxing the test particle approximation (i.e., none of the masses have insignificant mass) already allows for deviations from the nominal TPQ behavior. This is because now is no longer conserved and instead the total angular momentum is. Note that the outer potential is axisymmetric and . The system is still integrable and has well define maxima and minima for the eccentricity and inclination. The conservation of the total angular momentum, i.e., sets the relation between the maximum/minimum total inclination and inner orbit eccentricity.

(0)

Note that in the quadrupole-level approximation , and thus , are constant. The right hand side of the above equation is set by the initial conditions. In addition, , and [see eqs. (2) and (2)] are also set by the initial conditions. Using the conservation of energy we can write, for the minimum eccentricity/maximum inclination case (i.e., setting )

(0)

The left hand side of this equation, and the remainder of the parameters in Equation (2.2.2) are determined by the initial conditions. Thus solving equation (2.2.2) together with (2.2.2) gives the minimum eccentricity/maximum inclination during the system evolution as a function of the initial conditions. We find a similar equation if we set for the maximum eccentricity/minimum inclination :

(0)

Equations (2.2.2) and (2.2.2) give a simple relation between the total minimum inclination and the maximum inner eccentricity as a function of the initial conditions.

Figure 3: Comparison between the TPQ formalism (dashed blue lines) and the full quadrupole calculation (solid red lines). The system has an inner binary with and , and the outer body has mass . The orbit separations are AU and AU. The system was set initially with and , , and relative inclination . The panels show from top to bottom, the mutual inclination , and , which in the TPQ formalism is constant (dashed line). Figure adopted from Naoz et al. (2013a).

An interesting consequence of this physical picture is if the inner binary members are more massive than the third object. We adopt this example from Naoz et al. (2013a) and consider the triple system PSR B162026. The inner binary contains a millisecond radio pulsar of and a companion of (e.g., McKenna and Lyne, 1988). We adopt parameters for the outer perturber of (Ford et al., 2000a) and set (see the caption for a full description of the initial conditions). Note that Ford et al. (2000a) found , which means that the quadrupole level of approximation is insufficient to represent the behavior of the system. We choose, however, to set to emphasis the point that even an axisymmetric outer potential may result in a qualitative different behavior if the TPQ approximation is assumed. For the same reason we also adopt a higher initial value for the inner orbit eccentricity ( compared to the measured one ). The time evolution of the system is shown in Figure 3. In This Figure we compare the z-component of the angular momentum (solid red line) with (dashed blue line), which is the angular momentum that would be inferred if the outer orbit were instantaneously in the invariable plane, as found in the TPQ formalism.

Taking the outer body to be much smaller than the inner binary (i.e., ), as done in Figure 3, yields yet another interesting consequence for relaxing the test particle approximation. In some cases large eccentricity excitations can take place for inclinations that largely deviate from the nominal range of the Kozai angles of . The limiting mutual inclination that can result in large eccentricity excitations can be easily found when solving Equations (2.2.2) and (2.2.2), since they depend on mutual inclination, as noted by Martin and Triaud (2015b). This evolution is shown in Figure 4, where large eccentricity oscillation for the inner binary is achieved for an initial mutual inclination of . This behavior, as expected from the Equations, is sensitive to the eccentricity of the outer orbit.

Figure 4: Small mass outer perturber that induces large eccentricity excitation away from the nominal range of the Kozai angles of . We consider  M,  M,  M,  AU and  AU. Both outer and inner eccentricities are set initially to zero, and also set initially and . We show two examples, the first shows the eccentricity excitations for as expected initial mutual inclination of , where in this case and . This produces eccentricity excitation with . We also consider an example for which the mutual inclination is set initially to be . In this case and . The latter parameters are adopted from Martin and Triaud (2015b), which leads to maximum inner eccentricity of . Note that in both examples is close to the nominal Kozai angles range.

In the circular outer orbit case, the regular oscillations of the eccentricity and inclination yields a well defined associated timescale. This can be easily achieved by considering the equation of motion of the argument of periapsis [see the part that is proportional to in Equation (8)]. More precisely, , where is given in Eq. (2.2.2). Integrating between the well defined maximum and minimum eccentricity Antognini (2015) found a numerical factor , and got

(0)

This timescale is in a good agreement with the numerical evolution.

2.3 Eccentric outer orbit

2.3.1 Inner Orbit’s Test Particle Approximation

In this approximation we will allow for an eccentric outer orbit but will restrict ourselves to take the mass of one of the inner members to zero, which yields . In the test particle limit, the outer orbit is stationary and the system reduces to two degrees of freedom. The eccentric outer orbit yields the quadrupole level of approximation inadequate and thus we consider the test particle octupole (TPO) level here. This approximation is extremely useful in gaining an overall understanding of the general hierarchical system and the EKL mechanism. The hamiltonian of this system is very simple and can be written as (e.g., Lithwick and Naoz, 2011),

(0)

where

(0)

is defined in Equation (2.2.1), and we reiterate it here for completeness,

and

(0)
Figure 5: Time evolution example of the TPO (test partial octupole) approximation (red lines) and the TPQ (test particle quadrupole) approximation (blue lines). Left panels show high inclination flip while the right panels show the low inclination flip (see text). In this example we consider the time evolution of a test particle at  AU around a  M intermediate black hole located  pc from the massive black hole in the center of our galaxy ( M). In the left panels the system initially is set with , , , and . In the right panels the system is initially set with , , , and . In the top panels we show the inclination and in the bottom the inner orbit eccentricity as .

In this case the z component of the outer orbit is not conserved and the system can flip from to (Naoz et al., 2011, 2013a). The flip is associated with an extremely high eccentricity transition (see for example Figure 5). The octupole level of approximation introduces higher order resonances which overall render the system to be qualitatively different from a system at which the quadrupole level of approximation is applicable. We will begin by reviewing the different effects in the systems which can be divided into two main initial inclination regimes.

High initial inclination regime and chaos

When the system begins with in a high inclination regime the resonance arises from the quadrupole level of approximation can cause large inclination and eccentricity amplitude modulations. Recall that this angle range is associated with the TPQ seperatrix. The octupole-level of approximation is associated with high order resonances that result in extremely large eccentricity peaks, flips (see Figure 5) as well as chaotic behavior (as explained below). As can be seen from equation (2.3.1) these resonances arise from higher order harmonics of the octupole-level Hamiltonian: and . A useful tool to analyze this system is in the form of surface of section (see fore example Figure 6). For a two-degrees of freedom system, the surface of section projects a four-dimensional trajectory on a two-dimensional surface. The resonant regions are associated with fixed points and chaotic zones are a result of the overlap of the resonances between the quadrupole and the octupole resonances (Chirikov, 1979; Murray and Holman, 1997).

Figure 6: Surface of section for and . This initial configuration is associated with high initial inclination . The quadrupole level resonances can clearly be seen (the big islands) as well as the emergence of high order resonances (the small islands). Figure adopted from Li et al. (2014a).

Figure 6 shows the surface of section for and , which is associated with high initial inclination . In this Figure we can identify three distinct regions: Òresonant regions,Ó Òcirculation regions,Ó and Òchaotic regionsÓ. The resonant regions are associated with trajectories of which the momenta ( and ) and the angles and ) undergo bound oscillations. The system is classified in a liberation mode and the trajectories are quasi-periodic. The libration zones in the TPQ approximation are shown in Figure 2, and for the TPO in Figure 6. The circulation regions describes trajectories for which the coordinates are not constrained to a specific interval, and can take any value. Note that both resonant and circulatory trajectories map onto a one-dimensional manifold on the surface of section. On the contrary, chaotic trajectories map onto a two-dimensional manifold. In other words, while quasi-periodic trajectories form lines on the surface of section, chaotic trajectories are area-filling regimes. Embedded in the chaotic region, the small islands correspond to the higher, octupole order resonances, which are also quasi-periodic . The flip from to covers large parts of the parameter space as can be seen in Figure 7 right panel.

Figure 7: High inclination flip parameter space Right panel shows the results of numerical integrated systems associated maximum eccentricity (color coded as ) in the parameter space for , after . Systems above the black line flipped. Left panel shows the comparison with the analytical conditions derived by Katz et al. (2011). Open circles are the result of a numerical integration (red indicates systems that flipped and blue is for those that did not). Here the solid lines represent the flip conditions which for and is reduced to equation (2.3). The bottom left panel shows the case of , and note that it shows only part of the parameter space. Left panels are adopted from Katz et al. (2011) and right panel is adopted from Li et al. (2014b).

In some cases an analytical condition for the flip can be achieved by averaging over a quadrupole cycle (Katz et al., 2011). This averaging process yields a constant of motion

(0)

where the function is defined by:

(0)

where and are the complete elliptic functions of the first and second kind, respectively. For initial high inclination a flipping critical value for the octupole pre-factor is a function of the initial inclination and the approximations takes a simple form

(0)

where , was defined in Equation (2.2.1) and the subscript “0” marks the initial conditions. We note that in this TPO case is no longer constant (unlike the TPQ case). The parameter has the range . For cases where , i.e., and Equation (2.3) takes a simple form:

(0)

This approximation is valid for . The validity of this approximation for different initial values of and are shown in the left panels in Figure 7.

Figure 8: Flip timescales. We consider the following supermassive black hole binary system  M,  M (note that in this case ). The other parameters of this system are:  pc  pc and . The system is sent initially with , and for the left panels and for the right panels. The initial inclinations considered are colored labeled in the Figure. Note the difference in flip timescale as a function of initial inclinations. In the inset we show the inner orbit eccentricity as a function of time for the TPQ (brown line) and TPO (green) for the initial setting of and case, which emphasis the different short (quadruple) timescales between the TPQ and TPO level of approximation.

A timescale for the high inclination oscillation or flip is difficult to quantify since the evolution is chaotic. Furthermore, numerically it seems that the timescale for the first flip depends on the inclination (as can be seen in Figure 8). However, an approximate analytical condition, for the regular (none chaotic) mode was achieved recently by Antognini (2015), following Katz et al. (2011) formalism. This timescale has the following functional form:

(0)

where

(0)

and note that defined in Antognini (2015) is simply in the notation used here. The upper limit of the integral in Equation (2.3) is easy to find, since for the z component of the angular momentum is zero, thus

(0)

and the minimum limit of the integral is found from solving . This timescale takes a simple form, for setting initially and :

(0)

In the TPO level of approximation the short (quadrupole) timescales differ from the associated timescale at the TPQ level. In other words following the evolution of the same system, once by using the TPO and once using the TPQ yields different timescales, as depicted in the inset of Figure 8. This is because the Hamiltonian (i.e., the energy) is slightly different as the TPO includes the octupole term. Thus, the two calculations sample somewhat different values of the system energy. The difference is within a factor of a few as it represents the range of the phase space away from the seperatrix (See Figure 2 for the different oscillation’s amplitudes for given initial different energy.

Low initial inclination regime

Figure 9: Surface of section for and for , this associated with low initial inclination , Figure adopted from Li et al. (2014a). See similar plots in Petrovich (2015a), reproducing this analysis.

The octupole level of approximation yields an interesting behavior even beyond the Kozai angles. This is a result of the octupole level harmonics, i.e., and . Since the low order resonances are missing, the co-planer flip is not associated with chaotic behavior. Figure 9 shows the surface of section for two low inclination examples, specifically and for .

As can seen from Figure 5 (as well as Figures 6 and 9) the two inclination regimes exhibit qualitative differences. The high inclination flip is driven by the quadrupole level resonance with the actual flip arrises by accumulating effects from the high order resonates. Furthermore, this flip, many times, is associated with a chaotic behavior (Lithwick and Naoz, 2011; Li et al., 2014a). On the other hand, the low inclination flip is due to a regular trajectory. In addition, this flip takes place on a much shorter timescale than the high inclination flip.

Figure 10: Low inclination flip criterion: Comparison between the analytical expression Equation (2.3), solid line, and numerical integration, (green crosses mark no flip after , and blue crosses systems that flipped). The system’s parameters are:  M, ,  M,  AU,  AU. The outer orbit eccentricity was changed to match the values indicated on the vertical axis. The system was initially set with , , and as indicated in the figure. Figure adopted from Li et al. (2014b).

Similarly to the analytical approximation for the high inclination flip conditions, Li et al. (2014b) achieved an analytical condition for the low inclination flip, after averaging over the flip timescale

(0)

Comparing this condition to the high inclination condition Equation (2.3), also emphasis the qualitative difference between these two regimes.

The low inclination regime yields a flip timescale that can be easily found by setting . Li et al. (2014b) found a expression for the flip timescale:

(0)

where is the initial inner orbit eccentricity and is the energy that corresponds to and the rest of the initial conditions. The reason for the two integrals is because if initially , the inner orbit eccentricity, , decreases before it increases, otherwise if .

Figure 11: Flip and maximum eccentricity parameter space in two hierarchical planets configuration. The color describes the maximum eccentricity reached over integration time of (right top panel) and the flip ratio, defined as the time the total inclination spends over from the entire integration time (the other two panels). The top two panels show the phase space corresponding to (right) and the flip ratio (left) as a function of the initial outer orbit eccentricity () and the initial mutual inclination. Note that both follow exhibit interesting behavior at similar parts in the parameter space. However, for initial large inclination , the flip is suppressed. The system considered here has the following parameters:  M,  M,  M,  AU, and  AU. The bottom panels shows the flip ratio in the initial phase space. The system considered in this panel has the same parameters as the top two panels, but with and varying . The flip condition for the TPQ, following the condition in Equation (2.3) is shown in purple dots. The TPQ analysis for the top left (bottom) panel suggests that all systems above (below) the “TP” dotted line are expected to flip,. The solid line represents the stability condition, see Equation (3). Figure adopted from Teyssandier et al. (2013).

2.3.2 Beyond the Test Particle Approximation

Relaxing the test particle approximation leads to some qualitative differences. The first is that now one of the inner bodies can torque the outer body, and thus suppress the flip. This also causes a shift in the parameter space of the flip condition and the extreme eccentricity achieved compared to the TPQ case (see Figure 11). While the value of the maximum of is similar to that in the TPQ case, large eccentricity excitations may take place in different parts of the parameter space (compare Figure 11 to Figure 7). In particular, in the high inclination regime, the flips and the large eccentricity excitations of the TPQ case are concentrated around but in the full case they can shift to lower mutual inclinations and tap to larger range of inclinations (Figure 11). This is mainly because the outer orbit is being torqued by the inner orbit. Teyssandier et al. (2013) studied the effect of a similar mass companion and showed that if the outer body mass is reduced to below twice the smallest mass of the inner orbit, the flip and large eccentricity excitations are suppressed for large parts of the parameter space.

The system’s hamiltonian is (here again the nodes were eliminated for simplicity, but the z-component of the angular momenta are not conserved):

(0)

where is define in equation (2.2.2) and we copy it here for completeness

the octupole level approximation is:

where

(0)

and