A stroll through the kozai python package

Installation

The kozai package is available on PyPI and can be installed with pip like so:

pip install kozai

If you don't have the right permissions, try installing it like this:

pip install --user kozai

If you run into problems with the package tests, please let me know at antognini@astronomy.ohio-state.edu so we can work to fix them. Note that the kozai package only works on Python 2 at the moment.

Making your first triple

The kozai package makes it easy to create and evolve triples. The package comes with two classes which represent triples in different ways: TripleDelaunay and TripleVectorial. The TripleDelaunay class represents the triple using the Delaunay orbital elements and the TripleVectorial class represents the triple using the eccentricity and angular momentum vectors. At the moment, the TripleVectorial class only supports triples with a test particle secondary, whereas TripleDelaunay requires that all three stars have mass.

The TripleDelaunay class currently has more features than the TripleVectorial , so I recommend using the TripleDelaunay class wherever possible. If you need to evolve a triple in the test particle approximation you can simply set the secondary mass to be something very small (e.g., $10^{-5}$ $M_{\odot}$).

To start with we will use the TripleDelaunay class and explore the TripleVectorial class later. We first have to import the class from the kozai.delaunay module:

To create a triple we just call TripleDelaunay:

This triple has some default orbital parameters. We can see what the inner semi-major axis is:

Note that the semi-major axes are all in AU. Let's take a look at the inclination:

Note that the angles are given in degrees.

Evolving a triple

Now that we have a high-inclination triple we can evolve it in time and see Kozai-Lidov oscillations. To start with, let's turn the octupole term off so that we are just looking at quadrupole oscillations:

Now let's evolve the triple for $3 \times 10^4$ yr:

evolution is a numpy array where each entry is a different step in time. Each entry has the format

• $t$ (yr)
• $a_1$ (AU)
• $e_1$
• $g_1$ ($^{\circ}$)
• $a_2$ (AU)
• $e_2$
• $g_2$ ($^{\circ}$)
• $i$ ($^{\circ}$)]

Let's now see how the eccentricity changes in time.

We see here some nice Kozai-Lidov oscillations. Let's look at how the inclination changes in time:

We see that we start out at $80^{\circ}$ and then move to close to the Kozai critical angle of $39.2^{\circ}$.

The Eccentric KL Mechanism

Let's now take a look at what happens when we turn on the octupole order term. First we have to reset our triple:

Now the triple has been reset to its initial conditions. So the time is set to 0:

And the eccentricity and inclination are both at their original values:

Let's now turn the octupole term on:

Since we're looking for longer-term changes to the KL cycles we'll evolve this system for $10^5$ yr.

Now let's take a look at the eccentricity evolution:

It's exactly the same as before! Why is that? Let's take a look at the masses of our stars:

By default our triple has an equal mass inner binary, so the octupole order term is zero! We can see this explicitly by looking at the $C_3$ coefficient:

But the $C_2$ coefficient is not zero:

Note that the unit here is Joules.

In order to see the eccentric Kozai mechanism we need to have an unequal mass inner binary, so let's reset the triple and set the mass of the secondary to be something small, like $10^{-5}$ $M_{\odot}$.

And now we'll take a look at the eccentricity evolution:

We can see the characteristic increase in the maximum eccentricity during each KL oscillation as the orbit of the inner binary flips from prograde to retrograde. Let's see what's happening to the orbital inclination:

Some other features of the TripleDelaunay class

Setting the orbital parameters

The TripleDelaunay class comes with a number of other features to make it easy to work with. You can set the orbital parameters of the triple from the call to the class:

Printing the triple parameters

It is also easy to get general information about the triple. If print is called the class returns the various parameters of the triple in JSON format:

If you're working with many triples this makes it easy to keep track of the parameters of each one.

Getting the eccentricity extrema

If you're evolving a triple over a long period of time you may not be interested in most of the evolution. Instead you may just want the points at which the eccentricity is at its largest or smallest. To only get the eccentricity extrema the TripleDelaunay class has an extrema method:

Let's compare the extrema to the full evolution:

Finding flips

Over even longer timescales we may just be interested in the times that the inner orbit flips from prograde to retrograde or vice versa. To find these moments TripleDelaunay has the find_flip method:

Let's compare to the full evolution.

Note that this method just looks for a change in the sign of $\cos i$. For massive secondaries the "pole" in the Hamiltonian at which arbitrarily large eccentricities are possible occurs at inclinations slightly larger than $90^{\circ}$. find_flip will just find when the orbit crosses $90^{\circ}$, not when the arbitrarily large eccentricities occur.

The vectorial module

As mentioned at the beginning, the kozai package offers two classes to evolve hierarchical triples, of which TripleDelaunay is currently more powerful. But in cases where the test particle approximation is used the TripleVectorial class may be preferable. The two classes have been designed to have APIs as similar to each other as possible. Thus the methods between the two classes are mostly the same. The major differences are as follows:

• The mass of the secondary cannot be set in TripleVectorial --- it is always 0. Similarly, the mass of the secondary cannot be set to 0 in TripleDelaunay (although it can be set to some very small non-zero value).
• Post-Newtonian terms are not supported in TripleVectorial yet.
• The hexadecapole term is not included in TripleVectorial yet.
• vectorial has an additional method called flip_period which can be used to numerically calculate the time between flips. This will be described in more detail a little later.

Evolving a triple in vectorial

Triples can be made in TripleVectorial in much the same way as in TripleDelaunay:

We can examine the eccentricity evolution in the same way as before:

And similarly we can examine the evolution of the inclination:

And just as with TripleDelaunay we can just take the eccentricity maxima:

The flip_period method

The TripleVectorial class has one unique method: flip_period. This method will integrate along until some specified number of flips have occured (by default three) and will return the average time between flips. This method was not included in TripleDelaunay because arbitrarily large eccentricities do occur when the triple passes through some inclination slightly above $90^{\circ}$ which is dependent in a non-trivial way on the orbital and physical parameters of the system. In the test particle case this critical inclination is always $90^{\circ}$ so it is easy to know when arbitrarily large eccentricities are possible.

The method can be used very easily:

Let's compare this to evolving the triple directly.

So the flip period looks to have been calculated correctly!

The kl_period module

The kl_period module has a number of utilities to calculate the timescales in hierarchical triples. To start, let's import the module, create a triple, and calculate the orbital times:

We can also use the standard KL timescale to estimate the period of KL oscillations:

$$t_{\textrm{KL}} \sim \frac{8}{15 \pi} \left(1 + \frac{m_1}{m_3}\right) \left(\frac{P_{\textrm{out}}^2}{P_{\textrm{in}}}\right)\left(1 - e_2^2\right)^{3/2}$$

Let's compare this to the evolution of the triple.

So the timescale is off by a factor of a few. For TripleVectorial objects we can calculate the period exactly, however, using equation (33) of Antognini (2015), which is implemented in kl_period:

This matches the observed period much more closely. We can verify this by evolving the triple and calculating the mean time between eccentricity maxima using the numeric_kl_period function:

Evidently the numeric period does not too closely match the analytic period. This is because the length of a KL oscillation varies from one oscillation to the next due to the octupole order term. If the octupole term is turned off, these variations should vanish and we should match the result from kl_period much more closely:

One final feature of the kl_period module is the is_librating function. As its name implies, it returns True if the triple librates and False if it rotates:

The ekm module

The last module we'll explore is the ekm module. This module presents a new class called TripleOctupole. This class makes it possible to evolve a hierarchical triple in the test particle limit by averaging over individual KL oscillations so as to capture only the contribution from the octupole term.

We can evolve this triple in time. The evolve method returns the following:

• $t$ (in units of the secular timescale)
• $j_z$
• $\Omega$ (in degrees)
• $\left< f_j \right>$
• $\left< f_{\Omega} \right>$
• $x$

See Katz et al. (2011) for definitions of these quantities.

We can now see the change in the $z$-component of the angular momentum as the inner binary undergoes a few flips:

As with TripleDelaunay and TripleVectorial we can print out the properties of the triple in JSON format:

Finally, we can calculate the time between flips two different ways. We can either do it analytically using equation (67) of Antognini (2015) using the flip_period method:

We can also calculate the flip period numerically by evolving the triple directly using numeric_flip_period:

As we would hope, these two calculations yield very similar results.