Past exam of the mathematics course of the University of Cambridge / 2018 / iii / Paper 316 (Planetary System Dynamics)

Put ${\mu }_{\star }=G{M}_{\star }$. Conservation of specific orbital energy and specific angular momentum givesThe last comet expression is the vis-viva equation at the planet's orbit. The planet has only tangential velocity, while the comet's tangential component is $v_t=h/a_p$. Consequently its prograde relative velocity is exactly, before planetary deflection,In the near-parabolic encounter limit $Q\gg a_p$, this reduces toThe displayed approximation also requires $Q\gg a_p$. The assumptions $Q\gg q$ and $q\ll a_p$ alone do not guarantee it: an orbit with $Q=a_p$ encounters the planet near apoapsis and has $v_{\mathrm{rel}}\simeq v_p$, not $\sqrt{3}{v}_p$.

If all coplanar near-parabolic planet-crossing comets are admitted, so that $0\le q\le a_p$, the tangential component ranges from $-\sqrt{2}{v}_p$ to $+\sqrt{2}{v}_p$ when retrograde orbits are included. ThusFor prograde orbits alone the upper limit is $\sqrt{3}{v}_p$; for the stated deep-crossing limit $q\ll a_p$, both prograde and retrograde encounters approach $\sqrt{3}{v}_p$. These are the speeds outside the planet's gravitational well. The surface impact speed is $v_{\mathrm{imp}}=\sqrt{v_{\mathrm{rel}}^2+v_{\mathrm{esc},\mathrm{p}}^2}$, where $v_{\mathrm{esc},\mathrm{p}}$ is the planetary escape velocity; under negligible gravitational focusing the two speeds agree.

For an unperturbed Kepler orbit, successive revolutions advance the planet-relative orientation by $-n_p{T}_c$ modulo $2\pi$, where $T_c$ is the comet's orbital period. An irrational period ratio makes these orientations dense; sufficiently distant mean-motion resonance and adequate observation time justify the corresponding phase mixing approximation.

At each radius between its pericentre distance and apocentre distance, the comet eventually visits every azimuth in the rotating reference frame. The spatial projection is thereforeIt is an annulus, with a nonuniform radial residence probability, rather than a uniformly filled area. Strictly, the full position-velocity phase space is not this annulus: at each radius the energy and angular momentum constrain the velocity, with inward and outward branches. The annulus is its position-space projection. This description assumes the orbital elements have not yet been substantially changed by planetary scattering.

The radial velocity of a Kepler orbit follows from the vis-viva equation after subtracting the tangential component:Each radial interval is crossed twice per orbital period, so its fraction of time isThis is the phase-mixed radial probability of a Kepler orbit, normalized to one between $q$ and $Q$. Dividing by the annular area $2\pi rdr$ gives the phase-mixed comet surface density per comet:For $Q\gg a_p>q$,The printed density has the powers of $Q$ and $a_p$ interchanged. Its expression is too large by $Q/a_p$ in this limit. The normalized residence-time derivation fixes the corrected expression, which is also required for the later $Q^{1/2}$ ejection-time scaling. For a population of $N_c$ independent identical comets, multiply ${\Sigma }_c$ by $N_c$ to obtain a number surface density of a disk.

With negligible comet radius and mass, the planet's radius and escape velocity areGravitational focusing changes the limiting collision impact parameter from $R_p$ to $R_p\sqrt{1+v_{\mathrm{esc},\mathrm{p}}^2/v_{\mathrm{rel}}^2}$. Thus it is negligible when $v_{\mathrm{esc},\mathrm{p}}^2\ll v_{\mathrm{rel}}^2$.

Put $s=q/a_p$ and $\gamma =\sqrt{3-2\sqrt{2s}}$ in the $Q\gg a_p$ prograde approximation. Substituting $v_{\mathrm{rel}}=\gamma \sqrt{G{M}_{\star }/a_p}$ givesEquivalently $2(M_p/M_{\star })(a_p/R_p)\ll {\gamma }^2$. For $q\ll a_p$, use ${\gamma }^2\simeq 3$. If negligible gravitational focusing is required for every coplanar near-parabolic orbit, including nearly tangential prograde encounters, use the smaller value ${\gamma }^2=(\sqrt{2}-1{)}^2$.

A coplanar target sweeps a strip of width $2R_p$ through the comet's local position probability. Hence the relevant geometric collision cross-section is a length, and the rare collision rate per comet is ${\Gamma }_{\mathrm{col}}=2R_p{\Sigma }_c(a_p)v_{\mathrm{rel}}$.

Using the corrected phase-mixed comet surface density and the notation $s=q/a_p$, $\gamma =\sqrt{3-2\sqrt{2s}}$, the coplanar comet collision time isFor $q\ll a_p$, this becomes ${\pi }^2{Q}^{3/2}{a}_p/(2\sqrt{3}{R}_p\sqrt{G{M}_{\star }})$, so $t_{\mathrm{col}}\propto Q^{3/2}{M}_p^{-1/3}{a}_p$ at fixed stellar mass and planetary mass density. This is the mean waiting time for one orbiting comet; the mean interval between impacts from $N_c$ such independent comets is $t_{\mathrm{col}}/N_c$. It assumes rare encounters and orbital phase mixing, not a Poisson process derived from a single deterministic trajectory.

Suppose the inclined population samples its arguments of periapsis and nodal orientations, or that sufficiently slow precession supplies this sampling. Near the planet the vertical amplitude is $H\simeq a_{p}I$. For uniformly distributed vertical phase, $z=H\sin \psi$ has probability density function $[\pi \sqrt{H^2-z^2}{]}^{-1}$, so the midplane volume probability is ${\Sigma }_c/(\pi H)$.

In the genuinely three-dimensional regime $R_p\ll H\ll a_p$, the geometric collision cross-section is $\pi R_p^2$. The relative velocity differs from the coplanar value only at order $I^2$ for the specified high-speed encounters. Therefore the inclined comet collision time isAt fixed ${\rho }_p,M_{\star },I$ the leading scaling is $Q^{3/2}{M}_p^{-2/3}{a}_p^2$, with no leading $q$ dependence when $q\ll a_p$. The displayed factor retains the first dependence on $q/a_p$. Also $t_{\mathrm{col},\mathrm{I}}/t_{\mathrm{col}}\simeq 2a_{p}I/R_p$ in this averaging convention.

Inclination alone does not specify a collision time for one fixed orbit. If its orbital nodes miss the planet's orbit, a tilted Kepler orbit can have no collisions at all. The finite result is an orientation-averaged population result, or requires precession or scattering to sample those orientations. When $a_{p}I\lesssim R_p$, the three-dimensional estimate crosses over to the coplanar result instead of tending to zero.

In a planetary gravitational assist, the incoming and outgoing relative speeds have the same magnitude. If $\mathbf{u}$ is the planet-frame velocity, the heliocentric specific orbital energy change is consequentlynot an arbitrary addition of $|\Delta \mathbf{u}|$ to the heliocentric speed. The initial binding energy per unit mass is $G{M}_{\star }/(Q+q)\simeq G{M}_{\star }/Q$.

For a weak deflection, $b\gg b_{90}=G{M}_p/v_{\infty }^2$, the hyperbolic scattering formula gives $|\Delta \mathbf{u}|\simeq 2G{M}_p/(b{v}_{\infty })$. Its direction is perpendicular to the incoming relative velocity. At the planet, in the near-parabolic prograde limit,On the favorable side of the planet the energy gain is approximately $2G{M}_p{v}_p\kappa /(b{v}_{\infty })$. Equating this to the binding energy gives the ejection impact parameterOnly one sign of impact parameter gains energy for each radial branch. Noncolliding ejections have $R_p<|b|<b_{\mathrm{ej}}$, giving width $(b_{\mathrm{ej}}-R_p{)}_{+}$, where $x_{+}=\max (x,0)$. For $b_{\mathrm{ej}}\gg R_p$, the planetary ejection of a comet rate is approximately ${\Gamma }_{\mathrm{ej}}={\Sigma }_c{v}_{\infty }{b}_{\mathrm{ej}}$. With ${\mu }_p=M_p/M_{\star }$ this givesThe scaling applies when $b_{90}\ll b_{\mathrm{ej}}\ll r_H$ so that the local two-body encounter and weak-deflection approximations are consistent. The favorable-side coefficient assumes a coplanar orbit; averaging over other encounter geometries changes it. Cumulative diffusion by many subthreshold encounters is not included in this one-encounter estimate.

Exclude collisions from the ejection count. The effective collision width is $2R_p$, while the favorable noncolliding ejection width is $(b_{\mathrm{ej}}-R_p{)}_{+}$. Define the first flyby-ejection threshold and the equal-rate crossover byThere are no direct noncolliding ejections in this weak-deflection model for $Q\le Q_{\mathrm{hit}}$. Above that threshold, comparison of the two rates gives collision is more likely for $Q<Q_{\mathrm{crit}}$, while ejection is more likely for $Q>Q_{\mathrm{crit}}$. In the deep-crossing limit $Q_{\mathrm{crit}}\simeq [9/(2\sqrt{2})](M_{\star }/M_p)R_p$, so its scaling at fixed density is $M_{\star }{M}_p^{-2/3}{\rho }_p^{-1/3}$. Order-one factors depend on encounter geometry; the robust condition compares $(M_p/M_{\star })Q$ with $R_p$.

Adding the mutually exclusive rare-event loss rates gives the collision-ejection lifetime of a comet:For $Q\le Q_{\mathrm{hit}}$ this is $t_{\mathrm{col}}\propto Q^{3/2}$. For $Q\gg Q_{\mathrm{hit}}$ it approaches $t_{\mathrm{ej}}\propto Q^{1/2}$. A logarithmic sketch therefore bends from slope $3/2$ to slope $1/2$, with no maximum. At the equal-rate crossover the lifetime is half either individual loss time.

These asymptotic branches use the corrected residence density from part iii and require $Q\gg a_p$ and a local encounter scale below the Hill radius. They cannot be extrapolated indefinitely to arbitrarily distant weakly bound orbits.

For $0<e<1$, differentiate the logarithm of the proposed Poynting–Robertson drag invariant:Substitution of the orbit-averaged Poynting–Robertson drag rates gives, over the common denominator $a^2(1-e^2{)}^{3/2}$, the numerator $A[-(2+3e^2)+5e^2+2(1-e^2)]=0$. HenceEquivalently, eliminating time gives $d\log a/de=(2/5)(2+3e^2)/[e(1-e^2)]$, whose integral is $\log a=(4/5)\log e-\log (1-e^2)+\text{constant}$. An exactly circular orbit remains circular; the expression for $C$ is singular at $e=0$ and is replaced by that separate solution.

Using $q=a(1-e)$ and $Q=a(1+e)$, differentiate each apsidal distance. The apsidal evolution under Poynting–Robertson drag isBoth the pericentre distance and apocentre distance decrease. Since $a=(Q+q)/2$ and $e=(Q-q)/(Q+q)$, the same rates areTheir ratio isIn particular $Q\gg q$ gives $dQ/dq\simeq (5/3)(Q/q{)}^2$, whereas the circular limit gives slope $1$.

Parameterize the trajectory by the decreasing orbital eccentricity. The Poynting–Robertson drag invariant givesFor $Q_0\gg q_0$, $C\simeq 2q_0$. Initially $e\simeq 1$, and the orbit moves almost vertically down a plot of $Q/q_0$ against $q/q_0$: the apocentre distance shrinks rapidly while the pericentre distance changes little. Integrating the high-eccentricity slope givesOnce the orbit has moderate orbital eccentricity, both distances change appreciably. For example $e=1/2$ gives $q\simeq 0.766q_0$ and $Q\simeq 2.30q_0$ in the $Q_0/q_0\to \infty$ limit. Eventually $e\ll 1$, $q\simeq Q\simeq C{e}^{4/5}$, and the trajectory approaches the diagonal $Q=q$ before reaching the origin. Thus the two approximate phases are apocentre contraction at nearly fixed pericentre, followed by nearly circular inward migration. They are a smooth crossover, not two separate exact solutions.

For a circular Kepler orbit, $e=0$ remains an exact solution of the averaged Poynting–Robertson drag equations, and $\dot{a}=-2A/a$. Integrating gives $a^2(t)=a_0^2-4At$, so the point-star inspiral time under Poynting–Robertson drag isFor a star of radius $R_{\star }$, contact occurs instead at $(a_0^2-R_{\star }^2)/(4A)$ within this model. Sublimation or other forces can remove a real grain earlier.

Eliminate $a$ using the Poynting–Robertson drag invariant, $a=C{e}^{4/5}/(1-e^2)$. The eccentricity rate becomesSince the point-star limit corresponds to $a\to 0$ and $e\to 0$, the inspiral time under Poynting–Robertson drag isFor $e_0\to 0$, the integral is $(5/8)e_0^{8/5}+o(e_0^{8/5})$, recovering $a_0^2/(4A)$ despite the apparent singularity of $C$.

For $e_0\to 1$, use $C\simeq 2q_0$ and $1-e_0=2q_0/(Q_0+q_0)\simeq 2q_0/Q_0$. The endpoint asymptotic integral is $[2(1-e_0){]}^{-1/2}$, givingMost of this long lifetime is spent near the initially large apocentre distance, before substantial circularization. The result assumes weak drag and a valid orbital average; near-parabolic release can violate that assumption.

Assume negligible fragment release speed relative to a circular parent planetesimal. Its initial speed is $v_b=\sqrt{G{M}_{\star }/r_b}$, but radiation pressure reduces the grain's gravitational parameter to $(1-\beta )G{M}_{\star }$. Its specific orbital energy and specific angular momentum areFor $0\le \beta <1/2$, the resulting dust orbit released from a circular parent ring hasRelease is at periapsis: the inherited tangential speed exceeds the new circular speed. Low-$\beta$ grains remain near the dust birth ring; higher bound $\beta$ gives increasingly eccentric orbits extending well outside it. Uniform release longitude produces an approximately axisymmetric dust halo, with long residence times near apoapsis. At $\beta =1/2$ the conservative release orbit is parabolic; above this radiation-pressure blowout threshold it is initially unbound. For $\beta \ge 1$, zero or repulsive effective gravity also prevents a bound ellipse.

The left panel shows one release longitude; rotating it around the ring fills the halo. At longer times Poynting–Robertson drag also carries bound grains inward, so an initially outward-spanning release distribution does not imply an empty inner region in a continuously supplied, collision-free population.

For the bound release orbits, use $e_0=\beta /(1-\beta )$ andThe inspiral time under Poynting–Robertson drag is thenThe limiting branches areThus the formal bound-grain curve, shown in the right panel of the preceding figure, diverges at both ends and has a minimum at an intermediate $\beta$. At $\beta =0$ there is no Poynting–Robertson drag and the circular grain survives indefinitely in this collision-free model.

Initially unbound grains do not have a bound-orbit inspiral lifetime. For $\beta \ge 1/2$ the release calculation gives outward trajectories; their residence time in a specified finite observation region is a dynamical escape time, not the averaged $t_{\mathrm{PR}}$ above. A numerical escape lifetime cannot be assigned without specifying the outer boundary. Extremely close to $1/2$, energy lost on the first passage can alter even this bound/unbound classification: orbital averaging requires roughly $1-2\beta \gg v_b/c$. Accordingly the divergence is a formal limit of the supplied secular equations, not an exact prediction for marginal grains.

Differentiating the disturbing function in Lagrange planetary equations givesFor the complex eccentricity $z=e{e}^{i\varpi }$,Use $\sin \delta +i\cos \delta =i{e}^{-i\delta }$. Each planetary term becomes $i{A}_j{e}_j{e}^{i{\varpi }_j}=i{A}_j{z}_j$, soUnlike the polar orbital element equations, this linear equation is nonsingular at $e=0$.

The planetary expansion is a decomposition into secular eigenmodes of conservative linear Laplace-Lagrange secular theory. It assumes small orbital eccentricities, coplanarity, approximately constant semi-major axis values, and averaging over fast orbital phases away from important mean-motion resonances and close encounters. Dissipation, migration, and nonlinear secular effects are absent from this constant-coefficient representation.

The $g_i$ are real eigenfrequencies of the Laplace-Lagrange secular matrix, describing apsidal precession. For mode $i$, the signed real coefficients $e_{ji}$ give its amplitude in planet $j$; their relative magnitudes are fixed by an eigenvector, and their signs distinguish aligned and anti-aligned apsides. The common phase ${\beta }_i$ specifies that mode's orientation at the chosen time origin. A planet's actual orbital eccentricity is $|z_j(t)|$, which generally varies through interference of modes; it is not a sum of the positive mode amplitudes.

Put $B_i={\sum }_j{A}_j{e}_{ji}$. For $A\ne g_i$, the solution of the linear secular forcing of a test particle equation isThe initial condition gives $e_p=|z(0)-z_f(0)|$ and ${\beta }_p=\arg [z(0)-z_f(0)]$. The proper eccentricity $e_p$ is the constant amplitude of the homogeneous response; its proper longitude of periapsis is $At+{\beta }_p$. The forced eccentricity is the driven vector $z_f(t)$, not necessarily a constant magnitude.

On an Argand diagram, $z-z_f$ describes a circle of radius $e_p$. Its center $z_f$ can itself move under the planetary modes, so $z$ need not trace one fixed circle in the inertial complex plane. At a secular resonance $A=g_i$, that mode instead produces $i{B}_{i}t{e}^{i(At+{\beta }_i)}$, and the undamped linear response grows until the approximation fails.

For one isolated perturbing planet, its Kepler orbit has fixed longitude of periapsis, so $z_1$ is constant and its secular forcing frequency is zero. A constant forced eccentricity satisfies $0=iA{z}_f+i{A}_1{z}_1$. Since $A=-A_{1}f({\alpha }_1)$,This assumes $A\ne 0$. If an external process made the planet precess at frequency $g$, the particular solution would instead be $z_f=A_1{z}_1/(g-A)$; the stated formula relies on the isolated single-planet case.

To first order in orbital eccentricity, the semi-minor axis is $a+O(e^2)$, while the ellipse's center is displaced from the star by $-a\mathbf{e}$. Thus each orbit is a circle of radius $a$ about $-az$, in complex spatial coordinates.

Write $z=z_f+e_p{e}^{i{\varpi }_p}$ at one secular epoch. The centers of the individual circles lie on a circle of radius $a{e}_p$ about $C_f=-a{z}_f$. Their union is the proper-eccentricity annulusIts full radial width is $2a{e}_p$ and its center is displaced by $a{e}_f$ toward the forced apoapsis direction. Relative to the star, if the forced longitude of periapsis is zero,The middle panel of the preceding figure shows this construction. For fully sampled orbital phases and proper phases, the annulus is occupied throughout, but its density need not be uniform; turning points in the proper radial excursion give integrable edge enhancements. The construction requires $e_f,e_p\ll 1$ and is an instantaneous secular geometry.

If the proper eccentricity is negligible and the semi-major axis ranges over $a\pm \Delta a$, particles lie on nearly aligned forced Kepler orbits. For a locally constant forced vector, these are geometrically similar ellipses, represented to first order by circles of radius $a^{\prime }$ centered at $-a^{\prime }{\mathbf{e}}_f$.

The boundaries of the aligned eccentric ring are thereforeIt is narrower at forced periapsis and wider at forced apoapsis, unlike the constant-width proper-eccentricity annulus. Its centers shift slightly between the two edges, as shown in the right panel of the figure. The condition $e_p/e_f\ll \Delta a/a$ makes the proper radial excursion negligible compared with the variation in forced center location.

The forced vector need not actually be constant across the interval. More generallyFor aligned apsides this becomes $2\Delta a[1-(e_f+a{e}_f^{\prime })\cos \theta ]$, provided the orbits remain nested. The radial gradient of forced eccentricity matters for the density comparison even in a narrow ring.

For a phase-mixed orbit, probability in a short arc is $dt/T=ds/(Tv)$, so the line density on a Kepler orbit is inversely proportional to speed. At true anomaly $\theta$,Using the vis-viva equation, $v^2=(G{M}_{\star }/a)[1+2e\cos \theta ]+O(e^2)$ and hence $v=na[1+e\cos \theta ]+O(e^2)$. ThereforeFor one particle, the normalized line probability is $(2\pi a{)}^{-1}(1-e\cos \theta )+O(e^2)$. This is density per arc length; density per unit polar angle has a further factor $ds/d\theta \simeq a(1-e\cos \theta )$ and is proportional to $1-2e\cos \theta$.

In the proper-eccentricity annulus, averaging over random proper phases cancels the proper vector in the first-order angular variation. The mean line density is proportional to $1-e_f\cos \theta$, while the full ring width is independent of angle. Thus at corresponding radial locationsThere is a forced-apocentre enhancement in number surface density of a disk, together with the proper-excursion edge profile.

For a narrow aligned eccentric ring with locally constant $e_f$, its width has the same first-order factor $1-e_f\cos \theta$ as its line density. Their ratio is constant: the constant-eccentricity aligned ring has no first-order surface-density contrast. It still contains more particles per unit longitude at apoapsis, spread across a wider area.

The more general surface density contrast in an eccentric ring depends on the eccentricity gradient. With a smooth number distribution per unit semi-major axis and aligned apsides,Hence decreasing $e_f(a)$ can also make the aligned ring apocentre-dense; for $e_f\propto a^{-1}$ it has the same first-order contrast as the proper-phase case. A positive gradient gives a pericentre enhancement. The problem leaves this profile unspecified, so there is no universal observational distinction based solely on the two apse densities.

Resolved ring width and a measured forced-eccentricity profile can distinguish the constant-width and similar-ellipse constructions. Interpreting observed brightness additionally requires the grain temperature, illumination, emissivity, and viewing geometry; brightness is not itself particle number density. Inferences from one contrast alone can therefore be degenerate.

Measure angles from periapsis in the direction of motion. The true anomaly $f$ is the angle at the stellar focus between the pericentre ray and the star-particle radius vector. Draw the auxiliary circle of radius $a$ around the ellipse's center; the point on that circle with the same projected major-axis coordinate as the particle defines the eccentric anomaly $E$ at the center.

The mean anomaly $M=n(t-\tau )$ is a uniformly advancing time phase, with $M=0$ at pericentre passage $t=\tau$. It is the angle swept by a fictitious uniformly moving circular phase marker, not generally the actual position angle. The mean motion is the angular frequencyThe second formula is Kepler's equation, which converts the uniform time phase into geometric position.

With the ellipse center as origin and its major axis horizontal, the auxiliary-circle construction gives particle coordinates $X=a\cos E$, $Y=a\sqrt{1-e^2}\sin E$, while the stellar focus is at $(ae,0)$. ThereforeSquaring and adding gives $r^2=a^2(1-e\cos E{)}^2$. The Kepler orbit radius and the anomaly conversion are consequentlyEquivalently,The sine and cosine formulas, or a quadrant-aware angle function, fix the correct branch when using the half-angle relation.

Iterating Kepler's equation once gives $E=M+e\sin M+O(e^2)$. The exact radius expression then yieldsExpansion of the true anomaly conversion gives $f=E+e\sin E+O(e^2)$, soThe factor of two includes both the nonuniform relation between $M$ and $E$ and the change from center angle $E$ to focus angle $f$.

Choose the common initial longitude of the particle and $G$ as zero. Let the fixed inertial longitude of periapsis be $\varpi$, positive ahead of that ray. At release $f_0=-\varpi$, but the mean anomaly isThereafter $M(t)=nt+M_0$, $f(t)=M(t)+2e\sin M(t)+O(e^2)$, and the particle's longitude is $\lambda (t)=\varpi +f(t)$. The reference ray has longitude ${\lambda }_G=n_{g}t$.

In the rotating reference frame, the radial $x$ axis points outward along the star-$G$ ray and the $y$ axis points forward. Define $\delta \lambda =\lambda -{\lambda }_G$. The exact positional relations are

The diagram shows a later instant; initially $\delta \lambda =0$, so $y(0)=0$. Retaining the $O(e)$ correction in $M_0$ is what will preserve that initial condition in the approximate trajectory.

To first order, the initial-phase correction in $M_0$ can be omitted inside a term already multiplied by $e$. ThusFor $|a-a_g|/a\ll 1$ and $|\delta \lambda |\ll 1$, expand the exact rotating reference frame coordinates. Products of the small radial and angular displacements are second order. This givesAt $t=0$ the two sine terms cancel, as required. Equal mean motion gives a closed local epicyclic motion with radial-to-tangential amplitude ratio $1:2$; unequal mean motion adds along-track drift. The local expansion also requires $|(n-n_g)t|\ll 1$, so small frequency mismatch alone does not justify it for arbitrarily large times.