Two Body Problem
Two body Problem
From earlier equations, we have
\[ \ddot{\mathbf{R}}_2 - \ddot{\mathbf{R}}_1 = \frac{G m_1}{r_{21}^{3}} \mathbf{r}_{21} - \frac{G m_2}{r_{12}^{3}} \mathbf{r}_{12} = -G (m_1 + m_2)\,\frac{\mathbf{r}_{12}}{r_{12}^{3}} \]
which is clearly,
\[ \ddot{\mathbf{r}}_{12} + \frac{G (m_1 + m_2)}{r_{12}^{3}} \mathbf{r}_{12} = 0 \]
Dropping the subscript, we have
\[ \ddot{\mathbf{r}} + \frac{\mu}{r^{3}} \mathbf{r} = 0, \qquad \mu = G (m_1 + m_2) \]
Here \(\mu\) is called the standard gravitational parameter. It is a measure of the strength of the gravitational field of the two-body system. So to analyse the two-body problem, we don’t need to know the individual masses, but only this constant.
Note
In the early development of celestial mechanics, the standard gravitational parameter was often used instead of the individual masses. Also, since the mass of the Sun is so much larger than the mass of any planet, as long as they had a good estimate of the mass of the Sun, they could calculate the orbits of the planets without knowing their masses.
So we need to solve only this nonlinear second order ODE to get the relative motion of the two bodies. The individual motions can be recovered from the relative motion.
Angular Momentum Conservation
We will proceed as the previous section, and calculate all the integrals of motion.
Taking a cross product of the equation with \(\mathbf{r}\), we have
\[ \mathbf{r} \times \ddot{\mathbf{r}} = 0 \]
Since, \(\frac{d}{dt}(\mathbf{r} \times \dot{\mathbf{r}}) = \mathbf{r} \times \ddot{\mathbf{r}}\), we have
\[ \frac{d}{dt}(\mathbf{r} \times \dot{\mathbf{r}}) = 0 \]
Defining \(\mathbf{h} = \mathbf{r} \times \dot{\mathbf{r}}\), we have \(\frac{d\mathbf{h}}{dt} = 0\). Therefore, \(\mathbf{h}\) is a constant vector. This will only happen if the motion is planar. Therefore, the motion of the two bodies is planar.
This is a strong geometric statement: the trajectory in 3D collapses to a curve in a single 2D plane. From here on, we can do orbital geometry in two dimensions.
Now let us recall Kepler’s second law, which states that the line joining the two bodies sweeps out equal areas in equal times. This is a consequence of the conservation of angular momentum. Since \(\mathbf{h}\) is constant, the areal velocity is also constant.
Eccentricity Vector
So solve the two-body problem, let us now post multiply the equation with \(\mathbf{h}\) to get
\[ \ddot{\mathbf{r}} \times \mathbf{h} = -\frac{\mu}{r^3} \mathbf{r} \times \mathbf{h} = -\frac{\mu}{r^3} \mathbf{r} \times \mathbf{r} \times \dot{\mathbf{r}} \]
Since \(\mathbf{a} \times (\mathbf{b} \times \mathbf{c}) = (\mathbf{a} \cdot \mathbf{c})\mathbf{b} - (\mathbf{a} \cdot \mathbf{b})\mathbf{c}\), we have
\[ \ddot{\mathbf{r}} \times \mathbf{h} = -\frac{\mu}{r^3} \left[ (\mathbf{r} \cdot \dot{\mathbf{r}})\mathbf{r} - (\mathbf{r} \cdot \mathbf{r})\dot{\mathbf{r}} \right] = -\frac{\mu}{r^3} \left[ (\mathbf{r} \cdot \dot{\mathbf{r}})\mathbf{r} - r^2\dot{\mathbf{r}} \right] \]
Since \(\mathbf{r} = r \hat{\mathbf{r}}\), we have \(\dot{\mathbf{r}} = \dot{r} \hat{\mathbf{r}} + r \dot{\hat{\mathbf{r}}}\), and hence
\[ \begin{aligned} \mathbf{r} \cdot \dot{\mathbf{r}} & = & r \hat{\mathbf{r}} \cdot (\dot{r} \hat{\mathbf{r}} + r \dot{\hat{\mathbf{r}}})\\ & = & r \dot{r} \hat{\mathbf{r}} \cdot \hat{\mathbf{r}} + r^2 \hat{\mathbf{r}} \cdot \dot{\hat{\mathbf{r}}}\\ & = & r \dot{r} + 0\\ \end{aligned} \]
We have used, \(\frac{d}{dt}(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}) = 0\), which implies \(\hat{\mathbf{r}} \cdot \dot{\hat{\mathbf{r}}} = 0\).
Therefore,
\[ \ddot{\mathbf{r}} \times \mathbf{h} = -\frac{\mu}{r^3} \left[ r \dot{r} \mathbf{r} - r^2\dot{\mathbf{r}} \right] = \mu \left( \frac{\dot{\mathbf{r}}}{r} - \frac{\dot{r} \mathbf{r}}{r^2} \right) \]
Since, \[ \frac{d}{dt}\left( \frac{\mathbf{r}}{r} \right) = \frac{\dot{\mathbf{r}}}{r} - \frac{\dot{r}}{r^2}\mathbf{r} \]
we have,
\[ \ddot{\mathbf{r}} \times \mathbf{h} = \mu \frac{d}{dt}\left( \frac{\mathbf{r}}{r} \right) \]
Since \(\mathbf{h}\) is constant, integrating both sides gives
\[ \dot{\mathbf{r}} \times \mathbf{h} = \mu \left( \frac{\mathbf{r}}{r} + \mathbf{e} \right) \]
Here, \(\mathbf{e}\) is a dimensionless constant vector. We note that it is perpendicular to both the \(\mathbf{h}\) and \(\dot{\mathbf{r}}\) vectors. Since \(\mathbf{h}\) is a constant vector, \(\mathbf{e}\) clearly lies in the orbital plane.
Solving the last equation for \(\mathbf{e}\),
\[ \mathbf{e} = \frac{\dot{\mathbf{r}} \times \mathbf{h}}{\mu} - \frac{\mathbf{r}}{r} \]
Differentiating with respect to time, and using \(\dot{\mathbf{h}} = 0\),
\[ \dot{\mathbf{e}} = \frac{\ddot{\mathbf{r}} \times \mathbf{h}}{\mu} - \frac{d}{dt}\left( \frac{\mathbf{r}}{r} \right) \]
But we have already shown that \(\ddot{\mathbf{r}} \times \mathbf{h} = \mu \frac{d}{dt}\left( \frac{\mathbf{r}}{r} \right)\), and hence the two terms cancel exactly,
\[ \dot{\mathbf{e}} = \mathbf{0} \]
A vector is fixed only if both its magnitude and its direction are fixed, so let us extract the two statements separately.
Magnitude. Writing \(e = |\mathbf{e}|\),
\[ \frac{d}{dt}\left( e^2 \right) = \frac{d}{dt}(\mathbf{e} \cdot \mathbf{e}) = 2\, \mathbf{e} \cdot \dot{\mathbf{e}} = 0 \quad \Longrightarrow \quad e = \text{constant} \]
Direction. With \(\hat{\mathbf{e}} = \mathbf{e}/e\), and using \(\dot{e} = 0\) from above,
\[ \dot{\hat{\mathbf{e}}} = \frac{\dot{\mathbf{e}}}{e} - \frac{\dot{e}}{e^2}\mathbf{e} = \mathbf{0} - \mathbf{0} = \mathbf{0} \quad \Longrightarrow \quad \hat{\mathbf{e}} = \text{constant} \]
So \(\mathbf{e}\) picks out one fixed line in the orbital plane, for all time. This is a much stronger statement than the conservation of \(\mathbf{h}\) alone: \(\mathbf{h}\) fixes the plane of motion, while \(\mathbf{e}\) additionally fixes an orientation within that plane. We shall soon see that this fixed direction points towards the point of closest approach (the periapsis), and that \(e\) is precisely the eccentricity of the conic. Its constancy is therefore the statement that the orbit does not precess — the ellipse traced out is closed, and returns on itself.
Note
The vector \(\mu\,\mathbf{e}\) is known as the Laplace–Runge–Lenz vector, and its conservation is special to the inverse-square law. Any perturbation to \(1/r^2\) — oblateness of the central body, a third body, or the corrections of general relativity — makes \(\dot{\mathbf{e}} \neq \mathbf{0}\), and the apsidal line then slowly rotates. The residual perihelion precession of Mercury is the most famous instance of exactly this. Under certain conditions, it may change the shape of the orbit itself.
Two body motion
Figure 1: Two-body motion at two mass ratios. Both panels share the same total mass, the same initial separation and the same eccentricity (\(e = 0.70\)); only \(m_1/m_2\) differs. Green arrows are velocities (drawn to scale), purple arrows the gravitational forces (drawn to a compressed \(F^{1/3}\) scale, since \(F \propto 1/d^2\) spans too wide a range).
What the animation shows
The centre of mass never moves. Taking the origin at the barycentre, the two bodies stay diametrically opposite about it: \(m_1 \mathbf{R}_1 + m_2 \mathbf{R}_2 = \mathbf{0}\) at every instant.
With \(\mathbf{r} = \mathbf{r}_{12} = \mathbf{R}_2 - \mathbf{R}_1\) and \(M = m_1 + m_2\) as before, the relative coordinate obeys \(\ddot{\mathbf{r}} = -GM\,\mathbf{r}/r^3\), and the individual orbits are recovered from \[\mathbf{R}_1 = -\frac{m_2}{M}\,\mathbf{r}, \qquad \mathbf{R}_2 = \frac{m_1}{M}\,\mathbf{r}.\] So each body traces a similar conic about the centre of mass: same shape, same eccentricity, sizes in the inverse ratio of the masses. Note the weights are crossed — body 1 is scaled by \(m_2\) — so the heavier body gets the smaller orbit.
How to find \(R_1\) and \(R_2\)
Define the change of variables \((\mathbf{R}_1, \mathbf{R}_2) \mapsto (\mathbf{R}_{cm}, \mathbf{r})\):
\[ \begin{aligned} \mathbf{R}_{cm} &= \frac{m_1}{M}\mathbf{R}_1 + \frac{m_2}{M}\mathbf{R}_2, \\ \mathbf{r} &= -\mathbf{R}_1 + \mathbf{R}_2 . \end{aligned} \]
In matrix form, acting componentwise on each Cartesian component,
\[ \begin{pmatrix} \mathbf{R}_{cm} \\ \mathbf{r} \end{pmatrix} = \underbrace{\begin{pmatrix} m_1/M & m_2/M \\ -1 & 1 \end{pmatrix}}_{\displaystyle A} \begin{pmatrix} \mathbf{R}_1 \\ \mathbf{R}_2 \end{pmatrix}. \]
- It is invertible as \[ \det A = \frac{m_1}{M} + \frac{m_2}{M} = 1 \neq 0 . \]
Invert \(A\) directly:
\[ A^{-1} = \frac{1}{\det A}\begin{pmatrix} 1 & -m_2/M \\ 1 & \phantom{-}m_1/M\end{pmatrix} = \begin{pmatrix} 1 & -m_2/M \\ 1 & \phantom{-}m_1/M\end{pmatrix}, \]
that is,
\[ \boxed{\; \begin{aligned} \mathbf{R}_1(t) &= \mathbf{R}_{cm}(t) - \frac{m_2}{M}\,\mathbf{r}(t), \\[4pt] \mathbf{R}_2(t) &= \mathbf{R}_{cm}(t) + \frac{m_1}{M}\,\mathbf{r}(t). \end{aligned}\;} \]
- The two forces are equal and opposite (Newton’s third law), but the accelerations are not: the heavier body responds \(m_1/m_2\) times less.
Towards the restricted problem
As \(m_1/m_2\) grows, the orbit of the heavier body shrinks and the centre of mass sinks towards its centre — this is the approximation of the previous section, now visible.
For the Sun-Jupiter pair, \(m_1/m_2 \approx 1047\), so the Sun’s orbit about the barycentre is \(1/1047\) of Jupiter’s. It is small, but not zero: the barycentre lies about \(1.07\) solar radii from the centre of the Sun, i.e. just outside its surface.
A deliberately three-dimensional start
Nothing so far has assumed a plane. So let us start two equal bodies in a way that has nothing planar about it, and let the integrals do the work. Take \(G = 1\), so that \(\mu = G(m_1 + m_2) = 1\):
\[ \begin{aligned} m_1 = m_2 &= \tfrac{1}{2}, \\[4pt] \mathbf{R}_1(0) &= \left(+\tfrac{1}{2},\, 0,\, 0\right), &\qquad \mathbf{V}_1(0) &= \tfrac{1}{2}\,\hat{\mathbf{y}}, \\ \mathbf{R}_2(0) &= \left(-\tfrac{1}{2},\, 0,\, 0\right), &\qquad \mathbf{V}_2(0) &= \tfrac{1}{2}\,\hat{\mathbf{z}} . \end{aligned} \]
The bodies sit on the \(x\)-axis, equidistant from the origin; one is launched along \(y\), the other along \(z\). The separation lies along one axis and the two velocities along the other two — all three directions mutually perpendicular. Worse, the total linear momentum
\[ \mathbf{P} = m_1\mathbf{V}_1 + m_2\mathbf{V}_2 = \tfrac{1}{4}\left(\hat{\mathbf{y}} + \hat{\mathbf{z}}\right) \neq \mathbf{0} \]
does not vanish, so the barycentre drifts. This is not the tidy equal-and-opposite setup of the previous animation.
The orbit, before any integration
The relative state at \(t = 0\) is
\[ \mathbf{r}(0) = \mathbf{R}_2 - \mathbf{R}_1 = -\hat{\mathbf{x}}, \qquad \dot{\mathbf{r}}(0) = \mathbf{V}_2 - \mathbf{V}_1 = \tfrac{1}{2}\left(\hat{\mathbf{z}} - \hat{\mathbf{y}}\right), \]
so the integrals are fixed once and for all by these two vectors:
\[ \mathbf{h} = \mathbf{r} \times \dot{\mathbf{r}} = \tfrac{1}{2}\left(\hat{\mathbf{y}} + \hat{\mathbf{z}}\right), \qquad h = \tfrac{1}{\sqrt{2}}, \qquad \varepsilon = \frac{v^{2}}{2} - \frac{\mu}{r} = \frac{1}{4} - 1 = -\frac{3}{4}, \]
and — borrowing three results we shall derive in the next few slides (vis-viva, the period, and the second Casimir relation) — the whole geometry follows without solving anything:
\[ a = -\frac{\mu}{2\varepsilon} = \frac{2}{3}, \qquad e = \sqrt{1 + \frac{2\varepsilon h^{2}}{\mu^{2}}} = \frac{1}{2}, \qquad T = 2\pi\sqrt{\frac{a^{3}}{\mu}} = 3.4201 . \]
Since \(\mathbf{r}\cdot\dot{\mathbf{r}} = 0\) at \(t = 0\), the run starts at an apsis; and \(r(0) = 1 = a(1+e)\) identifies it as apoapsis.
Helices, not closed curves
One more feature of this choice of initial data, which the animation is about to make obvious. The barycentre travels at
\[ \mathbf{V}_{cm} = \frac{\mathbf{P}}{M} = \tfrac{1}{4}\left(\hat{\mathbf{y}} + \hat{\mathbf{z}}\right) \qquad\text{while}\qquad \mathbf{h} = \tfrac{1}{2}\left(\hat{\mathbf{y}} + \hat{\mathbf{z}}\right), \]
so \(\mathbf{V}_{cm} \parallel \mathbf{h}\): the drift is along the normal of the orbital plane. The plane therefore slides through space like a sheet pushed face-first, and neither body ever returns to a point it has already visited — each traces a helix, even though the relative orbit is a closed ellipse of period \(T = 3.4201\).
Two bodies in 3D
Figure 2: Left: the inertial frame — the barycentre (amber) runs off along a straight line, dragging the orbital plane with it, and each body traces a helix. Right: the barycentric frame — the drift subtracted, the two bodies run on similar (here congruent, as \(m_1 = m_2\)) ellipses in one plane that never moves. Shaded patch: the plane \(\perp \mathbf{h}\). Integrated with adaptive RKF7(8), \(rtol = 10^{-12}\).
Why the motion is planar
1. \(\mathbf{h}\) is a fixed vector. From the equation of motion, \(\dot{\mathbf{h}} = \mathbf{r} \times \ddot{\mathbf{r}} = -\frac{\mu}{r^{3}}\,\mathbf{r} \times \mathbf{r} = \mathbf{0}\).
2. \(\mathbf{r}\) never leaves the plane \(\perp \mathbf{h}\). At every instant
\[ \mathbf{r} \cdot \mathbf{h} = \mathbf{r} \cdot \left(\mathbf{r} \times \dot{\mathbf{r}}\right) = 0, \qquad \dot{\mathbf{r}} \cdot \mathbf{h} = \dot{\mathbf{r}} \cdot \left(\mathbf{r} \times \dot{\mathbf{r}}\right) = 0, \]
both being triple products with a repeated vector. Position and velocity lie in the plane through the origin with normal \(\mathbf{h}\).
3. It is the same plane for all time. Step 2 on its own would only give a plane at each instant; it is step 1 — \(\mathbf{h}\) fixed in both magnitude and direction — that stops the normal from rotating. One plane, fixed for all \(t\): the invariable plane of the pair.
Why the motion is planar (contd.)
4. From \(\mathbf{r}\) to the two bodies. Using the inversion derived earlier,
\[ \mathbf{R}_1 = \mathbf{R}_{cm} - \frac{m_2}{M}\,\mathbf{r}, \qquad \mathbf{R}_2 = \mathbf{R}_{cm} + \frac{m_1}{M}\,\mathbf{r} . \]
In the barycentric frame \(\mathbf{R}_{cm} \equiv \mathbf{0}\), so both bodies are scalar multiples of the same vector \(\mathbf{r}\) — and therefore lie in that one fixed plane through the barycentre. This is the right-hand panel of Figure 2.
5. In the inertial frame. The linear-momentum integral gives \(\mathbf{R}_{cm}(t) = \mathbf{R}_{cm}(0) + \mathbf{V}_{cm}t\), so the plane translates rigidly, always parallel to itself, its normal \(\hat{\mathbf{h}}\) never changing. Here \(\mathbf{V}_{cm} \parallel \mathbf{h}\), so it slides along its own normal and the trajectories are helices (left panel).
The one exception
If \(\mathbf{h} = \mathbf{0}\) — that is, \(\mathbf{r} \parallel \dot{\mathbf{r}}\), as when the pair is released from rest — there is no preferred normal and the argument says nothing. Correctly so: the motion is then rectilinear, the two bodies fall straight at each other, and every plane containing that line will do.
Does the integration agree?
The run of Figure 2 uses an adaptive Runge–Kutta–Fehlberg 7(8) (codes/two-body-3d-plane.py --verify), which knows nothing about any of the integrals above — it only advances \(\ddot{\mathbf{R}}_i\). Over two full periods, with \(158\) accepted steps:
| quantity | meaning | value |
|---|---|---|
| \(\max_t\left|\left(\mathbf{R}_i - \mathbf{R}_{cm}\right)\cdot\hat{\mathbf{h}}\right|\) | departure from the plane | \(1.7 \times 10^{-16}\) |
| \(\Delta h\) | drift in angular momentum | \(2.4 \times 10^{-13}\) |
| \(\Delta \varepsilon\) | drift in energy | \(5.9 \times 10^{-13}\) |
| \(\Delta e\) | drift in eccentricity | \(1.1 \times 10^{-12}\) |
| \(\max_t\left|\mathbf{R}_{cm} - \left(\mathbf{R}_{cm}(0) + \mathbf{V}_{cm}t\right)\right|\) | barycentre drift is uniform | \(1.1 \times 10^{-15}\) |
The planarity is exact to rounding, not merely to the integration tolerance: it is an algebraic identity, not a numerical accident.
What must be true of the initial data for both bodies to move in a single fixed plane in the inertial frame as well?
Energy Integral (Vis-Viva)
Important
We shall see, dotting the EOM with \(\dot{\mathbf{r}}\) produces a scalar conservation law — the energy integral, also known as the vis-viva (literally, “living force”) equation.
This is arguably the single most useful formula in mission design.
Dot-multiply the equation of motion by \(\dot{\mathbf{r}}\):
\[ \dot{\mathbf{r}} \cdot \ddot{\mathbf{r}} = -\frac{\mu}{r^3}\,\mathbf{r}\cdot\dot{\mathbf{r}}, \]
\[ \frac{d}{dt}\!\left(\frac{v^2}{2}\right) = -\frac{\mu\dot{r}}{r^2} = \frac{d}{dt}\!\left(\frac{\mu}{r}\right). \]
Integrating
\[ \frac{v^2}{2} - \frac{\mu}{r} = \text{const}, \]
This is the first form of the vis-viva equation.
Equation of the Orbit
We now want the explicit geometric shape of the orbit — \(r\) as a function of position angle. The cleanest route uses what we have already derived. We have \(\mathbf{h}\) (which fixes the orbital plane) and \(\mathbf{e}\) (which fixes the line of apsides). All we have to do is exploit a single dot-product identity to read off \(r(f)\).
Take the dot product of \(\mathbf{r}\) with the eccentricity-vector integral \(\dot{\mathbf{r}} \times \mathbf{h} = \mu(\hat{\mathbf{r}} + \mathbf{e})\):
\[ \mathbf{r} \cdot (\dot{\mathbf{r}} \times \mathbf{h}) = (\mathbf{r} \times \dot{\mathbf{r}}) \cdot \mathbf{h} = \mathbf{h} \cdot \mathbf{h} = h^2. \]
On the right side, let \(f\) be the angle between \(\mathbf{e}\) and \(\mathbf{r}\) (the true anomaly):
\[ \mathbf{r} \cdot \mu(\hat{\mathbf{r}} + \mathbf{e}) = \mu(r + re\cos f). \]
Equating and solving for \(r\) (Battin Eq. 3.20):
\[ \boxed{r = \frac{p}{1 + e\cos f}}, \qquad p = \frac{h^2}{\mu}. \]
This is the polar equation of a conic section with focus at the origin. The single number \(e\) — the magnitude of the eccentricity vector — determines the qualitative shape of the orbit:
| \(e\) | Conic | Orbit type |
|---|---|---|
| \(0\) | circle | bound |
| \(0 < e < 1\) | ellipse | bound |
| \(e = 1\) | parabola | escape, zero energy |
| \(e > 1\) | hyperbola | escape, positive energy |
The ellipse
Figure 3: Geometry of the elliptic orbit (\(e = 0.6\)). The attracting body sits at the occupied focus \(F\); the other focus is empty. The centre-to-focus distance is \(ae\), so the eccentricity is the fractional offset of the focus from the centre. The semi-latus rectum \(p\) is the radius at \(f = \pi/2\), i.e. \(p = a(1-e^{2}) = b^{2}/a\), and the true anomaly \(f\) is measured at the focus from the periapsis direction.
Various Conic sections
All four numbers below describe the same orbit — move any one of them and watch the other three, and the curve, follow. Units are non-dimensional, with \(\mu = G(m_1+m_2) = 1\); the attracting body sits at the amber focus and periapsis is to the right.
Things to try
Hold \(h\) and sweep \(e\) from \(0\) to \(2.5\). Every curve in the family keeps the same semi-latus rectum \(p = h^{2}/\mu\) — the conic simply unrolls from circle to ellipse to parabola to hyperbola, and \(\varepsilon\) changes sign exactly at \(e = 1\).
Hold \(a\) (or, identically, \(\varepsilon\)) and sweep \(e\). Now the size is frozen: every ellipse has the same major axis, the same energy and the same period \(T\), but the focus slides from the centre out towards the periapsis. Energy is blind to shape.
Notice that the \(a\) and \(\varepsilon\) buttons do the same thing. They must: \(\varepsilon = -\mu/2a\) is a one-to-one map, so \((a, e)\) and \((\varepsilon, e)\) are the same parametrisation, and only \((h, e)\) is genuinely different.
Try to cross \(e = 1\) while holding \(a\). You cannot — the slider stops. An ellipse cannot become a hyperbola without \(a\) passing through infinity, which is the Casimir relation \(\mu^{2}\left(e^{2}-1\right) = 2\varepsilon h^{2}\) forbidding \(\operatorname{sign}\varepsilon \neq \operatorname{sign}(e-1)\). Hold \(h\) instead and the passage is smooth.
Two spacecraft are on orbits with the same \(a\) but \(e = 0.1\) and \(e = 0.9\). Which one is moving faster at periapsis, and which has the larger \(h\)?
Polar form of the conic in \(a\) and \(e\)
The orbit equation \(r = p/(1 + e\cos f)\) is complete, but \(p\) is not the number anyone quotes for a mission — the semi-major axis \(a\) is. So let us trade \(p\) for \(a\).
For a bound orbit (\(e < 1\)) evaluate \(r\) at the two apsides, where the radius is stationary:
\[ f = 0: \quad r_p = \frac{p}{1 + e}, \qquad\qquad f = \pi: \quad r_a = \frac{p}{1 - e}. \]
The major axis runs from apoapsis to periapsis straight through the focus, so \(2a = r_p + r_a\):
\[ 2a = p\left(\frac{1}{1+e} + \frac{1}{1-e}\right) = p\,\frac{(1-e) + (1+e)}{1 - e^{2}} = \frac{2p}{1 - e^{2}}, \]
that is,
\[ \boxed{\;p = a\left(1 - e^{2}\right)\;} \qquad\Longrightarrow\qquad \boxed{\;r = \frac{a\left(1 - e^{2}\right)}{1 + e\cos f}\;} \]
This is the polar form of the conic in the variables we shall actually use. Three corollaries we need at once:
\[ r_p = a(1 - e), \qquad r_a = a(1 + e), \qquad h^{2} = \mu p = \mu a\left(1 - e^{2}\right). \]
Vis-viva Eqn
We left the energy integral as \(\dfrac{v^{2}}{2} - \dfrac{\mu}{r} = \text{const}\) without saying what the constant is. Call it \(\varepsilon\) — the specific orbital energy, energy per unit reduced mass — and let us now evaluate it.
Since \(\varepsilon\) is constant, we may evaluate it anywhere on the orbit. Periapsis is the cheapest point: there \(\dot r = 0\), the velocity is purely transverse, and hence \(h = r_p v_p\), i.e. \(v_p = h/r_p\):
\[ \varepsilon = \frac{v_p^{2}}{2} - \frac{\mu}{r_p} = \frac{h^{2}}{2 r_p^{2}} - \frac{\mu}{r_p}. \]
Substituting \(h^{2} = \mu a\left(1 - e^{2}\right)\) and \(r_p = a(1-e)\) from the previous slide:
\[ \begin{aligned} \varepsilon &= \frac{\mu a\left(1 - e^{2}\right)}{2 a^{2}(1-e)^{2}} - \frac{\mu}{a(1-e)} = \frac{\mu (1+e)}{2 a (1-e)} - \frac{\mu}{a(1-e)} \\[4pt] &= \frac{\mu}{a(1-e)}\left[\frac{1+e}{2} - 1\right] = \frac{\mu}{a(1-e)} \cdot \frac{e-1}{2} = -\frac{\mu}{2a}. \end{aligned} \]
The constant is therefore fixed by the size of the orbit alone:
\[ \boxed{\;\frac{v^{2}}{2} - \frac{\mu}{r} = -\frac{\mu}{2a}\;} \qquad\Longleftrightarrow\qquad \boxed{\;v^{2} = \mu\left(\frac{2}{r} - \frac{1}{a}\right)\;} \]
This is the vis-viva equation in its working form.
Eccentricity has vanished. Two orbits with the same \(a\) but wildly different shapes carry exactly the same energy — and, as we have just seen, the same period \(T = 2\pi\sqrt{a^{3}/\mu}\). Energy is blind to shape.
The sign of \(\varepsilon\) classifies the conic, matching the table of \(e\) we drew earlier:
| Conic | \(a\) | \(\varepsilon = -\mu/2a\) |
|---|---|---|
| ellipse | \(a > 0\) | \(\varepsilon < 0\) (bound) |
| parabola | \(a \to \infty\) | \(\varepsilon = 0\) (marginal) |
| hyperbola | \(a < 0\) | \(\varepsilon > 0\) (unbound) |
Two speeds fall out immediately. Put \(a = r\) (a circle) to get the circular speed \(v_c = \sqrt{\mu/r}\); put \(\varepsilon = 0\) (\(a \to \infty\)) to get the escape speed \(v_{esc} = \sqrt{2\mu/r} = \sqrt{2}\, v_c\).
Every \(\Delta v\) budget in mission design starts from this one equation: given where you are (\(r\)) and the orbit you want (\(a\)), it hands you the speed you need.
Period (\(T\)) of an orbit
Step 1 — the areal velocity. With \(\mathbf{r} = r\hat{\mathbf{r}}\) and \(\dot{\mathbf{r}} = \dot r\,\hat{\mathbf{r}} + r\,\dot{\hat{\mathbf{r}}}\) (used earlier in the eccentricity-vector derivation),
\[ \mathbf{h} = \mathbf{r} \times \dot{\mathbf{r}} = r\hat{\mathbf{r}} \times \left(\dot r\,\hat{\mathbf{r}} + r\,\dot{\hat{\mathbf{r}}}\right) = r^{2}\,\hat{\mathbf{r}} \times \dot{\hat{\mathbf{r}}} . \]
Because \(\mathbf{e}\) is a fixed vector, the apsidal line does not move, so the true anomaly \(f\) is the polar angle of \(\hat{\mathbf{r}}\): writing \(\dot{\hat{\mathbf{r}}} = \dot f\,\hat{\boldsymbol{\theta}}\) with \(\hat{\boldsymbol{\theta}} \perp \hat{\mathbf{r}}\) in the orbital plane,
\[ h = r^{2}\dot f . \]
The radius vector sweeps a thin sector of area \(dA = \tfrac{1}{2} r^{2}\,df\), hence
\[ \frac{dA}{dt} = \frac{1}{2} r^{2}\dot f = \frac{h}{2} = \text{const}, \]
which is Kepler’s second law once more. Constant rate is the whole point: the area is swept uniformly, so one revolution can be timed by simple division.
Step 2 — the area enclosed. For an ellipse \(A = \pi a b\), so we need \(b\); and \(b\) follows from the orbit equation alone. Put \(\cos f = -e\) in \(r = a(1-e^{2})/(1 + e\cos f)\):
\[ r = \frac{a\left(1 - e^{2}\right)}{1 - e^{2}} = a, \qquad r\cos f = -ae . \]
That point lies a distance \(ae\) to the left of the focus; but the centre is itself \(ae\) to the left of the focus (Figure 3), so the point sits directly above the centre — it is the end of the minor axis. By Pythagoras,
\[ b^{2} = a^{2} - (ae)^{2} \qquad \Longrightarrow \qquad b = a\sqrt{1 - e^{2}} = \sqrt{a p} . \]
Step 3 — divide. One revolution sweeps the whole ellipse at the constant rate \(h/2\):
\[ T = \frac{A}{dA/dt} = \frac{\pi a b}{h/2} = \frac{2\pi\, a \cdot a\sqrt{1 - e^{2}}}{\sqrt{\mu a \left(1 - e^{2}\right)}} = \frac{2\pi a^{2}}{\sqrt{\mu a}}, \]
\[ \boxed{\;T = 2\pi\sqrt{\frac{a^{3}}{\mu}}\;} \qquad\qquad n \equiv \frac{2\pi}{T} = \sqrt{\frac{\mu}{a^{3}}} \]
where \(n\) is the mean motion, the average angular rate — the quantity we shall need when we come to Kepler’s equation.
Note once again that \(e\) has cancelled: the period, like the energy, depends only on \(a\). Indeed \(T\) and \(\varepsilon = -\mu/2a\) are two readings of the same number. And \(T\) exists only for \(a > 0\): a parabola or hyperbola never closes.
Kepler’s third law Modification
Squaring, \(T^{2} = \dfrac{4\pi^{2}}{G(m_1 + m_2)}\,a^{3}\). Kepler asserted \(T^{2} \propto a^{3}\) with the same constant for every planet; Newton’s version shows the constant carries \(m_1 + m_2\). For the Sun and even its heaviest planet \(m_2/m_1 \approx 10^{-3}\), so the discrepancy is a few parts in \(10^{4}\) — comfortably below Kepler’s observational precision, which is why he never saw it.
Casimir relations
We have now squeezed the two-body equation of motion three different ways — cross product, BAC–CAB, dot product — and harvested three integrals:
\[ \mathbf{h} = \mathbf{r} \times \dot{\mathbf{r}}, \qquad \mathbf{e} = \frac{\dot{\mathbf{r}} \times \mathbf{h}}{\mu} - \frac{\mathbf{r}}{r}, \qquad \varepsilon = \frac{v^2}{2} - \frac{\mu}{r} \]
That is \(3 + 3 + 1 = 7\) scalar constants for a system whose state has only \(6\) degrees of freedom. They cannot all be independent. Exactly two relations tie them together, and these are called the Casimir relations:
\[ \boxed{\;\mathbf{h} \cdot \mathbf{e} = 0\;} \qquad\qquad \boxed{\;\mu^{2}\left(e^{2} - 1\right) = 2\,\varepsilon\, h^{2}\;} \]
Neither is a new integral. They are identities — constraints satisfied by the integrals we already have, at every instant, on every orbit.
First relationship is trivial. So we will prove the second.
Square the expression for \(\mathbf{e}\),
\[ \mu^{2} e^{2} = \underbrace{\left| \dot{\mathbf{r}} \times \mathbf{h} \right|^{2}}_{\text{(i)}} - 2\mu\, \underbrace{\hat{\mathbf{r}} \cdot \left( \dot{\mathbf{r}} \times \mathbf{h} \right)}_{\text{(ii)}} + \mu^{2}\, \underbrace{\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}}_{= \,1} \]
Term (i): \(\mathbf{h} = \mathbf{r} \times \dot{\mathbf{r}}\) is perpendicular to \(\dot{\mathbf{r}}\), so the cross product has magnitude \(v h\) with no \(\sin\) factor to worry about:
\[ \left| \dot{\mathbf{r}} \times \mathbf{h} \right|^{2} = v^{2} h^{2} \]
Term (ii): cycle the triple product so that \(\mathbf{h}\) meets \(\mathbf{r} \times \dot{\mathbf{r}}\):
\[ \hat{\mathbf{r}} \cdot \left( \dot{\mathbf{r}} \times \mathbf{h} \right) = \frac{1}{r}\, \mathbf{r} \cdot \left( \dot{\mathbf{r}} \times \mathbf{h} \right) = \frac{1}{r}\, \mathbf{h} \cdot \left( \mathbf{r} \times \dot{\mathbf{r}} \right) = \frac{h^{2}}{r} \]
Substituting both terms,
\[ \mu^{2} e^{2} = v^{2} h^{2} - \frac{2 \mu h^{2}}{r} + \mu^{2} = h^{2} \left( v^{2} - \frac{2\mu}{r} \right) + \mu^{2} \]
The bracket is precisely twice the specific energy, \(v^{2} - 2\mu/r = 2\varepsilon\). Hence
\[ \mu^{2} e^{2} = 2 \varepsilon h^{2} + \mu^{2} \qquad \Longrightarrow \qquad \mu^{2}\left( e^{2} - 1 \right) = 2 \varepsilon h^{2} \]
or, solved for the eccentricity,
\[ e = \sqrt{1 + \frac{2 \varepsilon h^{2}}{\mu^{2}}} \]
Since \(h^{2} \ge 0\) and \(\mu > 0\), the sign of \(\varepsilon\) alone decides the shape of the conic:
| Energy | Eccentricity | Orbit |
|---|---|---|
| \(\varepsilon < 0\) | \(e < 1\) | ellipse (bound) |
| \(\varepsilon = 0\) | \(e = 1\) | parabola (escape) |
| \(\varepsilon > 0\) | \(e > 1\) | hyperbola (unbound) |
This is why the vis-viva constant is the natural “mission design” variable: it tells you at a glance whether you are captured or escaping.
Substituting the vis-viva result \(\varepsilon = -\mu / (2a)\) into relation 2 gives
\[ \mu^{2}\left(e^{2} - 1\right) = -\frac{\mu h^{2}}{a} \qquad \Longrightarrow \qquad h^{2} = \mu\, a \left( 1 - e^{2} \right) = \mu\, p \]
where \(p = a(1 - e^{2})\) is the semi-latus rectum. So the orbit equation \(r = p / (1 + e\cos\theta)\) can be written with \(p = h^{2}/\mu\) — the geometry of the conic is fixed entirely by the two dynamical constants \(h\) and \(\varepsilon\).
Counting the integrals
The two Casimir relations settle the bookkeeping we started in the N-body section:
- We possess \(7\) scalar constants: \(\mathbf{h}\,(3)\), \(\mathbf{e}\,(3)\), \(\varepsilon\,(1)\).
- Relation 1 removes one; relation 2 removes another.
- That leaves \(5\) independent constants of the motion.
Five constants fix the orbit as a curve in space — its plane (\(i, \Omega\)), its shape (\(a, e\)) and its orientation in the plane (\(\omega\)). What they cannot tell us is where on that curve the body is. One more constant, the time of periapsis passage \(t_p\), supplies that.
\[ 5 \;(\text{geometry}) \;+\; 1 \;(\text{timing, } t_p) \;=\; 6 \;=\; \text{order of the reduced system} \]
and these six are exactly the classical orbital elements \((a, e, i, \Omega, \omega, t_p)\) that we will use for the rest of the course.
The N-body problem had only \(10\) integrals for a \(6N\)-order system. Why does the two-body problem manage to close the gap completely?
Why the name “Casimir”?
Optional aside — the hidden symmetry of the Kepler problem
Treat \(\mathbf{h}\) and the scaled eccentricity vector \(\mathbf{A} = \mu \mathbf{e}\) as generators and compute their Poisson brackets (unit reduced mass):
\[ \{h_i, h_j\} = \epsilon_{ijk} h_k, \qquad \{h_i, A_j\} = \epsilon_{ijk} A_k, \qquad \{A_i, A_j\} = -2\varepsilon\, \epsilon_{ijk} h_k \]
The first two say \(\mathbf{h}\) generates rotations and \(\mathbf{A}\) transforms as a vector under them — expected. The third is the surprise: the components of \(\mathbf{A}\) do not commute, and they close back onto \(\mathbf{h}\). For a bound orbit (\(\varepsilon < 0\)) put \(\mathbf{D} = \mathbf{A} / \sqrt{-2\varepsilon}\), and the algebra becomes \(\{D_i, D_j\} = \epsilon_{ijk} h_k\) — the Lie algebra \(\mathfrak{so}(4)\), the rotation group of four-dimensional space. (For \(\varepsilon > 0\) one gets \(\mathfrak{so}(3,1)\); for \(\varepsilon = 0\), the Euclidean algebra \(\mathfrak{e}(3)\).)
A Casimir invariant is a function of the generators that Poisson-commutes with all of them. The algebra \(\mathfrak{so}(4)\) has exactly two, and evaluated on the Kepler problem they are
\[ \mathbf{h} \cdot \mathbf{D} = 0, \qquad h^{2} + D^{2} = -\frac{\mu^{2}}{2\varepsilon} \]
which are relations 1 and 2 verbatim. So the two identities are not accidents of algebra — they are the statement that the inverse-square law has a symmetry larger than the obvious rotational one. This extra \(\mathfrak{so}(4)\) symmetry is the reason Kepler orbits close after exactly one revolution, and the reason the hydrogen atom’s energy levels are degenerate in \(\ell\). It is also why the \(1/r\) potential is special: perturb the exponent even slightly and \(\mathbf{e}\) stops being constant — the orbit precesses.
Solution of orbital equation
SFM, IIST 2026