Here's an animation I created showing the Kozai Mechanism at work in the Alpha Centauri system.

Two hypothetical planets orbit Alpha Centauri A at distances of 1 AU and 2 AU. Their orbits are inclined 60 degrees from the orbital plane of Alpha Centauri A and Alpha Centauri B. This shows about 1350 years, with a screen shot being generated every 2.7 years ~ the period of the outer planet. The time step used in the simulation was 2048.

Alpha Centauri B is not visible in this animation, but its presence is certainly felt. The orbit of the outer planet trembles about every 80 years in response to periastron, when the A/B distance decreases to only 11 AU.

The inner planet is doing the same thing as the outer planet, except much slower. I gave each planet 1 Earth Mass, so once they get out-of-plane, you've got a double Kozai going on. I'll try it with Mercury-Mars to see what happens.

Here's the graph of the full sim that produced the animation. The animation only takes you up to the outer planet's first max ecc. The graph covers 10x as much time. I stopped the graph there because The outer planet's eccentricity exceeded 1 and it escaped the system.

BTW, what happens if you have a low mass object in an external orbit (eg. swap the stellar companion and planet)? Does the planet get affected by the companion in the same way?

That's an interesting question. I'd have to modify the code to figure it out. Although Gravity Simulator gives you the ability to create an object with orbital elements with respect to a system barycenter, it only outputs orbital elements with respect to a single object. So the data wouldn't be any good. But at least the plot might yield a clue. I've got it running now.

BTW, how do the periods of the oscillations that you get compare to the Pkoz calculation? Do they match, or are they different?

Interesting question : out of the sims of the 5 earth-moon I always get a "match" , but not totally . Where the period should be rising with distance^3 according to the formula I get a period rising with distance^3.4 , regardless of the time-step I use ( steps between 64 and 512 ) . The factor I become ist always between 3.3 and 3.45 . This means a difference of about 10% .

I also noticed that the period is not exact or equal between two peaks . This doesn't surprise me because peaks (ecc) occur at slightly different incl and LAN .

BTW, how do the periods of the oscillations that you get compare to the Pkoz calculation? Do they match, or are they different?

Outer Planet: If you look at the graph, it takes 1730 years to reach max eccentricity from its initial 0 eccentricity. Multiplying this by 2 gives 3460 years as a period. The computed value is 3400. Not a bad match since I'm just eyeing the graph for dates.

But then following that, the period gets shorter. You'll notice that from the starting point where eccentricity was set to 0, it rose slowly at first. But after bottoming out at the first valley, it doesn't reach eccentricity = 0. And it rises faster. Consider the peaks centered on years 1736 to 6920. It does 3 complete cycles. (6920-1736)/20 = 1728 years per cycle, almost half of what is expected.

Inner Planet: From the graph, it takes about 5400 years to reach max ecc from its initial starting ecc of 0. Doubling this gives 12800 years. But again, eccentricity takes a long time to start picking up from ecc=0. Peak-to-peak is only about 6500 years, again about half the original value. The computed value is 9800 years.

I notice the same thing going on in Figure 2 (page 17) in Takeda's paper.

0.9 solar mass companion: Time from begining to first peak is 16 million years. Doubling gives 32 million years. But time from 1st peak to 2nd peak shortens considerably to 22 million years. Considering all 15 cycles on the graph, it has a time averaged period of 20.5 million years. The computed value is 25.6 million years.

0.08 solar mass companion: Time to first peak = 186 million years. Doubling this gives 372 million years. But the time from the first peak to the 2nd valley = 120 Myrs. Doubling this gives 240 Myrs. The computed value is 288 Myrs.

So it seems that the period is not the same from cycle to cycle and depends heavily on how close to ecc=0 the valley comes; the closer, the longer the period.

This sort of makes sense then the author would use the asymptoticlly equal to symbol. It's as if to say don't trust this formula for any one specific period, but as cycles approach infinity, the time-averaged period should approach the computed period.

BTW, what happens if you have a low mass object in an external orbit (eg. swap the stellar companion and planet)? Does the planet get affected by the companion in the same way?

This is the 2nd time I'm quoting this sentence. I started with no eccentricity in either the planet or the 2nd star. But other than that, the sims are the same with the planet and star reversed. The results are interesting enough for me to make an animation. Stay-tuned. This'll have to run overnight.

Well, I never did make the animation. But I can see what is going on. The planet, inclined 60 degrees from the barycenter, with initial ecc=0 remained 60 degrees inclined and ecc=0 for the dureation of the simulation. However, its longitude of ascending node advanced non-stop, completing 360 degrees in about 16000 years.

I placed Mercury, Venus, Earth & Moon, & Mars around Alpha Centauri A, and also around Alpha Centauri B. Since ACA it is 1.1 times as massive as the Sun, I also added sqrt(1.1) to their velocities so their orbits would retain their solar orbit properties. Likewise, since ACB is 0.907 times as massive as the Sun, I subtracted sqrt(0.907) from their velocities.

After only 1800 years at a nice slow time step of 64: Alpha Centauri A