Roxana Feier ‘12

Harvard University, Department of Astronomy

This paper aims to investigate the Hill stability of hierarchical three-body systems by taking into consideration the secular Kozai effect. We start by deriving a closed-form Hill stability criterion for circular, inclined orbits by expanding on previous analytical work. This establishes the first major result of our investigation. We then ask how the condition for stability changes when secular evolution is considered. Simulation results are shown to be largely in agreement with our conclusions. More general criteria, valid for eccentric orbits, are considered next. In the quadrupole approximation and for the test-particle limit, these stability conditions depend on the binary inclination and eccentricity through a quantity that is conserved during Kozai oscillations. We discuss the implications of this fact for our results in the circular case, and suggest that going beyond the quadrupole approximation is essential. In the second part of this paper, we go beyond the gravitational interactions studied in the classical three-body problem, and consider the role played by tidal friction in shaping the evolution of a three-body system. In highly eccentric orbits, tidal effects become noticeable during close encounters at the periastron. Using orbital data for 22 Kuiper belt objects, we investigate how Kozai cycles can drive these tidal interactions. Current orbital conditions, particularly as they relate eccentricity to orbital period, indicate the effectiveness of this combined Kozai-tidal friction effect. Furthermore, we show that currently circular orbits have a shorter tidal circularization timescale. This suggests that such systems originated in the high inclination/high eccentricity regime, but Kozai oscillations followed by tidal friction reduced these orbital parameters. In this proposed Kozai-tidal friction evolution, we suggest a final distribution of eccentricities that is expected after enough time has passed to allow the complete evolution to occur.


The three-body problem represents a very old topic in celestial mechanics, studied since gravity was first understood. This problem attempts to solve for the motion of three astronomical objects as they interact gravitationally with one another, given some specified initial conditions. Although this appears to be a simple problem, the general equations of motion cannot be solved analytically (Valtonen & Karttunen 2006; Georgakarakos 2008). Nonetheless, certain configurations of the three bodies are better understood than others. Of importance in our analysis are hierarchical triple systems, in which two objects that are close together in a binary orbit a more distant and massive perturbing object. Notable examples of this kind are Sun-orbiting binary asteroids and binary trans-Neptunian objects, the latter of which will be discussed extensively here.

One issue of interest for a three-body system is its stability, especially as this is easier to solve for than the full equations of motion. A system is considered to be stable if no collision or escape of one of the bodies occurs after a large number of orbital periods. In this paper, we would like to investigate what is called the Hill stability of hierarchical triple systems. For the inner binary, we say that the secondary component is within the Hill sphere of the primary if the acceleration induced by the latter dominates that of the more distant companion (Valtonen & Karttunen 2006). If this condition does not apply, then the secondary will escape from the system.

Stability is typically studied by considering the gravitational interactions between the three bodies at each given point in time, often with simplifying assumptions such as one of the objects being massless. In this paper, we start by deriving a closed-form Hill stability condition for circular, inclined orbits. In addition, we would like to investigate secular stability of three-body systems. Here, secular means that gravitational effects are averaged over the entire periods of the orbiting objects. An important secular effect is a periodic change in the inclination and eccentricity of the smaller binary, a process known as Kozai oscillations. In previous studies of Hill stability, both analytical and numerical, the impact of Kozai cycles is not fully addressed.

Our strategy is to consider previously established stability criteria and ask how these change when the Kozai effect is accounted for. After an overview of Kozai cycles, we first derive a Hill stability condition by expanding on work done by Innanen (1980) on circular, inclined three-body systems. Within the same section, we modify the criterion we established to account for Kozai oscillations. We show that our conclusions are largely confirmed by simulation results. We then ask whether a similar approach can be used for more general stability criteria. Specifically, we work with the results in Li et al. (2010), which give a closed-form expression for the critical separation for arbitrary eccentricities of the inner and outer orbits. We conclude with a discussion comparing what we deduced from the two papers and explaining why we think the results in Li et al. are problematic.

Following this theoretical analysis, we consider the way in which Kozai cycles affect the evolution of Kuiper belt binaries with known orbital parameters. This is an important discussion, since our current understanding of what Kuiper belt objects are—planetesimals disrupted into inclined, eccentric orbits by the outward migration of Neptune—does not fully explain the distribution of orbital parameters we see for Kuiper belt systems today. In this work, we go beyond considering just gravitational forces, as the classical three-body problem does, and take into account tidal interactions. In particular, we discuss how large eccentricity variations induced by Kozai oscillations lead to intensified tidal interactions within the binary.

Kuiper belt systems are next examined. We start by describing the data, and continue by assessing the Hill stability of these systems, based on the theoretical discussion above. Tidal effects are discussed in the next three sections, with a higher emphasis on Kozai oscillations. Lastly, we summarize our results and suggest further lines of research.

Secular stability criteria

The basic strategy for deriving secular stability criteria is taking known conditions for stability, putting an arbitrary system through a full Kozai evolution, and assessing whether all the intermediary orbital parameters satisfy the same stability condition. In that case, the system will be stable with Kozai oscillations accounted for. Before we go on and explore this idea further, a short overview of the Kozai effect is necessary.

Kozai cycles

Qualitatively, the most important aspect of Kozai oscillations is that they enable a long period oscillation between high inclination and high eccentricity, which can alter significantly the evolution of a three-body system. There are many interesting examples of astronomical bodies for which Kozai cycles are thought to play an important role. One such class of objects are Kreutz sungrazer comets. As the name suggests, these comets approach the Sun at perihelion within a very small distance, of about one solar radius (Heggie & Hut 2003). One hypothesis for their existence is that they started out on high-inclination, low-eccentricity orbits, but due to Kozai cycles induced by perturbations of the Solar System planets, they evolved into what we see today.

The Kozai effect typically arises in hierarchical triple systems, where the massive, distant object acts as the perturber on the motion of the inner binary. These perturbations affect several parameters of the inner orbit: the inner eccentricity ei, the inclination i taken relative to the outer orbit, the argument of pericenter ω, and the longitude of the ascending node Ω. We will primarily focus on the first three of these orbital elements—eccentricity, inclination, and argument of pericenter, with the latter defined as the angle between the orbit’s periapsis and its ascending node, the point where the secondary binary component crosses the plane of reference from South to North.

Under certain approximations that will be detailed later, two important quantities are preserved during Kozai oscillations, namely K=\sqrt{1-e_i^2}cos(i) and L=e_i^2-\sin^2{i}(5e^2_i\sin^2{\omega} + 1 - e_i^2) (Heggie & Hut 2003). The first of the two constants is proportional to the component of the inner angular momentum perpendicular to the outer orbit, and the second is the perturbation potential of the third object on the binary, averaged over the two orbits. Because these parameters remain constant, Kozai oscillations lie along the curves of intersection of the two surfaces described by the equations above. It turns out that two distinct configurations are possible. The first one, typical of small initial inclinations, corresponds to when the argument of pericenter ω circulates, meaning that it takes values in the entire range 0°-360°. For larger inclinations, ω librates around 90° or 270°. The latter situation is the standard definition of Kozai cycles, and occurs for i between about 39.2° and 140.8°.  [1]

The two constants K and L are determined by the initial configuration of the three-body system, and they also induce the limits within which the inclination, eccentricity, and argument of periapsis take values. In particular, the masses of the three objects and their mutual separations do not affect how much i, ei, or ω may change, but they affect the timescale on which this periodic variation occurs.

As hinted above, K and L are only constant under certain assumptions. Until recently, it was believed that a quadrupole approximation—a standard reduction in Hamiltonian mechanics which very roughly means that the eccentricity of the outer orbit does not change—was a sufficient condition for K and L to remain constant. However, new results by Naoz et al. (2011) show this not to be the case, unless the secondary object within the binary has negligible mass (the so-called “test particle” limit).

Beyond these simplifying assumptions, Kozai cycles become a more subtle matter, as K and L both incur some variability. We will discuss how this becomes relevant when considering the Hill stability of systems which undergo Kozai oscillations. However, all analytical derivations will assume that K (and occasionally L are constant), and so later results should be treated not as precise formulas, but rather as qualitative descriptions of stability.

In what follows, we make use of two previous works to adapt the stability criteria they establish to account for the Kozai effect. We work exclusively with hierarchical triple systems, in which a perturber (sometimes referred to as the “star”) is orbited around by a binary system consisting of a primary (“planet”) and a secondary (“satellite”).

The Hill stability of circular, inclined systems

The question of stability was initially addressed for planar, circular three-body systems. The first analytical work that investigates Hill stability for arbitrary inclinations is a paper by Innanen (1980). It thus seems natural that we start with this work in investigating how stability criteria are affected by Kozai cycles. We would like to establish the limiting orbital radius for the satellite to be retained by the planet, against the perturbing influence of the star. By setting the total relative acceleration of the satellite to 0, Innanen derives a formula which relates the critical radii for direct (rd) and retrograde (rr) motion corresponding to the same mutual inclination i, where 0° ≤ i ≤ 90°:

\frac{r_r}{r_d} = (1 + \frac{f}{2} + (f + (\frac{f}{2})^2)^{1/2})^{2/3} (1)

Here, f=4(cos2i)/3. For ease of notation, denote z = f/2 + [f + (f/2)2]1⁄2, which gives

\frac{r_r}{r_d} = (1+z)^{2/3} (2)

A priori, it is not obvious why rd and rr should be different for the same value of the inclination. The reason for this asymmetry is that the Coriolis acceleration points in different directions in the two cases, in a way which makes retrograde orbits more stable than prograde ones (Hamilton & Burns 1991). The goal now is to determine rd and rr separately as functions of the inclination i. For this, there needs to be another equation relating them. Again from Innanen (1980), we have the following relationship between the prograde and retrograde angular rates ωd and ωr:

\omega_d - \omega_r = 2\Omega \cos{i} (3)

where Ω is the angular velocity of the planet around the star. From the equations of uniform circular motion, we know that

\omega_d = \sqrt{\frac{Gm_{planet}}{r_d^3}} (4)

\omega_r = \sqrt{\frac{Gm_{planet}}{r_r^3}} (5)


\Omega = \sqrt{\frac{Gm_{star}}{r_{out}^3}} (6)

Here rout is the separation between the binary and the star. Using this along with equation (2), it follows that

\sqrt{\frac{Gm_{planet}}{r_d^3}} (1-\frac{1}{1+z}) = 2 \sqrt{\frac{Gm_{star}}{r_{out}^3}}\cos{i} (7)


r_d = r_{out} (\frac{m_{planet}}{4m_{star}})^{1/3}(\frac{z}{1+z})^{2/3}\frac{1}{\cos{i^{2/3}}} (8)


r_r = r_{out} (\frac{m_{planet}}{4m_{star}})^{1/3}(\frac{z}{\cos{i}})^{2/3} (9)

Note that with the notation used by Innanen, the inclination varies between 0 and 90 degrees, but for each i in this range we get two radii, depending on whether the motion of the satellite is prograde or retrograde relative to the primary orbit. We can set up a different framework, where i varies between 0 and 180 degrees, with 0° ≤ i ≤ 90° corresponding to prograde motion and 90° ≤ i ≤ 180° corresponding to retrograde orbits. Note that with this “change of coordinates” the equation for rr stays the same, since the radius depends on inclination only through the quantity cos2i, which is the same as cos2(180°−i).

For a given system, the quantity rout·(mplanet/4mstar)1/3 is constant, and remains so even as the satellite undergoes Kozai oscillations. Because here we are not concerned with a particular system, and so the value of this constant is not important, we will just scale rd and rr so that rd is 1 at i = 0°. Figure 1 plots the limiting radius as a function of inclination, across the whole domain from 0 to 180 degrees (continuous line).

Figure 1. Stability criterion for circular, inclined three body systems with Kozai (dashed line) and without (continuous line). The stability region lies below each of the two curves. When accounting for secular effects, a small (essentially 0) initial eccentricity is assumed. The area between the curves consists of quasistable systems, which only become unstable after an amount of time comparable to the Kozai timescale has passed.

This establishes the first major result of our work. We have derived a closed-form Hill stability condition valid for hierarchical triple systems of arbitrary inclinations. Previous closed-form criteria have been limited to direct and retrograde planar orbits. More general work which accounts for nonzero eccentricities has been done, for example, by Donnison (2010), but the results are not in a closed form. Similarly, general work by Li et al. (2010) provides a closed-form stability condition, but the results are problematic in ways that will be addressed later on. In addition, the criterion in Li et al. is weaker in the circular case than what we have derived. This is discussed more extensively later.

Next, we would like to address how the stability condition derived above is changed by Kozai oscillations. A simple analysis shows that systems which appear stable over a small number of orbits may be unstable over long timescales, when secular effects are accounted for. Outside the regime where Kozai effects are strong, meaning for inclinations between roughly 0°-40° and 140°-180°, the result obtained above remains largely unchanged. Now, consider a three-body system with initial inclination i0, and whose outer and inner orbits start out circular. This inclination can be exchanged for eccentricity through Kozai oscillations. In particular, for an initially circular orbit, and in the test particle limit where the satellite is massless, [2] the maximum eccentricity achieved through Kozai is e_{max} = \sqrt{1-\frac{5}{3}\cos^2{i}} (Koza 1962). For the test particle limit in the quadrupole approximation, the minimum inclination corresponding to the maximum eccentricity is given by the constraint that K = \sqrt{1-e^2}\cos{i} is preserved, so

i_{min} = \arccos{(\frac{\cos{i_0}}{\sqrt{1-e^2_{max}}})}

At this point of maximum eccentricity and minimum inclination, we would like to establish whether the system still lies within the non-secular stable regime of Figure 1. Because the orbit is now eccentric, there is the question of how to adapt this to the circular orbit considered by our stability condition. One reasonable choice is to regard this now eccentric orbit as a circular one with the radius equal to the apocenter distance. If this separation still renders the system as stable, then we know Hill stability to hold everywhere else in the orbit, where the physical distance is smaller.

This analysis can be done for a system with specified initial conditions, but there is also a more systematic approach. Again, consider a system with initial inclination i0 and initial orbital radius r0, but this time take i0 and r0to be variables rather than fixed values. When this system achieves its maximum Kozai eccentricity, its coordinates in the inclination-separation domain from Figure 1 are (imin, r0 (1+emax)). If the goal is to plot the initial conditions at which a system is borderline stable with Kozai accounted for, we want to ask for the most unstable point to lie on the graph shown in Figure 1. Therefore, for each i, we want to solve for r0 such that r0 (1+emax)=rcritical.

The results of this approach are shown in Figure 1. In this plot, everything below the lower curve is Hill stable even with Kozai changing the eccentricity and inclination of the system. The area in between the two curves describes a quasistable regime: systems which are stable for long periods of time, but eventually become unstable after a significant fraction of the Kozai timescale has passed.

Because the setup in this section does not account for nonzero eccentricities, which are very important for Kozai cycles, a more general theoretical framework is preferred. This discussion follows in a later section.

Comparison with simulation data

Before we consider systems with arbitrary eccentricities, we would like to validate our results from the previous section with simulation data. Hamilton & Burns (1991) use numerical methods to study the Hill stability of a massless particle orbiting an asteroid which, in turn, moves around the Sun. They determine retrograde orbits to be more stable, which agrees with our findings. In addition, they find a region of instability for inclinations around 90°, and note that particles in this region undergo oscillations between the direct and retrograde regimes. Although the authors do not mention it themselves, this is strongly indicative of Kozai cycles. In Figure 2, we show their data on our Hill-Kozai plot. Specifically, for inclination values taken in increments of 10°, we took the largest separation in the Hamilton & Burns data with the property that all smaller separations give stable orbits. Note that in the Kozai regime, the black dots lie above the dashed curve. This is because the simulation is time-limited to five asteroid orbits, or about twenty years. The points that we see in Figure 2 become Kozai unstable on timescales longer than this simulation time, and thus the true critical value for the separation in the simulation data is smaller than established. Letting the simulation run for a longer time, comparable to the Kozai timescale, would adjust the points in Figure 2 by placing them closer to the Kozai curve.

Figure 2. Earlier plot showing Hill-Kozai stability, now with simulation data from Hamilton & Burns (1991). Hill stability curve is a continuous line, with the dashed line indicating how the critical separation changes when Kozai cycles are accounted for. Data points from the simulation are shown in dots for 10° increments of the inclination. Each dot indicates the critical distance with the property that all smaller separations are stable. The simulation data contains several stable outliers about this critical distance; these points were disregarded, as they could be caused by too short of a simulation time and this could potentially be unstable. In the Kozai regime, we note that while the points do not exactly follow the dashed curve (see text for explanation), they fall significantly below the solid curve and hence suggest that a source of instability peculiar to this regime. We believe that this instability is caused by the Kozai mechanism.

There is another feature to be seen in the comparison between our analytical result and the simulation data. The region of instability around 90°, which we postulated to be induced by Kozai cycles, is shifted towards higher inclinations for the simulated data. The reason for this is unclear, because the simulated orbits start out on initially circular orbits, which is consistent with our assumption of small initial eccentricity. This discrepancy could be due to some outliers in the simulated data. See original figure in Hamilton & Burns (1991).

An alternative explanation of this trend comes from considering the octupole-level approximation, where the z-component of the angular momentum is no longer conserved. In such situations, it has been shown that the orbital motion can flip between retrograde to prograde, particularly for systems with high initial inclinations (Read et al. 2006). This kind of oscillation, due to tidal forces, is also observed in the simulations performed by Hamilton & Burns. In Figure 2, what we observe is an asymmetry in the flipping mechanism due to retrograde orbits being more stable. Specifically, prograde orbits that would otherwise be unstable become retrograde and thus more stable, shifting the Kozai regime to higher values of the inclination i. Confirming this hypothesis requires a better understanding of the oscillatory behavior between the prograde and retrograde states, which will not be addressed in this paper.

The Hill stability of eccentric systems

While Innanen’s analytical derivation, which was based on a two-body approximation, does not appear easily generalizable, a more recent paper by Li et al. (2010) investigates the problem of Hill stability for arbitrary eccentricities. They arrive at a closed-form solution for the critical semimajor axes ratio, which provides a sufficient, although not necessary condition for stability. Although, as mentioned before, some of the conclusions of this paper contradict what is widely regarded as true in the literature, we believe the analysis in this section serves as a good proof of concept for how closed-form Hill stability conditions can be generalized to include secular effects.

Consider a hierarchical triple system consisting of two masses m1 and m2 forming the inner binary, which together orbit a more massive body m3. Denote by e1, a1 the eccentricity and semimajor axis of the inner binary, and by e2, a2 the eccentricity and semimajor axis of the outer orbit. Let i be the mutual inclination between the two orbits. Then the formula established by Li et al. is the following:

(\frac{a_2}{a_1})_{cr} = \frac{\zeta}{\frac{9}{4}(1-e_2^2)x_0} \cdot \cos^2{(\frac{1}{3}\arccos{[\frac{3\sqrt{3}\cos{i}\sqrt{1-e_1^2}(3-4e_2^2)}{\zeta^{3/2}}]}) - \frac{\pi}{3}} (10)

where ζ = 9 + λ¹(1 + λ)2 e2² (x0² − 3x0) + 2(λ − 1)(1 + λ)-1x0, λ is the mass ratio of the binary elements m1/m2≥ 1, and x_0 = (\frac{m_1+m_2}{3m_3})^{1/3}. A sample plot for what this looks like as a function of the inclination when all the other parameters are set is shown in Figure 3. In this diagram, the stability region lies above the curve described by (a2/a1)cr. Having this very general analytical result, one would expect that Kozai could be added to the plot using the same procedure employed in the circular case. Note that the dependence on the inclination and eccentricity of the inner binary in the formula above is only through the quantity K = \cos{i}\sqrt{1-e_1^2}. In the test particle limit for the quadrupole approximation, this expression is constant. Therefore, with these simplifying assumptions, Kozai cycles alone cannot drive a system to instability.

Figure 3. Critical semimajor axes ratio as a function of inclination, for a triple system with m1/m2 = 1, (m1+m2)/m3 = 0.1, e1 = 0.2 (eccentricity of the inner orbit), and e2 = 0 (outer eccentricity). The Hill stable region lies above the curve. Up to a quadrupole approximation and for m2 small, this stability is valid even as the systems undergo Kozai oscillations.

It becomes clear, then, that going beyond these reductions is necessary. Indeed, even small variations in K may be able to drive specific systems to instability. In particular, close to the stability limit, the perturbations from the star are strong as the planet becomes unable to retain its satellite, and the quadrupole approximation likely becomes incorrect. This may explain why simulation results show instability in the Kozai regime, where the analytical derivations suggest otherwise, and why the stability condition we established for the circular case is different when secular effects are considered. Work has been done in characterizing secular perturbations at the octupole level (Ford et al. 2000). Using those results with stability criteria valid for arbitrary eccentricities may prove fruitful, although that is beyond the scope of our work.


The result established by Li et al. is, of course, more general than what we obtained based on Innanen’s work. However, the two disagree on whether retrograde orbits are more or less stable. Note that in the general case, since the stability region in Figure 3 lies above the curve for (a2/a1)cr, we see that the semimajor axis of the inner binary, a1, can go up to higher values in the case of prograde motion. This means that prograde orbits are more stable than retrograde ones: the Hill sphere where motion of the secondary binary component is permitted has a larger volume in the first case than in the second. However, for the results by Innanen we see the opposite to be the case. Based on the graph in Figure 1, the (normalized) Hill radius is larger for retrograde orbits than direct ones. This renders retrograde motion more stable, which is also explicitly pointed out in the original paper by Innanen.

One aspect that needs to be mentioned here is that like the majority of stability criteria, the two discussed above are sufficient, but not necessary conditions for Hill stability to occur. [3] This means that if one such criterion is satisfied, the system is Hill stable, though systems not satisfying the condition may be stable in the Hill sense as well. In Figure 4, we show the unscaled critical values for rin/rout from the Innanen paper, as well as a1/a2 from the Li et al. paper. Note that we reversed the a2/a1 ratio from Figure 3 so that the stable regions in both cases lie under the two curves. For this particular figure, we took (m1+ m2)/m3=0.1, and m1/m2=100 to be consistent with Innanen’s hypothesis that the secondary binary component is massless. The choice of these parameters does not affect the relative positions of the two curves, as the one corresponding to the criterion from Li et al. always lies entirely below the other one. This means that the results we derived for the circular case are stronger. This is somewhat expected, since the method developed in Innanen’s work was targeted to just one case, while the derivation in Li et al. aimed to give accurate results for any choice of eccentricities.

Figure 4. Comparing Innanen (top curve) and Li et al. (bottom curve) for e1 = e2 = 0. We took (m1 +m2)/m3 = 0.1, and m1 /m2 = 100, to be the consistent with Innanen’s hypothesis that the secondary component of the inner binary has negligible mass. The aspect of this figure doesn’t change if the parameters above are modifed, as long as m1/m2 >> 1. Hill stable regions lie under each of the two curves. We see that the two curves disagree on whether the prograde or retrograde regime is more stable.

While the two results do not directly contradict each other, it is baffling that they disagree on which kind of motion—prograde or retrograde—is more stable. The consensus in the literature is that retrograde motion is more stable, partly due to the Coriolis force. See for example Henon (1970), Zhang & Innanen (1988), or Hamilton & Krivov (1997). Furthermore, based on their results, Li et al. conclude that the observed Kuiper belt binaries are Hill unstable, which seems implausible. Lastly, the simulation results in Hamilton & Burns (1991) described in the previous section give additional credibility to Innanen over Li et al.

A different starting point for future work on secular stability criteria may reside in a different paper, by Donnison (2010). Although this work does not provide the condition for stability in a closed form, it does at least render retrograde orbits as more stable, making the results more reliable.

Applications to Kuiper belt binaries

Kuiper belt binaries (henceforth referred to as KBBs) provide excellent applications for the theory of Kozai cycles. First and foremost, these observed trans-Neptunian systems exhibit a large variety of inclinations and eccentricities (see Figure 5), and these orbital parameters allow for Kozai oscillations to occur (with the Sun acting as the perturbing object). The fact that they live within the Solar System is crucial for allowing us to observe highly inclined systems, thus removing a selection bias that would appear, for instance, in the case of exoplanets. Furthermore, it seems that highly eccentric or highly inclined Kuiper binaries are not an exotic species, but rather the norm across these type of systems. This, of course, suggests a formation mechanism different from the rest of the Solar System, which may involve successive collisions supplemented by Kozai evolution and tidal dissipation (Naoz et al. 2010). Interestingly enough, it is precisely these orbital parameters that give the most clue into the formation of these objects, and hence understanding how the parameters relate to each other is key.

Figure 5. The inclination-eccentricity distribution of the 22 Kuiper belt binaries used in our analysis. High inclinations and eccentricities are quite typical, making these objects ideal for the study of Kozai cycles.


In this analysis, we use data for 22 trans-Neptunian binaries published in Grundy et al. (2011). For each system, the data includes values for the semimajor axis, inclination, and eccentricity of the heliocentric orbit, as well as the inclination and eccentricity of the inner orbit. The total mass in the binary is given as well. As such, we have most of the necessary data for applying the theoretical results in the previous section.

Out of these 22 systems, 13 have a mirror ambiguity in inclination, meaning that it is not possible to tell whether the motion within the binary is prograde or retrograde relative to the heliocentric orbit. Specifically, this results in two possible values for the inclination, of the form i and 180° − i. The other orbital parameters are well-constrained despite this ambiguity.

For one part of the analysis, we use data compiled in Naoz et al. (2010). While similar to the Grundy data for the most part, it also supplies values for the diameter ratios of the binary components, which can be converted into mass ratios assuming that the two bodies have similar densities.

Hill stability of KBBs

Because of their high inclinations, most of the Kuiper binaries from the data we are using undergo Kozai oscillations. Note that the semimajor axes of these binary orbits don’t change through Kozai. However, an increase in eccentricity does lead to closer approaches at periastron, given by q = a(1 − e), and wider separations at the apocenter, equal to a(1 + e). Because of this, the secondary may potentially be pushed outside the Hill sphere of the primary, thus making the binary unstable.

Figure 6 shows this not to be the case, as the ratio a(1 + emax)/RHill never exceeds 0.1, and is below 0.04 for the majority of the considered systems. Here we used the values of RHill published by Grundy et al. Of course, it is not surprising that the systems observed today are well inside their corresponding Hill stable regions, because all unstable systems would have been ejected a long time ago. It is also interesting that the binaries do not fill the entire available phase space for a/RHill. As mentioned above, this ratio is typically well below 0.1, even if in theory any system with a/RHill <1 should be stable. In practice, this is not quite true, since KBBs are under the perturbing influence of other objects in the Solar System besides the Sun. This can reduce the effective Hill radius, making binaries with a/RHill close to 1 unstable.

Figure 6. Maximum separation during Kozai, relative to the Hill radius. The blue and green markers correspond to systems with a mirror ambiguity in inclination. Despite this uncertainty, the locations computed with one set of values for the inclination match very well those determined using the other possible values. The red markers indicate binaries with no ambiguity in inclination. All values are smaller than 1, indicating that the systems remain stable during their Kozai cycles. Furthermore, the lack of data points with values of a(1+emax)/RHill above 1.0 suggests that the actual stability domain is only a fraction of the entire Hill sphere. 0is could be because the Sun is not the only perturber acting on these binaries. 0e planets and other objects in the Solar System can also induce instabilities, and thus restrict the Hill stable regime.

In addition, note that as a(1 + emax)/RHill can never exceed 2a/RHill. Values of a/RHill above 0.2 are rarely observed—Petit et al. (2008) discusses one extreme example with a/RHill at roughly 0.29. Typically, observed systems remain inside their Hill spheres throughout the entire orbit, avoiding instability.

Tidal interaction in Kuiper belt binaries

A very simple plot that indicates the presence of tidal effects in the evolutions of these systems is shown in Figure 7. Specifically, the plot gives the period-eccentricity distribution of the 22 binaries considered. Note how the data points fall inside an envelope which is very close to the x axis for small eccentricities and then encompasses the other points inside a convex hull. Equivalently, for a given period interval, the maximum observed eccentricity is proportional to the period. This pattern can be explained by tidal interactions. Consider, for example, a system with a small orbital period. This correlates to a small semimajor axis, which in turn yields shorter tidal circularization timescales (Pont et al. 2011). Hence the system would have become circular by now, which is why it appears on the plot with eccentricity close to 0. In contrast, a Kuiper binary with a much longer orbital period has a correspondingly longer circularization timescale, which is why the eccentricity can go out to higher values.

Figure 7. Eccentricity-period distribution of Kuiper belt binaries. Systems with eccentricities close to 0 are all in the short period regime. Also note the general lack of systems with low eccentricities and large orbital periods. Such systems would typically not be affected by tidal friction, so their periods would not change significantly in their lifetimes. Because we don’t see these systems at the present, we conclude that when KBBs formed, higher eccentricities were more common than lower ones.

There is another feature to be noted in this figure. At large values of P, the distribution of eccentricities within the allowed envelope is biased towards higher rather than lower eccentricities. Indeed, looking at the entire dataset, all low eccentricity binaries have small periods, and thus they probably arrived in these nearly circular orbits due to tidal circularization. This suggests that the primordial eccentricity distribution of these systems may have been biased towards higher values.

There is a more rigorous way to determine the shape of the envelope encasing all these points, by considering the (scaled) semimajor axis in place of the period. Of course, the two quantities are very much correlated. Suppose we have a system with some fixed value for a/Rp, where Rp is the radius of the primary component. Then there is some minimum critical eccentricity e for which the pericenter distance a(1 − e) is small enough for tidal interactions to occur. The meaning of “close enough” here can be roughly approximated by the condition a(1 − e) = cRp, where c is some constant. This means that in a scatter plot showing a/Rp versus the eccentricity, the envelope is described by (a/Rp)(1 − e) = c, for some empirically determined constant c. A possible curve that suggests the general shape of the envelope is shown in Figure 8. For values of a/Rp greater than about 30, the corresponding systems fall inside the envelope. For zero eccentricity, all objects have a correspondingly small a/Rp, which is consistent with the idea that these orbits were circularized and brought closer together from their initial semimajor axis separation, due to tidal interactions. Lastly, note the two systems with significant eccentricities which do not fall within the envelope. Most likely, these orbits are, too, on the way to being circularized, but it could be that their Kozai timescales are longer. For instance, if the initial semimajor axis a was small (and comparable to what we see today—note how both systems are in the low a/Rp regime), then the Kozai period, which is proportional to 1/a3/2, would be long enough for the orbits to still be eccentric.

Figure 8. Envelope constraining possible eccentricities for a specified value for the semimajor axis. Rp is the radius of the primary component. The existence of this envelope means that there is a lack of binaries in the high eccentricity, low period domain. This is because a short orbital period correlates with a small physical separation within the binary, which leads to close approaches at the pericenter when high eccentricities are attained during Kozai cycles. In turn, this places the binary in the regime where tidal interactions become significant, which eventually circularize the orbit. Here, we also included an empirically determined fit for the envelope, given by a(1 - e) = 28Rp. This fit suggests tidal circulation to typically become effective at a separation of at most 28Rp. The outliers of this envelope could have somewhat different properties (e.g. mass ratio) from the rest, which affects the separation at which tidal circulation becomes important.

It should be noted that the envelope described above does not give a universally valid separation at which tidal effects become important. This separation also depends on various system parameters, such as the mass ratio of the two objects. Our goal here was to describe qualitatively how tidal effects may drive systems into the configurations we see today; investigating this idea more quantitatively could be a worthwhile direction for future work. It should additionally be noted that the idea of plotting eccentricity against some quantity that correlates with the circularization timescale has been explored before in the literature. Grundy et al. (2011) looks at the dependence of eccentricity on Msys17/6 γ4/ a13/2 (1+ γ3)7/3 where Msys=Mprimary +Msecondary, γ = Rsecondary/Rprimary, and a is the binary semimajor axis. By considering how the eccentricity depends on a simpler quantity, namely the scaled semimajor axis, our goal was to provide a more intuitive explanation for why tidal effects are important in the evolution of Kuiper binary systems.

A more theoretical overview of Kozai cycles and tidal friction is provided in Fabrycky & Tremaine (2007).

Timescale for tidal circularization

As previous work suggests, tidal dissipation plays a role in the evolution of other observed systems. Prior analysis of other such systems can be adapted to trans-Neptunian objects, which is what we will attempt to do here.

In the field of exoplanets, it is believed that the low eccentricities of hot Jupiters can be explained by a tidal circularization process following the inwards migration of these massive planets. This is addressed, for instance, in recent work by Pont et al. (2011), which is what we will use here to guide our discussion of tidal effects in Kuiper belt systems.

In the case of exoplanets, the planetary orbit is circularized due to tidal effects acting from the star on the planet, and not the other way around (Hansen et al. 2010, as cited in Pont et al.). Kuiper belt binaries can be made analogous to extrasolar systems by considering the primary object in the binary to act as the star, which exerts tidal forces on the secondary/planet. In Figure 9, we mimic the plot shown by Pont et al., which displays the mass ratio Mplanet/Mstar versus the scaled semimajor axis a/Rplanet. Now, the relevant quantities are Msecondary/Mprimary and a/Rsecondary.

Figure 9. Mass ratio versus scaled orbital distance, illustrating the tidal circularization timescale. Red circles correspond to orbits that are very close to circular (e < 0.01). The blue marker corresponds to an eccentricity of 0.1. Green squares correspond to higher eccentricities (e > 0.2). Shorter tidal circularization timescales correspond to smaller values of a/R secondary and Msecondary / Mprimary . The bottom-left of the figure therefore correlates to shorter timescales for orbits to become circular. It is in this region that the circular systems appear (red circles). This is because they already had enough time to circularize. Next in the progression from the faster to slower circularization times come to blue markers. 0ese systems are very close to the end of the tidal circularization process, but still have nonzero eccentricities. The other binaries (green squares) correspond to the slowest circularization times. The time evolution of the Kuiper systems on this plot is parallel to the x axis, since the mass ratio doesn’t change, and towards smaller values of a/Rsecondary. This evolution is exponentially faster for smaller initial values of the scaled semimajor axis, due to the dependence τD (a/Rsecondary)^5.

Theoretical models for tidal interactions provide a timescale for circularization to occur, given by

\tau = \frac{4}{63}Q(\frac{a^3}{GM_p})^{1/2}\frac{M_s}{M_p}(\frac{a}{R_s})^5 (11)

where Q is a factor measuring how responsive the secondary is to the tide raised by the primary component, Ms and Mp are the secondary and primary masses, a is the orbital semimajor axis, and G is the gravitational constant (Goldreich & Soter 1966, as cited in Pont et al.). This formula implies that the circularization timescale decreases with smaller values of the mass ratio Ms/Mp, and decreases steeply with smaller values of the semimajor axis, scaled to the secondary radius.

Figure 9 in our paper shows the mass ratios and scaled semimajor axis for the Kuiper binaries where Msecondary/Mprimary is known. The quantities along both axes carry information about the tidal circularization timescale, as per the formula above. Indeed, when we look at how systems in different eccentricity bins are distributed in this plot, there is pattern suggestive of tidal evolution. Orbits which are very close to circular, meaning they have eccentricities below 0.01, appear in the lower circularization timescale domain of the plot (red circles). This is because all such systems had enough time to become circularized by tidal interactions. In the next bin we placed the only system with eccentricity between 0.01 and 0.1 that has a known mass ratio (blue rhombus). Lastly, all systems with e > 0.2 are shown in green squares. They are visibly separated from the lower eccentricity binaries, appearing at higher values of Msecondary/Mprimary and a/Rsecondary. What this means is that the tidal circularization timescale is longer for all of these systems, and their nonzero eccentricities indicate that not enough time has passed for the orbits to become circular.

One observation that might better explain this figure is that the formation mechanisms are believed to be different for high and low mass ratio Kuiper belt binaries. These mechanisms are collisions for the former and dynamical capture for the latter. Although this will not be discussed here, it is something to keep in mind for future work on the topic.

We chose not to perform additional binning of eccentricities higher than 0.2 due to the limited data. However, with the inclusion of additional systems in the plot above, it is possible that one would observe a gradient in eccentricities going from the low to the high circularization timescale regime of the plot (i.e., from the bottom-left to the top-right). Additionally, one could look for patterns in the distribution of systems with similar eccentricities. With the current data, the three low eccentricity binaries with determined mass ratios (red circles) lie roughly on the same line, which may correspond to a specific tidal circularization isochrone. Although the current statistical significance of a pattern like this is low, with more data one could imagine fitting such isochrones and searching systems with comparable tidal circularization timescales for similar histories.

Kozai oscillations and tidal friction

At the beginning of our discussion on trans-Neptunian objects, we mentioned that Kozai cycles may yield small physical separations between the binary elements at the pericenter, due to increases in eccentricity.

This effect may yield close enough encounters between the binary elements to cause collisions or, in most cases, tidal interactions between the two objects. It is known that these tidal forces are able to circularize orbits, with the added effect that Kozai oscillations are dampened. Indeed, with Kozai alone, the quantity K = \sqrt{1-e_i^2}\cos{i} is conserved, leading to periodic changes in inclination and eccentricity between the same limits iminimax and eminemax; however, when tidal dissipation reduces the eccentricity to some etide without also increasing the inclination to keep K constant, the eccentricity and inclination will now Kozai oscillate with amplitudes dictated by the new K. This process repeats until the orbit is circularized and the system keeps some constant inclination ifinal, dictated by the last value of the quantity K.

In Figure 10, we plot Kozai oscillations for 18 systems in Grundy et al. (2011). The x axis displays cos i where i is the current inclination. On the y axis, we wanted a measure of how close the pericenter distance is to a separation at which tidal interactions are significant. For this reason, we scaled the quantity q = a(1 − e) to the Roche radius, which is the sphere surrounding the primary in the binary within which the tide raised on the secondary component would disintegrate it. The expectation is that tidal interactions become significant within several Roche radii from the primary. The formula used for the Roche limit is RRoche= (4Mprimary/ρsecondary)1/3 (Kutner 2003). Although the density of the secondary is not known, it is a reasonable assumption that ρsecondary = ρprimary, since two objects currently in a tight orbit are likely to have formed in similar conditions; additionally, photometric observations of typical trans-Neptunian binaries reveal matching colors for the primary and secondary, which further suggests that they have similar compositions. With this simplifying assumption, the Roche limit is R_{Roche} = R_{primary}\sqrt[3]{\frac{16\pi}{3}} since Mprimary = 4πRprimary3 / 3. This is what we’re using for our plot to scale the pericenter distance both for the current parameters of the systems, and throughout their Kozai evolution.

Figure 10. Kozai oscillations based on current values for the orbital parameters. Gray lines correspond to systems with no mirror ambiguity in inclinations. Lines of the same color indicate oscillations for the same system. Note that these lines are roughly symmetric with respect to the y axis, which is not surprising. We see that binaries at inclinations close to 90° have Kozai paths which intersect the line a/RRoche = 1. Even before the systems approach these points, tidal effects become noticeable. This should yield a higher concentration of systems close to i = 40° and i = 140° (see text). Although the figure does not contradict this claim, there are not enough data points to support it with a very high confidence.

The gray lines correspond to systems with uniquely determined inclinations, while the colored lines indicate binaries with a 180° uncertainty in inclination; the two solutions for such systems are shown in the same color, and are roughly symmetric with respect to the y axis. It is important to note that the plotted Kozai cycles are based on currently observed orbital parameters; any system which has interacted tidally at some point in its history may have started out with larger amplitude Kozai oscillations, due to the dampening effect described above.

This plot suggests a very specific evolution pattern of trans-Neptunian objects due to the combined effect of Kozai cycles and tidal interactions, which will be described in what follows. Consider first a binary whose initial inclination is below 40° (or above 140° for retrograde motion). In this regime, Kozai oscillations are ineffective, meaning that the system will keep its initial eccentricity and inclination, up to some small variations. Now, suppose instead the initial inclination is close to 90°. In this case, the amplitude of Kozai oscillations is very large, due to the quantity  being very close to 0. Then eccentricities close to 1 are probable, and these eccentricities drive the binary to periapsis separations where tidal dissipation becomes significant. At this point, subsequent Kozai oscillations are dampened due to the effect described earlier, and so the system now oscillates around a value for the inclination which is smaller than 90° (or larger, for retrograde motion). Lastly, systems with inclinations well between 40° and 90° (or 90° and 140°, respectively) will undergo Kozai oscillations with smaller amplitudes. Although some could cross inside the tidal dissipation regime, the rate of this happening is not as high as for systems with initial i around 90°. Indeed, it can be seen in our plot that binaries with higher initial inclination cross the line q/RRoche = 1 more readily than lower initial inclination systems.

We would like to consider how this evolutionary pattern produced the distribution of orbital parameters observed today. Assuming an initial uniform distribution of systems across different inclinations (at varying semimajor axes values), the discussion above implies the following: (1) binaries starting out with i < 40° or i > 140° will stay in the same regime due to the small amplitudes of Kozai cycles for these inclinations; (2) the majority of systems starting out with i around 90° (i.e. 90°−ε < i < 90°+ε, although good fits for ε are not addressed here) will be driven to lower inclinations (higher for retrograde motion) due to large amplitude Kozai oscillations, and will remain in that lower inclination regime due to dampening of Kozai cycles caused by tidal dissipation; (3) some binaries with 40° < i < 90°−ε or 90°+ε < i < 140° will be driven to lower inclinations by Kozai and kept there by tidal interactions, although many will continue to oscillate with little reduction in the Kozai amplitude. Therefore, after a time K = \sqrt{1-e^2}\cos{i} long enough for all these Kozai oscillations to occur, the initially uniform distribution in cos i will now present a depletion around i = 90° and a surplus in the regions where i approaches 40° from above, respectively 140° from below, while the regimes 0°-40° and 140°-180° keep their initial distribution.

Unfortunately, the current data sample we have is fairly small, which makes it difficult to test the theoretically predicted pattern above. Although more systems lie in the intermediate inclination regime, described by 0.5 ≤ |cosi| ≤ 0.7, than around i = 90°, the difference is not statistically significant. It is interesting to note, however, a general lack of systems at low inclinations (cos i close to -1 and 1). Since Kozai cycles and tidal circularization are not effective at removing binaries from this area, one could speculate that the primordial distribution of systems was not uniform, but rather skewed towards higher inclinations. However, given the small number of systems considered, it is not a low-probability event that few low-inclination binaries are observed. Additionally, binaries that are not edge-on, and thus have a considerable inclination, are easier to detect, introducing an observational bias.


In this paper, we started by deriving a closed-form Hill stability criterion by expanding on previous work by Innanen (1980). We then added secular effects to our results, and showed how our conclusions are in agreement with numerical simulations. More general conditions for stability, which allow for arbitrary eccentricities, were considered next. We concluded that a more involved analysis is necessary in this case, since Kozai cycles only affect the critical semimajor axes ratio in the octupole approximation. We also pointed out the disagreement of the two different criteria over which kind of motion is more stable—direct or retrograde.

In the second part of this paper, we worked with existent data on 22 Kuiper belt binaries to see what their current orbital parameters indicate about their evolution. We first showed all systems to be well within their corresponding Hill spheres even in their most unstable configuration during a Kozai cycle. We then explored how the Kozai effect can drive systems to tidal interactions, which in turn have the important effect of circularizing the orbit. This process explains an envelope observed in the period-eccentricity or semimajor axis-eccentricity plane. With additional data, a precise calculation of this envelope could establish the critical separation between the binary components at which tidal effects become important. The importance of Kozai cycles and tidal circularization in KBBs is further confirmed by showing a relationship between current eccentricity and the tidal circularization timescale for each system. We concluded this section by exploring what happens to the orbital parameters of KBBs as they undergo Kozai cycles and tidal circularization. A specific distribution of inclinations with a depletion around 90° was predicted, assuming a large, unbiased initial sample of binaries. Additional data is necessary to give more confidence to this prediction. Better constrained data will also allow for a more rigorous and quantitative analysis of these ideas.


I would like to thank my mentor, Dr. Hagai Perets for providing invaluable advice, guidance and support during the occasionally difficult, always intriguing research process. I am equally grateful to Dr. Ruth Murray-Clay for the constructive feedback on the initial draft. Lastly, I feel fortunate to have been part of an incredible learning experience in Astronomy 98, in no small part due to the enthusiasm and dedication of Prof. David Charbonneau.


Donnison, J. R. 2010, Planet. Space Sci., 58, 1169

Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298

Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385

Georgakarakos, N. 2008, Celestial Mechanics and Dynamical Astronomy, 100, 151

Grundy, W. M., et al. 2011, ArXiv e-prints

Hamilton, D. P., & Burns, J. A. 1991, Icarus, 92, 118

Hamilton, D. P., & Krivov, A. V. 1997, Icarus, 128, 241

Heggie, D., & Hut, P. 2003, The gravitational million-body problem (Cambridge University Press), 238-244

Henon, M. 1970, A&A, 9, 24

Innanen, K. A. 1980, AJ, 85, 81

Kozai, Y. 1962, AJ, 67, 591

Kutner, M. 2003, Astronomy: A Physical Perspective (Cambridge University Press)

Li, J., Fu, Y., & Sun, Y. 2010, Celestial Mechanics and Dynamical Astronomy, 107, 21

Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2011, ArXiv e-prints

Naoz, S., Perets, H. B., & Ragozzine, D. 2010, ApJ, 719, 1775

Petit, J.-M., et al. 2008, Science, 322, 432

Pont, F., Husnoo, N., Mazeh, T., & Fabrycky, D. 2011, MNRAS, 378

Read, J. I., Wilkinson, M. I., Evans, N. W., Gilmore, G., & Kleyna, J. T. 2006, MNRAS, 366, 429

Valtonen, M., & Karttunen, H. 2006, The Three-Body Problem (Cambridge University Press)

Zhang, S., & Innanen, K. A. 1988, Icarus, 75, 105

  1. [1]These numbers are actually determined assuming a small initial eccentricity. When this eccentricity is larger, the limits also change slightly. Nevertheless, the range 40°-140° is a good indicator for the domain where Kozai is effective.
  2. [2]Naoz et al. (2011) derives a more general expression for the maximum eccentricity, when the satellite has nonzero mass.
  3. [3]Note that the fact that a system is Hill stable does not mean that it is stable against other disruptive forces. In fact, prograde orbits are only stable up to half the Hill radius when everything else is accounted for.