Assignment 02 — Solutions

Companion to Assignment 02. Notation as in the lecture slides.

Part A — The N-body Problem

Solution 1 — The tenth integral

(a) Collapsing the double sum. Write the force on \(i\) due to \(j\) as \[ \mathbf{F}_{ij} = G\frac{m_i m_j}{r_{ij}^{2}}\hat{\mathbf{r}}_{ij} = G m_i m_j \frac{\mathbf{r}_{ij}}{r_{ij}^{3}}, \qquad \mathbf{r}_{ij} = \mathbf{R}_j - \mathbf{R}_i . \] Dotting the equation of motion with \(\dot{\mathbf{R}}_i\) and summing over \(i\), \[ \underbrace{\sum_i m_i \ddot{\mathbf{R}}_i\cdot\dot{\mathbf{R}}_i}_{\text{LHS}} = \sum_i \sum_{j\neq i} \mathbf{F}_{ij}\cdot\dot{\mathbf{R}}_i . \] The left side is already exact: \[ \text{LHS} = \sum_i \frac{d}{dt}\left(\tfrac12 m_i \dot{\mathbf{R}}_i\cdot\dot{\mathbf{R}}_i\right) = \frac{dT}{dt}. \] On the right, the ordered double sum contains, for each unordered pair \(\{i,j\}\), the two terms \(\mathbf{F}_{ij}\cdot\dot{\mathbf{R}}_i\) and \(\mathbf{F}_{ji}\cdot\dot{\mathbf{R}}_j\). Newton’s third law \(\mathbf{F}_{ji}=-\mathbf{F}_{ij}\) lets them be combined: \[ \mathbf{F}_{ij}\cdot\dot{\mathbf{R}}_i + \mathbf{F}_{ji}\cdot\dot{\mathbf{R}}_j = \mathbf{F}_{ij}\cdot\left(\dot{\mathbf{R}}_i - \dot{\mathbf{R}}_j\right) = -\,\mathbf{F}_{ij}\cdot\dot{\mathbf{r}}_{ij}. \] This is the crucial step: the pair interaction does work only through the rate of change of the separation of that pair. Now use \(\mathbf{r}\cdot\dot{\mathbf{r}} = r\dot r\) (proved on the eccentricity-vector slide): \[ \mathbf{F}_{ij}\cdot\dot{\mathbf{r}}_{ij} = G m_i m_j \frac{\mathbf{r}_{ij}\cdot\dot{\mathbf{r}}_{ij}}{r_{ij}^{3}} = G m_i m_j \frac{\dot r_{ij}}{r_{ij}^{2}} = -\frac{d}{dt}\left(\frac{G m_i m_j}{r_{ij}}\right). \] Hence the pair contributes \[ -\mathbf{F}_{ij}\cdot\dot{\mathbf{r}}_{ij} = \frac{d}{dt}\left(\frac{Gm_im_j}{r_{ij}}\right) = -\frac{d\Phi_{ij}}{dt}, \qquad \Phi_{ij} = -\frac{Gm_im_j}{r_{ij}} . \]

(b) The integral. Summing over the \(\binom{n}{2}\) unordered pairs, \[ \frac{dT}{dt} = -\frac{d}{dt}\sum_{i<j}\Phi_{ij} = -\frac{dU}{dt} \qquad\Longrightarrow\qquad \frac{d}{dt}\left(T + U\right) = 0, \] so \(E = T + U = C_4\) is constant. Together with \(\mathbf{C}_1\) (3), \(\mathbf{C}_2\) (3) and \(\mathbf{C}_3\) (3) this completes the ten classical integrals.

(c) Why \(i<j\). Each pair carries one interaction energy. The ordered double sum \(\sum_i\sum_{j\ne i}\Phi_{ij}\) visits the pair \(\{i,j\}\) twice and therefore equals \(2U\). The forces are not double counted (body \(i\) genuinely feels \(\mathbf{F}_{ij}\) and body \(j\) feels \(\mathbf{F}_{ji}\)), but the energy is a property of the pair, not of a body. Getting this wrong doubles \(U\); in an N-body integrator this shows up immediately as an energy-conservation diagnostic that is off by a factor of two, and it corrupts any virial or escape test built on \(E\) (Problem 2). A useful check: for two bodies the correct \(U\) must reproduce \(\varepsilon = v^2/2 - \mu/r\) per unit reduced mass.

Solution 2 — Lagrange–Jacobi and Jacobi’s criterion

Work in the barycentric frame, where \(\ddot{\mathbf{R}}_{cm} = \mathbf{0}\) implies \(\ddot{\boldsymbol{\rho}}_i = \ddot{\mathbf{R}}_i\). Note also that the mutual separations are unchanged by the shift, \(\boldsymbol{\rho}_j - \boldsymbol{\rho}_i = \mathbf{r}_{ij}\), so \(U\) is the same function as before.

(a) The identity. Differentiating \(I = \sum_i m_i\,\boldsymbol{\rho}_i\cdot\boldsymbol{\rho}_i\) twice, \[ \dot I = 2\sum_i m_i\,\boldsymbol{\rho}_i\cdot\dot{\boldsymbol{\rho}}_i, \qquad \ddot I = 2\sum_i m_i\,\dot{\boldsymbol{\rho}}_i\cdot\dot{\boldsymbol{\rho}}_i + 2\sum_i \boldsymbol{\rho}_i\cdot\left(m_i\ddot{\boldsymbol{\rho}}_i\right) = 4T + 2\sum_i \boldsymbol{\rho}_i\cdot\mathbf{F}_i , \] where \(\mathbf{F}_i = \sum_{j\ne i}\mathbf{F}_{ij}\) is the total force on body \(i\). Reorganise the last sum over unordered pairs exactly as in Problem 1: \[ \sum_i \boldsymbol{\rho}_i\cdot\mathbf{F}_i = \sum_{i<j}\left(\boldsymbol{\rho}_i - \boldsymbol{\rho}_j\right)\cdot\mathbf{F}_{ij} = -\sum_{i<j}\mathbf{r}_{ij}\cdot\mathbf{F}_{ij} = -\sum_{i<j} Gm_im_j\frac{\mathbf{r}_{ij}\cdot\mathbf{r}_{ij}}{r_{ij}^{3}} = -\sum_{i<j}\frac{Gm_im_j}{r_{ij}} = U . \] (This quantity is called the virial; the calculation just performed is Euler’s theorem \(\sum_i \mathbf{R}_i\cdot\nabla_i U = -U\) for a potential homogeneous of degree \(-1\), and note it is origin-independent because \(\sum_i\mathbf{F}_i = \mathbf{0}\).) Therefore \[ \boxed{\;\ddot I = 4T + 2U\;} \qquad\text{and, using } T = E - U, \qquad \ddot I = 4E - 2U . \]

(b) Jacobi’s criterion. \(U<0\) always, so \(-2U = 2|U| > 0\): the identity says \(\ddot I = 4E + 2|U|\). Suppose the system stayed bounded, i.e. \(r_{ij}\le D\) for all pairs and all \(t\). Then \[ |U| = \sum_{i<j}\frac{Gm_im_j}{r_{ij}} \;\ge\; \frac{1}{D}\sum_{i<j}Gm_im_j \;\equiv\; c > 0 , \] a positive constant. Hence, if \(E\ge 0\), \[ \ddot I \;\ge\; 4E + 2c \;\ge\; 2c > 0 \quad\text{for all } t, \] and integrating twice, \[ I(t) \;\ge\; I(0) + \dot I(0)\,t + c\,t^{2} \;\longrightarrow\; \infty . \] But a bounded configuration has \(I = \sum m_i\rho_i^2 \le M D^2\), bounded. This is a contradiction, so no N-body system with \(E \ge 0\) can remain bounded: at least one body must escape (or the configuration must expand without limit). Note the converse is false — \(E<0\) does not guarantee that every body stays bound, since one body may escape while the rest form a tighter binary.

(c) Virial theorem. For bounded motion \(I\) and \(\dot I\) remain bounded. Average the identity over \([0,\tau]\): \[ \frac{1}{\tau}\int_0^\tau \ddot I\,dt = \frac{\dot I(\tau) - \dot I(0)}{\tau} \xrightarrow[\tau\to\infty]{} 0 , \] because the numerator is bounded. Hence \(\langle 4T + 2U\rangle = 0\), i.e. \[ \boxed{\;\langle 2T\rangle = -\langle U\rangle\;} \] and since \(E = \langle T\rangle + \langle U\rangle = \langle T\rangle - 2\langle T\rangle\), \[ E = -\langle T\rangle = \tfrac12\langle U\rangle . \] So a gravitationally bound system has half its potential energy in kinetic form — and, notably, adding energy to it makes it hotter is false: removing energy (\(E\) more negative) increases \(\langle T\rangle\). This negative heat capacity is the reason star clusters and stars contract and heat up as they radiate.

(d) Zwicky’s argument. With \(T = \tfrac12 M v_{rms}^2\) and, for a system of size \(R\), \(U = -\alpha\,GM^2/R\) with \(\alpha\) a dimensionless structure constant of order unity (\(\alpha = 3/5\) for a uniform sphere), \[ M v_{rms}^{2} = \frac{\alpha G M^{2}}{R} \qquad\Longrightarrow\qquad \boxed{\;M = \frac{v_{rms}^{2} R}{\alpha G}\;} \] Applying this to the Coma cluster in 1933, Zwicky found a mass some two orders of magnitude larger than the mass visible as luminous galaxies — the first evidence for dark matter. The entire inference rests on the identity proved in (a), which requires no solution of the N-body problem at all.

Solution 3 — The invariable plane

(a) König decomposition. Substitute \(\mathbf{R}_i = \mathbf{R}_{cm} + \boldsymbol{\rho}_i\) and \(\dot{\mathbf{R}}_i = \dot{\mathbf{R}}_{cm} + \dot{\boldsymbol{\rho}}_i\): \[ \mathbf{H} = \sum_i m_i\left(\mathbf{R}_{cm} + \boldsymbol{\rho}_i\right)\times\left(\dot{\mathbf{R}}_{cm} + \dot{\boldsymbol{\rho}}_i\right) = M\,\mathbf{R}_{cm}\times\dot{\mathbf{R}}_{cm} + \mathbf{R}_{cm}\times\underbrace{\sum_i m_i\dot{\boldsymbol{\rho}}_i}_{=\,\mathbf{0}} + \underbrace{\left(\sum_i m_i\boldsymbol{\rho}_i\right)}_{=\,\mathbf{0}}\times\dot{\mathbf{R}}_{cm} + \sum_i m_i\boldsymbol{\rho}_i\times\dot{\boldsymbol{\rho}}_i , \] the two middle terms vanishing by the definition of the barycentre. So \(\mathbf{H} = M\mathbf{R}_{cm}\times\dot{\mathbf{R}}_{cm} + \mathbf{H}_{cm}\). The orbital term is constant on its own, since \[ \frac{d}{dt}\left(M\mathbf{R}_{cm}\times\dot{\mathbf{R}}_{cm}\right) = M\dot{\mathbf{R}}_{cm}\times\dot{\mathbf{R}}_{cm} + M\mathbf{R}_{cm}\times\ddot{\mathbf{R}}_{cm} = \mathbf{0}, \] using \(\ddot{\mathbf{R}}_{cm}=\mathbf{0}\). Since \(\mathbf{H}\) is constant (slides) and the orbital term is constant, the difference \(\mathbf{H}_{cm}\) is constant too.

(b) Transformation properties. Under a shift of origin \(\mathbf{R}_i\mapsto\mathbf{R}_i-\mathbf{b}\) (with \(\mathbf{b}\) constant, so velocities unchanged), \[ \mathbf{H} \mapsto \sum_i m_i(\mathbf{R}_i - \mathbf{b})\times\dot{\mathbf{R}}_i = \mathbf{H} - \mathbf{b}\times\sum_i m_i\dot{\mathbf{R}}_i = \mathbf{H} - \mathbf{b}\times\mathbf{C}_1 . \] Under a boost \(\dot{\mathbf{R}}_i\mapsto\dot{\mathbf{R}}_i - \mathbf{u}\), \[ \mathbf{H}\mapsto \mathbf{H} - \left(\sum_i m_i\mathbf{R}_i\right)\times\mathbf{u} = \mathbf{H} - M\,\mathbf{R}_{cm}\times\mathbf{u}. \] Both changes are generally nonzero, and — crucially — they change the direction of \(\mathbf{H}\), not merely its magnitude. But \(\boldsymbol{\rho}_i = \mathbf{R}_i - \mathbf{R}_{cm}\) and \(\dot{\boldsymbol{\rho}}_i\) are differences, hence unchanged by both a shift (\(\mathbf{R}_{cm}\) shifts identically) and a boost (\(\dot{\mathbf{R}}_{cm}\) boosts identically). Therefore \(\mathbf{H}_{cm}\) is invariant under both.

(c) Conclusion. “The plane through the origin perpendicular to \(\mathbf{H}\)” depends on where you put the origin and how fast you are moving — two astronomers with different conventions would disagree about it. The plane through the barycentre perpendicular to \(\mathbf{H}_{cm}\) is built entirely from internal quantities and is therefore the same for every inertial observer; it is a property of the solar system itself. That is what earns it the name invariable. (In the barycentric frame, \(\mathbf{C}_1 = \mathbf{0}\) and \(\mathbf{R}_{cm}=\mathbf{0}\), so \(\mathbf{H} = \mathbf{H}_{cm}\) and the distinction disappears — which is exactly why we choose that frame.)

Solution 4 — Lagrange’s equilateral solution, unequal masses

(a) The acceleration is linear in \(\mathbf{R}_i\). Put the origin at the barycentre, so \(\sum_j m_j\mathbf{R}_j = \mathbf{0}\). The acceleration of body \(i\) is \[ \ddot{\mathbf{R}}_i = G\sum_{j\ne i} \frac{m_j}{r_{ij}^{3}}\left(\mathbf{R}_j - \mathbf{R}_i\right). \] Here is the only place the geometry enters: in an equilateral configuration all three separations are equal, \(r_{12}=r_{23}=r_{31}=d\), so the factor \(1/r_{ij}^3\) is common to every term and comes out of the sum: \[ \ddot{\mathbf{R}}_i = \frac{G}{d^{3}}\left[\sum_{j\ne i} m_j\mathbf{R}_j - \mathbf{R}_i\sum_{j\ne i}m_j\right] = \frac{G}{d^{3}}\Big[\underbrace{\Big(\textstyle\sum_j m_j\mathbf{R}_j\Big)}_{=\,\mathbf 0} - m_i\mathbf{R}_i - \left(M - m_i\right)\mathbf{R}_i\Big], \] and the \(m_i\mathbf{R}_i\) terms cancel: \[ \boxed{\;\ddot{\mathbf{R}}_i = -\frac{GM}{d^{3}}\,\mathbf{R}_i\;},\qquad M = m_1+m_2+m_3 . \] Every body is accelerated exactly towards the barycentre, with magnitude proportional to its own distance \(\rho_i\) from it — even though the three distances \(\rho_i\) are different and the three masses are different.

(b) The rotation rate. Uniform circular motion of body \(i\) about the origin at rate \(\omega\) requires \(\ddot{\mathbf{R}}_i = -\omega^{2}\mathbf{R}_i\). Comparing, \[ \boxed{\;\omega^{2} = \frac{GM}{d^{3}}\;} \] The same \(\omega\) for all three bodies — which is what makes rigid rotation consistent — and it depends only on the total mass and the side length. Note the configuration is not symmetric (the heaviest body sits closest to the barycentre, on the smallest circle), yet all three complete a revolution in the same time, which is precisely the condition for the triangle to stay equilateral.

Check (Assignment 01): for \(m_1=m_2=m_3=m\) the barycentre is the centroid, \(\rho = d/\sqrt3\), and the formula gives \(\omega_0^2 = 3Gm/d^3\), i.e. \(\omega_0 = \sqrt{3Gm/d_0^{3}}\) — as found there.

(c) Non-circular equilateral solutions. Since \(\ddot{\mathbf{R}}_i = -(GM/d^3)\mathbf{R}_i\) holds instantaneously for any equilateral configuration, it is not restricted to circles. Write \(\mathbf{R}_i = \sigma_i\,\mathbf{q}(t)\) for a common vector \(\mathbf{q}\) scaled by fixed factors: each body’s equation reduces to the same Kepler-type equation. Consequently the triangle may pulsate — expanding and contracting while remaining equilateral — with each body tracing a similar Keplerian conic about the barycentre, all with the same eccentricity and period. These are Lagrange’s homographic solutions; the circular case is the special \(e = 0\) member. Combined with Euler’s collinear solutions, they are the only known exact solutions of the general three-body problem, and the \(m_3 \to 0\) limit of the equilateral family gives the L4 and L5 points — where the Trojan asteroids sit.

Part B — The Two-body Problem

Solution 5 — Binet’s equation

(a) Change of independent variable. Angular-momentum conservation gives \(\dot f = h u^{2}\) with \(u = 1/r\). Then, writing \(u' = du/df\), \[ \dot r = \frac{d}{dt}\left(\frac1u\right) = -\frac{1}{u^{2}}\,u'\,\dot f = -\frac{1}{u^{2}}\,u'\,\left(hu^{2}\right) = -h\,u' , \] \[ \ddot r = \frac{d}{dt}\left(-h u'\right) = -h\,u''\,\dot f = -h^{2}u^{2}\,u'' . \] Note how \(h\) absorbs every trace of time: the \(r^{-2}\) non-linearity of the problem has been traded for constant coefficients.

(b) The equation. The radial component of \(\ddot{\mathbf{r}} = -\mu\mathbf{r}/r^3\) in plane polars is \(\ddot r - r\dot f^{2} = -\mu/r^{2}\). Since \(r\dot f^{2} = u^{-1}(hu^{2})^{2} = h^{2}u^{3}\), \[ -h^{2}u^{2}u'' - h^{2}u^{3} = -\mu u^{2} \quad\xrightarrow{\;\div\,(-h^{2}u^{2})\;}\quad \boxed{\;u'' + u = \frac{\mu}{h^{2}}\;} \] A linear, constant-coefficient, inhomogeneous oscillator equation — the whole Kepler problem, in one line. (That it is linear at all is a peculiarity of the inverse square: any other exponent leaves a power of \(u\) on the right.)

(c) Solution and identification of the constants. The general solution is the particular solution \(\mu/h^{2}\) plus the homogeneous solution \(A\cos(f - f_0)\): \[ u = \frac{\mu}{h^{2}} + A\cos\left(f - f_0\right) \qquad\Longrightarrow\qquad r = \frac{h^{2}/\mu}{1 + \dfrac{Ah^{2}}{\mu}\cos\left(f - f_0\right)} . \] This is a conic with semi-latus rectum \(p = h^{2}/\mu\) and eccentricity \(e = Ah^{2}/\mu\); \(r\) is smallest when \(f = f_0\), so \(f_0\) is the direction of periapsis. The two constants of integration are thus \((A, f_0)\) = (size of the radial oscillation, orientation of the apse line).

To identify them with \(\mathbf{e}\), dot the slides’ result \(\dot{\mathbf{r}}\times\mathbf{h} = \mu\left(\hat{\mathbf{r}} + \mathbf{e}\right)\) with \(\mathbf{r}\): \[ \mathbf{r}\cdot\left(\dot{\mathbf{r}}\times\mathbf{h}\right) = \mathbf{h}\cdot\left(\mathbf{r}\times\dot{\mathbf{r}}\right) = h^{2} = \mu\left(r + \mathbf{e}\cdot\mathbf{r}\right) = \mu\left(r + e\,r\cos\theta\right), \] where \(\theta\) is the angle from \(\mathbf{e}\) to \(\mathbf{r}\). Hence \(r = (h^{2}/\mu)/(1 + e\cos\theta)\). Comparing with the Binet solution term by term: \[ e = \frac{Ah^{2}}{\mu} \;=\;|\mathbf{e}|, \qquad f_0 = \text{direction of } \hat{\mathbf{e}} , \] so \(\theta = f - f_0\) and the constant vector \(\mathbf{e}\) is nothing but the pair \((A, f_0)\) packaged as a vector. The LRL vector is therefore not extra information; it is the two integration constants of Binet’s equation, written covariantly.

(d) The sixth constant. Binet’s equation gives the orbit as a curve, \(r(f)\), and knows nothing about time — \(t\) was eliminated in step (a). Time re-enters only when one integrates \(dt = r^{2}df/h\), whose own constant of integration is \(t_p\). Counting: \(\mathbf{h}\) (3: orientation of the plane, plus \(h\) which fixes \(p\)), \(e\) and \(f_0\) within that plane, i.e. \(5\) constants for the geometry, plus \(t_p\) for the timing — exactly the “\(5 + 1 = 6\)” of the Counting the integrals slide. The Casimir relations are automatically respected here, since we never introduced a redundant constant in the first place.

Solution 6 — Apsidal precession

(a) Perturbed Binet equation. The perturbation is still central (\(\propto\hat{\mathbf{r}}\)), so \(\mathbf{r}\times\ddot{\mathbf{r}} = \mathbf{0}\) and \(\mathbf{h} = \mathbf{r}\times\dot{\mathbf{r}}\) is conserved exactly as before: the motion stays planar and \(h = r^2\dot f\) still holds, so the whole substitution of Problem 5 goes through unchanged. The radial equation is now \[ \ddot r - r\dot f^{2} = -\left(\mu u^{2} + \delta u^{3}\right) \quad\Longrightarrow\quad -h^{2}u^{2}u'' - h^{2}u^{3} = -\mu u^{2} - \delta u^{3} . \] Dividing by \(-h^{2}u^{2}\), \[ u'' + u = \frac{\mu}{h^{2}} + \frac{\delta}{h^{2}}u \qquad\Longrightarrow\qquad \boxed{\;u'' + k^{2}u = \frac{\mu}{h^{2}}, \qquad k^{2} = 1 - \frac{\delta}{h^{2}}\;} \] The \(1/r^{3}\) force is exactly the one perturbation that keeps the equation linear — it shifts the frequency of the radial oscillation instead of distorting it.

(b) Precession. With \(f\) measured from a periapsis, \[ u = \frac{\mu}{h^{2}k^{2}}\left(1 + e\cos kf\right) \qquad\Longrightarrow\qquad r = \frac{p'}{1 + e\cos kf}, \qquad p' = \frac{h^{2}k^{2}}{\mu} . \] The radius is periodic in \(f\) with period \(2\pi/k\), not \(2\pi\): successive periapsis passages occur at \(kf = 0, 2\pi, 4\pi,\ldots\), i.e. \(\Delta f = 2\pi/k\). The apse line therefore rotates by \[ \Delta\varpi = \frac{2\pi}{k} - 2\pi = 2\pi\left(k^{-1} - 1\right) \;\approx\; 2\pi\cdot\frac{\delta}{2h^{2}} = \frac{\pi\delta}{h^{2}}, \] expanding \(k^{-1} = (1-\delta/h^2)^{-1/2} \approx 1 + \delta/(2h^{2})\). For \(\delta > 0\) (extra attraction) \(k<1\), so \(\Delta\varpi > 0\): the orbit must turn through more than \(2\pi\) to return to periapsis and the apse line advances in the direction of motion. For \(\delta<0\) it regresses. Only \(\delta = 0\) gives \(k=1\) and a closed ellipse — which is the quantitative form of the slides’ statement that \(\dot{\mathbf{e}} = \mathbf{0}\) is special to the inverse square.

(c) The relativistic term. Write \(u = u_0 + u_1\), with \(u_0 = (\mu/h^{2})(1 + e\cos f)\) the Kepler solution and \(u_1\) small. To first order, \[ u_1'' + u_1 = \frac{3\mu}{c^{2}}u_0^{2} = \frac{3\mu^{3}}{c^{2}h^{4}}\left(1 + 2e\cos f + e^{2}\cos^{2}f\right). \] Of the three terms on the right, the constant and the \(\cos^2 f = \tfrac12(1+\cos 2f)\) terms force the oscillator off resonance and produce bounded, periodic corrections — tiny changes in the size and shape of the orbit, of no long-term consequence. The term \(2e\cos f\) is resonant with the left-hand side, and its particular solution grows without bound: \[ u_1'' + u_1 = C\cos f \;\Longrightarrow\; u_1 = \frac{C}{2}\,f\sin f, \qquad C = \frac{6e\mu^{3}}{c^{2}h^{4}} . \] Keeping only this secular piece, \[ u \approx \frac{\mu}{h^{2}}\left[1 + e\cos f + \alpha\,e\,f\sin f\right], \qquad \alpha = \frac{3\mu^{2}}{c^{2}h^{2}} . \] Since \(\alpha \ll 1\), use \(\cos\left[(1-\alpha)f\right] = \cos f\cos\alpha f + \sin f \sin\alpha f \approx \cos f + \alpha f\sin f\) to resum: \[ u \approx \frac{\mu}{h^{2}}\Big[1 + e\cos\big((1-\alpha)f\big)\Big] , \] which is case (b) with \(k = 1-\alpha\). Hence, using \(h^{2} = \mu p\), \[ \Delta\varpi = 2\pi\left(k^{-1}-1\right) \approx 2\pi\alpha = \frac{6\pi\mu^{2}}{c^{2}h^{2}} = \boxed{\;\frac{6\pi\mu}{c^{2}p}\;}\;\text{per revolution.} \]

Mercury. \(p = a(1-e^{2}) = 5.79\times10^{10}\times\left(1-0.2056^{2}\right) = 5.545\times10^{10}\,\mathrm{m}\), so \[ \Delta\varpi = \frac{6\pi\left(1.327\times10^{20}\right)}{\left(2.998\times10^{8}\right)^{2}\left(5.545\times10^{10}\right)} = 5.02\times10^{-7}\ \mathrm{rad/rev}. \] Revolutions per Julian century: \(36525/87.97 = 415.2\). Therefore \[ \Delta\varpi_{\text{century}} = 5.02\times10^{-7}\times415.2 = 2.08\times10^{-4}\,\mathrm{rad} = 2.08\times10^{-4}\times 206265'' \approx 43''\ \text{per century}, \] the famous residual left over after all Newtonian planetary perturbations were accounted for, and the first confirmed prediction of general relativity.

Solution 7 — The hyperbolic Kepler equation

(a) The parametrisation. The attracting focus is at the origin of the perifocal frame and periapsis lies on the \(+x_P\) axis at \(x_P = \bar a(e-1) = p/(1+e)\), which is what the given formulae return at \(H = 0\). The centre of the hyperbola lies a further \(\bar a\) along, at \(x_P = \bar a e\). Relative to the centre, \[ X = x_P - \bar a e = -\bar a\cosh H, \qquad Y = y_P = b\sinh H, \qquad b = \bar a\sqrt{e^{2}-1}, \] so that \(X^{2}/\bar a^{2} - Y^{2}/b^{2} = \cosh^{2}H - \sinh^{2}H = 1\): the standard hyperbola, and since \(X \le -\bar a\) it is the branch on the far side of the centre from the focus — the branch that wraps around the attracting body, as gravity requires. (\(H<0\) is the inbound leg, \(H>0\) the outbound.)

The radius follows exactly as for the ellipse: \[ \begin{aligned} r^{2} &= \bar a^{2}\left(e - \cosh H\right)^{2} + \bar a^{2}\left(e^{2}-1\right)\sinh^{2}H\\ &= \bar a^{2}\left[e^{2} - 2e\cosh H + \cosh^{2}H + \left(e^{2}-1\right)\left(\cosh^{2}H - 1\right)\right]\\ &= \bar a^{2}\left[e^{2}\cosh^{2}H - 2e\cosh H + 1\right] = \bar a^{2}\left(e\cosh H - 1\right)^{2}, \end{aligned} \] giving \(\;\boxed{r = \bar a\left(e\cosh H - 1\right)}\;\) — again linear in the new variable, which is the entire point. Consistency check with the orbit equation: \(\cos f = x_P/r\) gives \[ 1 + e\cos f = \frac{e\cosh H - 1 + e^{2} - e\cosh H}{e\cosh H - 1} = \frac{e^{2}-1}{e\cosh H - 1}, \] so \(p/(1+e\cos f) = \bar a(e^2-1)\cdot\frac{e\cosh H-1}{e^2-1} = r\). ✓

(b) The time equation. With \(dx_P = -\bar a\sinh H\,dH\) and \(dy_P = \bar a\sqrt{e^{2}-1}\cosh H\,dH\), \[ \begin{aligned} r^{2}df = x_P\,dy_P - y_P\,dx_P &= \bar a^{2}\sqrt{e^{2}-1}\left[\left(e-\cosh H\right)\cosh H + \sinh^{2}H\right]dH\\ &= \bar a^{2}\sqrt{e^{2}-1}\left(e\cosh H - 1\right)dH , \end{aligned} \] using \(\sinh^2 H - \cosh^2 H = -1\). With \(h = \sqrt{\mu p} = \sqrt{\mu\bar a}\sqrt{e^{2}-1}\), the \(\sqrt{e^{2}-1}\) cancels exactly as \(\sqrt{1-e^{2}}\) did in the elliptic case: \[ dt = \frac{r^{2}df}{h} = \sqrt{\frac{\bar a^{3}}{\mu}}\left(e\cosh H - 1\right)dH . \] Integrating from periapsis (\(H = 0\) at \(t = t_p\)), \[ \boxed{\;M_h \equiv \sqrt{\frac{\mu}{\bar a^{3}}}\left(t - t_p\right) = e\sinh H - H\;} \] Unlike \(M = E - e\sin E\), this is monotonic and unbounded — the body passes periapsis once and never returns — but it is equally transcendental, and is solved by the same Newton–Raphson iteration \(H_{k+1} = H_k - \dfrac{e\sinh H_k - H_k - M_h}{e\cosh H_k - 1}\).

(c) Half-angle relation. Use \(\tan\frac f2 = \dfrac{\sin f}{1 + \cos f} = \dfrac{r\sin f}{r + r\cos f}\): \[ r\sin f = y_P = \bar a\sqrt{e^{2}-1}\,\sinh H, \qquad r + r\cos f = \bar a\left(e\cosh H - 1\right) + \bar a\left(e - \cosh H\right) = \bar a\left(e-1\right)\left(\cosh H + 1\right). \] Dividing, and using \(\tanh\frac H2 = \dfrac{\sinh H}{1+\cosh H}\), \[ \tan\frac f2 = \frac{\sqrt{(e-1)(e+1)}}{e-1}\cdot\frac{\sinh H}{1 + \cosh H} = \sqrt{\frac{e+1}{e-1}}\,\tanh\frac{H}{2} \qquad\Longleftrightarrow\qquad \boxed{\;\tanh\frac{H}{2} = \sqrt{\frac{e-1}{e+1}}\tan\frac f2\;} \] the elliptic relation with \(\tan\to\tanh\) and \(e\pm1\) interchanged — as it must be, since \(E \to iH\) maps one to the other.

(d) Asymptotic quantities. \(r\to\infty\) requires the denominator of \(r = p/(1+e\cos f)\) to vanish: \[ \cos f_\infty = -\frac1e . \] (Consistently, \(\tanh(H/2) \to \pm1\) as \(H\to\pm\infty\) gives the same \(f_\infty\).) From \(\varepsilon = -\mu/(2a) = +\mu/(2\bar a)\) and vis-viva as \(r\to\infty\), \[ \frac{v_\infty^{2}}{2} = \varepsilon \qquad\Longrightarrow\qquad v_\infty = \sqrt{\frac{\mu}{\bar a}} . \] The incoming and outgoing asymptotes make angles \(\pm f_\infty\) with the apse line, so the velocity direction is turned through \(\delta = 2f_\infty - 180^{\circ}\), whence \[ \sin\frac{\delta}{2} = \sin\left(f_\infty - 90^\circ\right) = -\cos f_\infty = \frac{1}{e} . \] This is the master equation of the gravity-assist manoeuvre: a low-\(e\) (close, slow) flyby turns the velocity sharply; \(e\to\infty\) leaves it undeflected.

Solution 8 — Barker’s equation

(a) The integral. With \(e = 1\) and \(1 + \cos f = 2\cos^{2}\frac f2\), \[ t - t_p = \frac{p^{2}}{h}\int_0^f \frac{df'}{4\cos^{4}\frac{f'}{2}} . \] Substitute \(D = \tan\frac{f'}{2}\), so \(dD = \tfrac12\sec^{2}\frac{f'}2\,df'\), i.e. \(df' = 2\cos^{2}\frac{f'}{2}\,dD\): \[ \frac{df'}{4\cos^{4}\frac{f'}{2}} = \frac{2\cos^{2}\frac{f'}2\,dD}{4\cos^{4}\frac{f'}2} = \frac{1}{2}\sec^{2}\frac{f'}{2}\,dD = \frac{1}{2}\left(1 + D^{2}\right)dD . \] Hence \(\int_0^f = \frac12\left(D + \frac{D^{3}}{3}\right)\) with \(D = \tan\frac f2\). The prefactor simplifies with \(h = \sqrt{\mu p}\): \[ \frac{p^{2}}{h} = \frac{p^{2}}{\sqrt{\mu p}} = \sqrt{\frac{p^{3}}{\mu}} , \] giving Barker’s equation \[ \boxed{\;t - t_p = \frac{1}{2}\sqrt{\frac{p^{3}}{\mu}}\left(D + \frac{D^{3}}{3}\right), \qquad D = \tan\frac f2 \;} \]

(b) Inversion. Write \(B = 6\left(t-t_p\right)\sqrt{\mu/p^{3}}\), so the equation is the depressed cubic \[ D^{3} + 3D = B . \] The function \(D\mapsto D^{3}+3D\) has derivative \(3D^{2}+3 > 0\) everywhere, so it is strictly increasing from \(-\infty\) to \(+\infty\): exactly one real root for every \(B\), and the map \(t \leftrightarrow f\) is a bijection. (Equivalently, Cardano’s discriminant \(\left(B/2\right)^{2} + 1\) is strictly positive.) Cardano’s formula gives \[ D = \sqrt[3]{\frac{B}{2} + \sqrt{\frac{B^{2}}{4} + 1}} \;+\; \sqrt[3]{\frac{B}{2} - \sqrt{\frac{B^{2}}{4} + 1}} , \] or, more robustly in floating point, using \(\sinh 3x = 3\sinh x + 4\sinh^{3}x\) with \(D = 2\sinh x\) (which turns the cubic into \(2\sinh 3x = B\)), \[ \boxed{\;D = 2\sinh\left[\frac{1}{3}\operatorname{arcsinh}\frac{B}{2}\right]\;} \qquad\Longrightarrow\qquad f = 2\arctan D . \] No iteration, no Newton–Raphson, no initial guess.

(c) Why \(e = 1\) is different. For \(e\ne1\) the change of variable that linearises the integrand (\(E\) or \(H\)) buys a linear integrand at the price of a transcendental relation between the new variable and time: \(M = E - e\sin E\) mixes an algebraic term with a trigonometric one, and no finite combination of elementary functions inverts it (a theorem of Liouville-type; the slides’ aside). At \(e = 1\) the auxiliary circle degenerates — the ellipse’s centre runs off to infinity, \(a\to\infty\) and \(n\to0\), so \(E\) and \(M\) both become useless — but the half-angle variable \(D = \tan(f/2)\) takes their place, and the integral produces a polynomial in \(D\) rather than a mixture of a polynomial and a sinusoid. A polynomial of degree three is solvable in radicals; that is the whole story. The parabola is thus the only conic on which position-versus-time is a closed-form function — which is why near-parabolic comet orbits are traditionally handled with Barker’s equation plus a correction, rather than with Kepler’s equation.

Solution 9 — Rectilinear fall (\(h = 0\))

(a) Geometry of the degenerate orbit. Relative motion obeys \(\ddot{\mathbf{r}} = -\mu\mathbf{r}/r^{3}\) with \(\mu = G(m_1+m_2)\). At \(t=0\), \(\dot{\mathbf{r}}(0)=\mathbf{0}\), so \[ \mathbf{h} = \mathbf{r}(0)\times\dot{\mathbf{r}}(0) = \mathbf{0}, \] and since \(\mathbf{h}\) is constant it stays zero. As \(\ddot{\mathbf{r}}\) is always parallel to \(\mathbf{r}\), no transverse velocity can ever be generated: the motion is confined to the fixed line joining the bodies. The energy is \[ \varepsilon = \frac{v^{2}}{2} - \frac{\mu}{r} = 0 - \frac{\mu}{r_0} = -\frac{\mu}{r_0} < 0 \quad\Longrightarrow\quad a = -\frac{\mu}{2\varepsilon} = \frac{r_0}{2}, \] finite, while the Casimir relation gives \[ e = \sqrt{1 + \frac{2\varepsilon h^{2}}{\mu^{2}}} = 1 \quad (\text{since } h = 0), \qquad p = \frac{h^{2}}{\mu} = 0 . \] So \(e = 1\) but \(a\) finite: not a parabola (which has \(a\to\infty\)) but a degenerate ellipse — an ellipse of zero width, i.e. the straight segment of length \(2a = r_0\) traversed from the far end (\(r_a = a(1+e) = r_0\) ✓) into the focus. This is the rectilinear ellipse.

(b) Direct quadrature. From vis-viva, \[ v^{2} = 2\left(\varepsilon + \frac{\mu}{r}\right) = 2\mu\left(\frac1r - \frac1{r_0}\right), \qquad \dot r = -\sqrt{2\mu\left(\frac1r - \frac1{r_0}\right)}\;\;(\text{infalling}). \] Hence \[ t_{\text{fall}} = \int_0^{r_0}\frac{dr}{\sqrt{2\mu\left(\frac1r - \frac1{r_0}\right)}} = \frac{1}{\sqrt{2\mu}}\int_0^{r_0}\sqrt{\frac{r\,r_0}{r_0 - r}}\;dr . \] Put \(r = r_0\sin^{2}\theta\), so \(dr = 2r_0\sin\theta\cos\theta\,d\theta\), \(r_0 - r = r_0\cos^{2}\theta\), and \(\dfrac{r r_0}{r_0-r} = r_0\tan^{2}\theta\); the limits \(r: 0\to r_0\) become \(\theta: 0\to\pi/2\): \[ t_{\text{fall}} = \frac{1}{\sqrt{2\mu}}\int_0^{\pi/2}\sqrt{r_0}\tan\theta\cdot2r_0\sin\theta\cos\theta\,d\theta = \frac{2r_0^{3/2}}{\sqrt{2\mu}}\int_0^{\pi/2}\sin^{2}\theta\,d\theta = \frac{2r_0^{3/2}}{\sqrt{2\mu}}\cdot\frac{\pi}{4} , \] that is \[ \boxed{\;t_{\text{fall}} = \frac{\pi}{2}\sqrt{\frac{r_0^{3}}{2\mu}}\;} \] (The integrand diverges as \(r\to0\), but only as \(r^{-1/2}\), which is integrable — the collision happens at finite time even though the speed blows up.)

(c) One-line Kepler argument. The motion is half a revolution of an ellipse with \(a = r_0/2\): the body starts at apoapsis \(r_a = r_0\) and reaches the focus, which is the \(e\to1\) limit of periapsis. By Kepler’s third law, \[ T = 2\pi\sqrt{\frac{a^{3}}{\mu}} = 2\pi\sqrt{\frac{r_0^{3}}{8\mu}}, \qquad t_{\text{fall}} = \frac{T}{2} = \pi\sqrt{\frac{r_0^{3}}{8\mu}} = \frac{\pi}{2}\sqrt{\frac{r_0^{3}}{2\mu}} \;\;\checkmark \] identical to (b). Equivalently, in Kepler’s equation with \(r = a(1-\cos E)\) and \(M = E - \sin E\): start at apoapsis \(E = \pi\)… measuring instead from periapsis, collision corresponds to \(E = \pi\), \(M = \pi\), so \(t = \pi/n = T/2\). The agreement of the two routes is the useful check — the degenerate conic really is a member of the Kepler family, not a separate case.

(d) Earth into the Sun. With \(r_0 = 1\,\mathrm{AU} = 1.496\times10^{11}\,\mathrm{m}\) and \(\mu \simeq \mu_\odot = 1.327\times10^{20}\,\mathrm{m^{3}s^{-2}}\), \[ t_{\text{fall}} = \frac{\pi}{2}\sqrt{\frac{\left(1.496\times10^{11}\right)^{3}}{2\times1.327\times10^{20}}} = \frac{\pi}{2}\left(3.55\times10^{6}\right)\,\mathrm{s} = 5.58\times10^{6}\,\mathrm{s} \approx 64.6\ \text{days}. \] Neatly, \(t_{\text{fall}} = T_\oplus/(4\sqrt2) = 365.25/5.657\) days — stop the Earth in its orbit and it takes about two months to fall into the Sun.

Solution 10 — Three averages of \(r\)

Throughout, \(r = p/(1+e\cos f) = a(1-e\cos E)\), \(n = 2\pi/T = \sqrt{\mu/a^{3}}\).

(a) Average over true anomaly. \[ \left\langle r\right\rangle_f = \frac{1}{2\pi}\int_0^{2\pi}\frac{p\,df}{1+e\cos f} . \] Evaluate with the Weierstrass substitution \(s = \tan\frac f2\), \(\cos f = \frac{1-s^{2}}{1+s^{2}}\), \(df = \frac{2\,ds}{1+s^{2}}\). Over the half-revolution \(f: 0\to\pi\) (\(s: 0\to\infty\)), \[ \int_0^{\pi}\frac{df}{1+e\cos f} = \int_0^{\infty}\frac{2\,ds}{\left(1+e\right)+\left(1-e\right)s^{2}} = \frac{2}{\sqrt{\left(1+e\right)\left(1-e\right)}}\left[\arctan\left(s\sqrt{\tfrac{1-e}{1+e}}\right)\right]_0^{\infty} = \frac{\pi}{\sqrt{1-e^{2}}} , \] and by symmetry the full revolution gives \(2\pi/\sqrt{1-e^{2}}\). Hence \[ \left\langle r\right\rangle_f = \frac{p}{\sqrt{1-e^{2}}} = \frac{a\left(1-e^{2}\right)}{\sqrt{1-e^{2}}} = a\sqrt{1-e^{2}} = b , \] the semi-minor axis.

(b) Average over eccentric anomaly. \[ \left\langle r\right\rangle_E = \frac{1}{2\pi}\int_0^{2\pi} a\left(1-e\cos E\right)dE = a , \] since \(\int_0^{2\pi}\cos E\,dE = 0\). So \(a\) is a mean distance — but only with respect to \(E\), a variable with no direct physical meaning.

(c) Time average. From Kepler’s equation, \(dM = n\,dt\) and \(M = E - e\sin E\) give \(dt = \left(1-e\cos E\right)dE/n\). Therefore \[ \left\langle r\right\rangle_t = \frac{1}{T}\int_0^{T} r\,dt = \frac{1}{T}\int_0^{2\pi} a\left(1-e\cos E\right)\cdot\frac{\left(1-e\cos E\right)}{n}\,dE = \frac{a}{2\pi}\int_0^{2\pi}\left(1-e\cos E\right)^{2}dE , \] using \(nT = 2\pi\). Expanding, \(\int_0^{2\pi}\left(1 - 2e\cos E + e^{2}\cos^{2}E\right)dE = 2\pi + \pi e^{2}\), so \[ \boxed{\;\left\langle r\right\rangle_t = a\left(1 + \frac{e^{2}}{2}\right)\;} \]

(d) Ordering and physical reading. Since \(\sqrt{1-e^{2}} \le 1 \le 1 + e^{2}/2\), \[ \left\langle r\right\rangle_f = a\sqrt{1-e^{2}} \;\;\le\;\; \left\langle r\right\rangle_E = a \;\;\le\;\; \left\langle r\right\rangle_t = a\left(1+\tfrac{e^{2}}{2}\right), \] with equality only for the circle. The reason is what each average weights:

• \(df\) weighting counts each degree of sweep equally. But by Kepler’s second law the body sweeps angle fastest near periapsis, so equal weight per degree over-represents the part of the orbit where \(r\) is small — the smallest average.
• \(dt\) weighting counts each second equally, and the body spends most of its time crawling near apoapsis — the largest average. (For a comet with \(e\to1\) the time average tends to \(3a/2\), dominated by the long slow outer arc.)
• \(dE\) sits between the two, which is exactly why \(a\) is so often quoted, and so often misdescribed as “the mean distance”.

(e) \(\langle 1/r\rangle_t\) and the virial theorem. The factor \((1-e\cos E)\) from \(dt\) now cancels the one in \(1/r\): \[ \left\langle \frac1r\right\rangle_t = \frac{1}{T}\int_0^{2\pi}\frac{1}{a\left(1-e\cos E\right)}\cdot\frac{\left(1-e\cos E\right)}{n}\,dE = \frac{1}{2\pi a}\int_0^{2\pi}dE = \frac{1}{a} . \] Then vis-viva \(v^{2} = \mu\left(\frac2r - \frac1a\right)\) time-averages to \[ \left\langle v^{2}\right\rangle_t = \mu\left(2\cdot\frac1a - \frac1a\right) = \frac{\mu}{a} . \] Now check against Problem 2(c). Per unit reduced mass, the two-body kinetic and potential energies are \(T = v^{2}/2\) and \(U = -\mu/r\), so \[ \left\langle 2T\right\rangle = \left\langle v^{2}\right\rangle = \frac{\mu}{a}, \qquad -\left\langle U\right\rangle = \mu\left\langle \frac1r\right\rangle = \frac{\mu}{a} . \qquad\checkmark \] And \(E = -\langle T\rangle = -\mu/(2a) = \varepsilon\), which is precisely the vis-viva constant. The virial theorem, proved for \(n\) bodies without solving anything, is therefore exactly consistent with the closed-form two-body solution — a satisfying consistency check on Part A.