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)
results.append(answer.T.copy())
resultz = []
for x0ydot0, thing in zip(x0ydot0s, results):
y = thing[1]
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], [0], 'ok')
plt.plot([x2], [0], '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()
