# What variables are needed to calculate simple horseshoe orbit times?

EDIT This was NOT a duplicate of Horseshoe orbit cycle times.

But that other question has been deleted, regardless. My original question asked for the answer to the equation(s), and provides many of the variables. But since no one was apparently willing/able to answer it, this new question asks only for the form of the equation(s) itself and a list of variables involved, so that I can either provide more information on a future question, or try and solve the equations myself (though I doubt my math ability to do that, hence the original question, if I know what equations and variables are involved I can at least make an attempt myself).

What variables need to be known in order to calculate a horseshoe orbit cycle time? In other words: What is the equation, and what do the variables in that equation represent (speed? mass? semi major axes? etc.), to calculate the time it takes for an object to go from Point A on that image, through Points B, C, D, and E and back to A, along the light blue line in the image above (image also available here)?

For example, Earth and Cruithne complete their cycle in 770 years, while Janus and Epimetheus complete theirs in 8 years. I'm confident that semi-major axes are key factors in the equation, but I don't know what other variables are included, or how the variables relate to each other in the calculation.

Also, I understand that these orbits are not stable, in the long term, and that my example of Cruithne is especially complex, as it will actually change orbit types periodically, and will likely be ejected from the solar system entirely, or impact the Sun or Jupiter, at some point. But those complexities are for another time. For this question, I just want to know the equation in its simplest form.

• That's not really a plot of a horseshoe orbit, it's a equipotential line for a pseudopotential in a rotating frame for an object with zero velocity in that frame. If the object were in a real horseshoe orbit, those potential lines would change. I wonder if you can find an image of an actual horseshoe orbit?
– uhoh
Feb 9, 2019 at 2:26
• See also Horseshoe orbit cycle times which also incorrectly uses that image and currently has no answer. Also unanswered is Horseshoe orbits and integration in C
– uhoh
Feb 9, 2019 at 2:28
• Wait, Horseshoe orbit cycle times is your question, which you've also posted in Physics SE. How is this question different than those?
– uhoh
Feb 9, 2019 at 2:31
• that makes sense. I think the issue of the cycle time of horseshoe orbits is really interesting. I'd like to propose a course of action. Why don't you temporarily delete your other question to address any concerns of duplication, and let me do a bit of a rewrite on this question. Then I'll vote to re-open it. I'm planing on visiting the library tomorrow and I can do a search on this and see what analytical models there have been published.
– uhoh
Feb 11, 2019 at 15:11
• @peterh If/when I figure out how to do these calculations, I'll use that information in a way that doesn't require long term stability in geological timescales. I'm aware of their inherent instability and have taken that in to consideration for my final goals. But thanks for the concern and pointing it out. Feb 11, 2019 at 19:15

That depends on the accuracy you want to work with. To zeroth order, as outlined in Murray & Dermott, "Solar System dynamics", Chapt 3., you can do the following:

• The zero-velocity contours that are plotted in your image will not be coinciding with particle orbits to infinite precision, but they're a good zeroth order approximation for objects with low eccentricity w.r.t the star ($$v_{r}/v_{\theta} \ll 1$$)
• A particle on those orbits is on a regular Keplerian orbit of radius $$r_{\rm H}$$, out of reach of the gravitational influence of the planet at $$r_{\rm P}$$. Thus, to get the largest "chunk" of traveltime of one part of the horseshoe orbit, either inner or outer part, you can work out with relative velocities and the assumption of keplerian velocities.
• Caution is required as to which orbital time you're interested in: If $$a$$ is Earth's distance to the sun, and the horseshoe is $$d$$ away from a perfectly circular orbit, thus alternating between distances $$a\pm d$$, then the orbital time in the rest-frame of the sun will be $$v_{-}=\sqrt{\frac{GM}{a-d}}$$, for a large number of orbits until close encounter, and after that $$v_{+}=\sqrt{\frac{GM}{a+d}}$$.
• Therefore, the simplified case that the satellite is on its own Keplerian orbit is true for most of the time. Armed with this knowledge, we can approximate the Planetocentric recurrence time as $$t_{\rm rec}=\frac{2\pi a}{v_{K} - v_{-}}$$ as a simple catch-up time between objects on different orbits.
• The relative velocity $$v_{\rm rel} = v_{K} - v_{-} = \sqrt{\frac{GM}{a}}-\sqrt{\frac{GM}{a-d}}$$ can be expanded in the limit of $$d/a\ll1$$ into $$v_{\rm rel} \approx \frac{1}{2} \sqrt\frac{GM}{a} \frac{d}{a}$$ and thus I derive from this $$t_{rec}\approx \frac{4\pi}{\sqrt{GM}} \frac{a^{5/2}}{d}$$. As expected, the cycle time diverges for $$d=0$$, as in the co-orbiting case with Earth this time must be infinite. An upper limit on $$d$$ cannot be derived from this, one would need to turn for the full solution for that.

Just out of curiosity, I've plugged in a few values into this formula and written something quick in python:

import numpy as np

#Basic physics quantities
G      = 6.678e-8 #cgs units
pi     = 3.141592
navo   = 6e23 # particles per mole
sigma  = 5.67e-5   #erg cm-2 s-1 K-4
kb     = 1.38e-16  #erg/K
km     = 1e5 #kilometers in cm
mearth = 5.98e27  #g
msun   = 2.0e33   #g
au     = 1.49e13  #cm
yr   = 365*24*3600
rearth   = 6370e5
rjupiter = 74000*km

#
# Returns the approximate horseshoe-cycle time in the reduced 3body problem
# Masses of bodies: m0>>m1>>m2
# Semimajor axis distance is from m0 to m1, radial distance is a(m0->m1)-a(m0->m2)
#

#
# https://en.wikipedia.org/wiki/(419624)_2010_SO16 around the Sun
#
# Quoted cycle time ~350 years, with d=0.004 AU
#
print("Predicted 2010_SO16 cycletime [years] = ", hs_cycle(msun, 1*au,0.004*au), " predicted = 350 yrs")

#
# Janus/Epimetheus around Saturn
#
# a = 151410 km, d = 25 km, as stated in https://en.wikipedia.org/wiki/Epimetheus_(moon)
# Quoted cycle time = 8 years (from comments)
#
print("Predicted Janus/Epimetheus cycletime [years] = ", hs_cycle(95*mearth, 151410*km,50*km), " predicted = 4 yrs")

#
# 3753 Cruithne
#
# a = 1 AU and semimajor axis difference from https://en.wikipedia.org/wiki/3753_Cruithne
# Quoted cycle time = 770 years
#
print("Predicted 3753 Cruithne cycletime [years] = ",hs_cycle(msun, 1*au, (1.0-0.99774)*au), " predicted = 770 yrs")


and the results I get are

Predicted 2010_SO16 cycletime [years] =  495.7747141830971  predicted = 350 yrs
Predicted Janus/Epimetheus cycletime [years] =  11.542076781209305  predicted = 8 yrs
Predicted 3753 Cruithne cycletime [years] =  877.4773702355546  predicted = 770 yrs


So the formula may be off by up to a factor of ~2. This is for sure simply because reality is more complex than a simple circular-orbit approximation, but also due to the quality of the values used. Wikipedia isn't well known for researching particular values well. I took those that I found there. For SO16 that was particularly confusing to select, so I took those two that were mentioned in the same line of text, hoping they would come from the same source.

Anyone finding more consistent values, is free to comment.

• While I understand that accuracy is an issue, to some extent, I'm having trouble reconciling this equation with the examples I provided of Earth+Cruithne and Janus+Epimetheus. I'd understand a deviation of a few decades on a 770 year cycle, or several months or a year on an 8 year cycle. But I get values of hundreds of thousands of years on one, and a couple of weeks on the other. What am I missing here? Feb 12, 2019 at 18:12
• @Dalila: I'm not getting those values. Mine seem much more reasonable, but it's still an approximation. You have to be careful on what is the semimajor axis and the value d, which is the difference of semimajor axis and the 2nd massive body. The central object is the Sun for Cruithne and Saturn for JEpi. But I've added my code, so you can check for yourself, and of course dig into the literature if you wish. Feb 12, 2019 at 19:10
• I think this is where I ran in to trouble: "You have to be careful on what is the semimajor axis and the value d, which is the difference of semimajor axis and the 2nd massive body." Can you clarify that value? did you mean it's the difference between the semi-major axes of the two co-orbital bodies? I just want to make sure I understood that last comment correctly Feb 12, 2019 at 19:35
• @Dalila: If we denote the bodies by 0,1,2, and the masses are $m_2 \gg m_1\ gg m_0$, then d is the difference of semimajor axes $a_{01}-a_{02}$. Note that this shouldn't exclude Janus and Epimetheus, as they obey the same physics as a regular Horseshoe-orbit: Each moon individually orbits Saturn unperturbed, for most of the time. So I'm not sure what's wrong with the value I get for JEpi. Feb 12, 2019 at 20:05
• Thanks for that clarification. I understood the JEpi cycle to be 8 years, total, not 4. 4 years was the time between the two approaches in the cycle, so 11 vs 8 is actually a reasonable margin of error, at least for my purposes, on that calculation. Did I misunderstand something else on that one? Feb 12, 2019 at 20:11

Ernest W. Brown's On a New Family of Periodic Orbits in the Problem of Three Bodies: (Plates 6, 7.) in MNRAS, 71, (5), pp 438–454 published on 10 March 1911 seems to be where horseshoe orbits were first proposed. (Available here as well). It begins:

There are four known asteroids which appear to oscillate about one or other of the vertices of the two equilateral triangles which have the line joining Jupiter and the Sun as base. These vertices are the well-known positions of relative equilibrium. The heliocentric vector of one of these asteroids can apparently move as far as 17° away from its equilibrium position.* The oscillations cannot therefore be considered very small. One naturally asks whether oscillations of this kind in arcs of still greater extent are possible ; and if so, in what manner the orbits may he most conveniently obtained.

*L. J. Linders, Arhivfor Mat., Ast. och Fys., So. Vet. Ak. i Stockholm, Bd. 4, No. 20.

I'll make some horseshoe orbits in the Circular Restricted Three-Body Problem formalism and plot them in Python, then compare to the synodic period estimation described in @AtmosphericPrisonEscape's answer.

tl;dr: There's good qualitative agreement, no surprises!

A brief summary of CR3BP math in dimensionless units. The distance between the two bodies is equal to 1, as is the gravitational constant. They orbit around a common center of mass in circular orbits, with a period of $$2 \pi$$. It's easier to visualize and calculate if you do it in a rotating frame, so the two masses are fixed. The third body at position $$x, y, z$$ is considered to have no gravitational effect on the first two,

$$\mu = \frac{m_2}{m_1 + m_2}$$

$$x_1 = -\mu$$ $$x_2 = 1-\mu$$

$$r_1 = \sqrt{(x-x1)^2 + y^2 + z^2}$$ $$r_2 = \sqrt{(x-x2)^2 + y^2 + z^2}$$

The Jacobi Energy $$C$$ is a conserved quantity in this rotating frame:

$$C = x^2 + y^2 + 2\frac{1-\mu}{r_1} + 2\frac{\mu}{r_2} - (\dot{x}^2 + \dot{y}^2 + \dot{z}^2)$$

where the $$x^2 + y^2$$ is the pseudopotential. If you set the velocity dependent terms $$(\dot{x}^2 + \dot{y}^2 + \dot{z}^2)$$ to zero, you get a zero velocity surface, that surface that's pasted into many/most questions about three body orbits. These plots do not apply when an object is moving, and so you can't superimpose orbits on top of them! The acceleration felt by the third body in this rotating frame has both the expected $$1/r^2$$ forces and a velocity-dependent pseudoforce which is not real, but accounts for the fact that the frame is rotating and not inertial.

$$\ddot{x} = x + 2\dot{y} - \frac{(1-\mu)(x+\mu)}{r_1^3} - \frac{\mu(x-1+\mu)}{r_2^3}$$ $$\ddot{y} = y - 2\dot{x} - \frac{(1-\mu)y}{r_1^3} - \frac{\mu y}{r_2^3}$$ $$\ddot{z} = -\frac{(1-\mu) z}{r_1^3} - \frac{\mu z}{r_2^3}$$

Here are some calculations. I chose $$\mu = 0.001$$ which is pretty close to the situation Jupiter and the Sun. I chose an array of starting points at the opposite point from $$m_2$$ at about $$x=-1$$ but that's not what I really did. What I really did is choose a bunch of starting velocities $$-0.08 < \dot{y} < 0.08$$ and for each I calculated the position on the $$x$$ axis near $$x=-1$$ where the acceleration in the $$x$$ direction was zero.

That gives the solutions a tiny bit of starting symmetry, but halo orbits are bumpy and wiggly and not always so stable, so this effort isn't really necessary.

I propagated each orbit until it came back around to the same area and stopped it when it crossed the x-axis, producing a family of half-cycles.

To make a long story short, the method shown in @AtmosphericPrisonEscape's answer of estimating the cycle time by calculating the synodic period in the inertial frame is in pretty good agreement with these halo orbits, and that shouldn't be very surprising! above: half-cycles of some wobbly horseshoe orbits above: times to first x-axis crossings of the same wobbly horseshoe orbits, used to calculate half-cycle times. above: cycle times from this calculation (black dots) versus from the synodic period estimation method (red dots). Good qualitative agreement. Also the starting y velocities at each starting point in x.

below: Python script for these plots.

def x_acc(x, ydot):
r1    = np.abs(x-x1)
r2    = np.abs(x-x2)
xddot = x + 2*ydot  -  ((1-mu)/r1**3)*(x+mu) - (mu/r2**3)*(x-(1-mu))
return xddot

def C_calc(x, y, z, xdot, ydot, zdot):
r1 = np.sqrt((x-x1)**2 + y**2 + z**2)
r2 = np.sqrt((x-x2)**2 + y**2 + z**2)
C = (x**2 + y**2 + 2.*(1-mu)/r1 + 2.*mu/r2 - (xdot**2 + ydot**2 + zdot**2))
return C

def deriv(X, t):
x, y, z, xdot, ydot, zdot = X
r1 = np.sqrt((x-x1)**2 + y**2 + z**2)
r2 = np.sqrt((x-x2)**2 + y**2 + z**2)
xddot = x + 2*ydot  -  ((1-mu)/r1**3)*(x+mu) - (mu/r2**3)*(x-(1-mu))
yddot = y - 2*xdot  -  ((1-mu)/r1**3)*y      - (mu/r2**3)*y
zddot =             -  ((1-mu)/r1**3)*z      - (mu/r2**3)*z
return np.hstack((xdot, ydot, zdot, xddot, yddot, zddot))

# http://cosweb1.fau.edu/~jmirelesjames/hw4Notes.pdf

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
from scipy.optimize import brentq

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

mu = 0.001

x1 = -mu
x2 = 1. - mu

x = np.linspace(-1.4, 1.4, 1201)
y = np.linspace(-1.4, 1.4, 1201)

Y, X = np.meshgrid(y, x, indexing='ij')
Z    = np.zeros_like(X)

xdot, ydot, zdot = [np.zeros_like(X) for i in range(3)]

C = C_calc(X, Y, Z, xdot, ydot, zdot)
C[C>8] = np.nan

if True:
plt.figure()
plt.imshow(C)
plt.colorbar()
levels = np.arange(2.9, 3.2, 0.04)
CS = plt.contour(C, levels,
origin='lower',
linewidths=2)
plt.show()

ydot0s   = np.linspace(-0.08, 0.08, 20)
x0ydot0s = []
for ydot0 in ydot0s:
x0, infob =  brentq(x_acc, -1.5, -0.5, args=(ydot0), xtol=1E-11, rtol=1E-11,
maxiter=100, full_output=True, disp=True)
x0ydot0s.append((x0, ydot0))

states = [np.array([x0, 0, 0, 0, ydot0, 0]) for (x0, ydot0) in x0ydot0s]

times  = np.arange(0, 150, 0.01)

results = []
for X0 in states:
answer, info = ODEint(deriv, X0, times, atol = 1E-11, full_output=True)

resultz = []
for x0ydot0, thing in zip(x0ydot0s, results):
y     = thing
check = y[2:]*y[1:-1] < 0
zc    = np.argmax(y[2:]*y[1:-1] < 0) + 1
if zc > 10:
resultz.append((thing, zc, x0ydot0))

if True:
plt.figure()
hw = 1.6
for j, (thing, zc, x0ydot0) in enumerate(resultz):
x, y = thing[:2,:zc]
plt.plot(x, y)
plt.xlim(-hw, hw)
plt.ylim(-hw, hw)
plt.plot([x1], , 'ok')
plt.plot([x2], , 'ok')
plt.show()

if True:
plt.figure()
for j, (thing, zc, x0ydot0) in enumerate(resultz):
x, y = thing[:2]
plt.plot(times[:zc], y[:zc])
plt.show()

if True:
plt.figure()
for j, (thing, zc, x0ydot0) in enumerate(resultz):
x0, ydot0 = x0ydot0
cycle_time = 2. * times[zc] / twopi
ratio = abs(x0/x2)
T_simple_model = twopi * abs(x0/x2)**1.5
T_synodic_simple_model = 1. / (1. - twopi/T_simple_model) # https://astronomy.stackexchange.com/a/25002/7982
plt.subplot(2, 1, 1)
plt.plot(x0, cycle_time, 'ok')
plt.plot(x0, abs(T_synodic_simple_model), 'or')
plt.subplot(2, 1, 2)
plt.plot(x0, ydot0, 'ok')
plt.subplot(2, 1, 1)
plt.xlabel('x0', fontsize=16)
plt.ylabel('cycle times (periods)', fontsize=16)
plt.subplot(2, 1, 2)
plt.xlabel('x0', fontsize=16)
plt.ylabel('ydot0', fontsize=16)
plt.show()

• Hah, that's nice. Yesterday I was also working on my own three-body solver to test the formula, but you beat me to it. Feb 14, 2019 at 13:50
• @AtmosphericPrisonEscape I've been dabbling in this for a very long time but have a long way to go.
– uhoh
Feb 14, 2019 at 14:29