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

Expand the secular Poynting–Robertson drag rates using $(1-e^2{)}^{-3/2}=1+\frac{3}{2}{e}^2+O(e^4)$ and $(1-e^2{)}^{-1/2}=1+O(e^2)$. Through second order in orbital eccentricity,The second formula follows from $d(e^2)/dt=2e\dot{e}$. The cubic correction in $\dot{e}$ contributes only at fourth order in $d(e^2)/dt$.

In resonant trapping of dust, the resonant torque must compensate the inward Poynting–Robertson drag. At the nearly fixed resonant semi-major axis, ${\dot{a}}_{\mathrm{res}}+{\dot{a}}_{\mathrm{pr}}=0$ requiresSince $|\sin \varphi |\le 1$, a necessary torque-balance condition is $|C|e^q\ge A(1+3e^2)/[a^2(p+q)]$. When the limiting orbital eccentricity is small, this givesEquality allows only the extreme phase $|\sin \varphi |=1$, with no margin for restoring torque; a robust trapped state needs $e>e_{\min }$. This is a necessary condition for resonance capture, not a sufficient probability-of-capture criterion for a migrating particle.

Write $u=e^2$ and use the torque balance from part (ii) to eliminate the resonant argument. The resonant contribution isAdding the Poynting–Robertson drag contribution gives the dissipative resonant eccentricity equilibrium equationwhereFor $5p>q$, integrating with $u(0)=e_0^2$ givesHere $a$ is the fixed resonant semi-major axis and the expression describes averaged evolution, neglecting the small resonant-argument libration. The physical bound-orbit interpretation also requires $B<1$; a quantitatively controlled second-order calculation needs $B\ll 1$.

At $e_{\max }=\sqrt{B}$, torque balance givesThe magnitude of the right-hand side must not exceed one. This fixes the dissipative offset of the resonant-argument libration centre.

For $C>0$, put $S=5A/[C{a}^2(5p-q)B^{q/2}]$. The two torque-balanced phases are ${\varphi }_0=-\arcsin S$ and ${\varphi }_0=\pi +\arcsin S$, modulo $2\pi$. In the usual weak-perturbation pendulum approximation, $\delta \dot{\varphi }\simeq -3(p+q)n\delta a/(2a)$ and $\delta \dot{a}\simeq -2(p+q)Ca{B}^{q/2}\cos {\varphi }_0\delta \varphi$. Restoring libration thus selects $\cos {\varphi }_0<0$, givingIt approaches $\pi$ when drag is weak compared with resonant torque. For $C<0$, the corresponding restoring branch has $\cos {\varphi }_0>0$ and approaches zero. The sine equation is the direct result independent of this leading-order stability approximation.

At periapsis, the particle's mean longitude equals its longitude of periapsis: $\lambda =\varpi$. Hence the pericentre geometry of an exterior resonance givesThe resonant argument therefore fixes the angular separation from the planet whenever the particle reaches periapsis. With $r_{\mathrm{peri}}=a(1-e)$ and planetary radius $a_p$, the actual separation isThus a restricted libration phase can prevent close encounters even when the radial ranges of the two orbits overlap: this is resonance protection.

The specific angular momentum is $h=\sqrt{G{M}_{\star }a(1-e^2)}$, so the periapsis angular speed isAt the exterior mean-motion resonance, neglecting radiation pressure, $n/n_p=p/(p+q)$. Thus at the limiting orbital eccentricity $e=\sqrt{2q/(5p-q)}$,The condition for overtaking the planet in angle isFor the formal second-order equilibrium, it holds for every physical bound value $0<q<5p/3$. To see this, rewrite $q/p=5e^2/(2+e^2)$, so $p/(p+q)=(2+e^2)/(2+6e^2)$. For $0<e<1$, $\sqrt{1+e}/(1-e{)}^{3/2}>1+2e$, whileThe instantaneous angular motion is therefore prograde relative to the planet near periapsis, although its mean relative motion is retrograde. The range $q\ll p$ is the regime where the eccentricity expansion is most reliable.

For the $5:3$ exterior mean-motion resonance, $p=3$, $q=2$. The truncated dissipative resonant eccentricity equilibrium givesFor weak drag and $C>0$, take the resonant-argument libration centre ${\varphi }_0\simeq \pi$. The pericentre geometry of an exterior resonance then places the three pericentres at angles $\pi /3$, $\pi$, and $5\pi /3$ relative to the planet. The trajectory closes after three particle orbits and five planet orbits. It moves retrograde in the rotating frame on average but develops a prograde loop near each periapsis, because the periapsis angular speed exceeds the planet's angular speed.

For a finite drag offset ${\varphi }_0=\pi +\arcsin S$, rotate the whole particle pattern through $\arcsin S/3$. The sketch uses the eccentricity requested by the second-order model; $e\simeq 0.555$ is not small, so the curve is an illustration of that approximation rather than a quantitatively accurate integration of the full dissipative resonance.

Assume spherical bodies of common mass density $\rho$, a well-mixed belt, and that $n(D)dD$ counts bodies in the entire belt rather than per unit volume. A body of diameter $D$ has mass $\pi \rho D^3/6$. The mass normalization of a power-law size distribution isConsequentlyThe thin-torus volume is circumference times radial width times full vertical height,The local number-density distribution needed for a planetesimal collision rate is $n(D)/V$. Thus belt geometry does not enter this total-number $K$ but does enter collision rates. This approximation treats the stated height as uniform, independent of size.

Let $d$ be an impactor diameter and assume it is much smaller than target diameter $D$, with equal material mass density. Neglect gravitational focusing and take the encounter speed to be the stated target-size law $v_{\mathrm{rel}}=v_p{D}^p$. The specific impact energy isA catastrophic planetesimal collision requires $Q\ge Q_D^{\ast }=Q_b{D}^b$. Therefore the catastrophic projectile diameter isUse $\max (D_{\mathrm{cc}},D_{\min })$ as the actual lower size cutoff when integrating a rate. If this exceeds $D_{\max }$, the specified population contains no destructive impactors. The small-impactor approximation itself requires $D_{\mathrm{cc}}\ll D$.

In the gravitational focusing regime, take $d\ll D$ and $v_{\mathrm{rel}}\ll v_{\mathrm{esc}}$ at large separation. Mutual gravitational acceleration then makes the impact speed approximately the target's escape velocity:Setting the specific impact energy equal to the catastrophic disruption threshold givesThe size-dependent far-field encounter speed drops out of this leading threshold because gravitational acceleration dominates at impact. Gravitational focusing also increases the collision cross-section by approximately $1+v_{\mathrm{esc}}^2/v_{\mathrm{rel}}^2$, but this concerns the rate rather than the energy threshold. These formulas assume the disruption-energy convention remains per unit target mass and that gravitationally accelerated impact energy couples to disruption as specified by $Q_D^{\ast }$.

Neglect gravitational focusing and assume a fixed well-mixed volume $V$, $\alpha >1$, and $D_{\min }\ll D_{\mathrm{cc}}\ll D\ll D_{\max }$. For the small impactors that dominate the rate, the geometric collision cross-section is approximately $\pi D^2/4$. ThusEvaluating and neglecting the upper-end term gives the catastrophic planetesimal collision rateUsing part (ii), its size dependence isIf the destructive cutoff is near an endpoint or comparable to $D$, retain the full integral and the cross-section $\pi (D+d{)}^2/4$. The simple power law then need not apply.

For equal logarithmic width $d\log D$, the number of bodies is $n(D)Dd\log D$. Their mass is thereforeThe mass removed by catastrophic planetesimal collisions per logarithmic bin per unit time scales asAssume a steady collisional cascade with supply only at the large end, loss only at the small end, negligible reassembly or cratering contributions, and self-similar fragmentation that transfers mass through comparable logarithmic intervals. Mass conservation requires the downward mass flux across every intermediate size boundary to be the same; otherwise mass would accumulate or deplete there. Under the scale-invariant bin-transfer assumption, this is proportional by a size-independent factor to ${\dot{M}}_{\mathrm{loss}}$, soThis is constant mass flux in a collisional cascade. Steady state alone fixes the net flux, while the identification with the destruction rate of each bin additionally uses the stated fragmentation assumptions.

Combining constant mass flux in a collisional cascade with the rate in part (iv) givesSet the exponent to zero and multiply by three:or $21+b+p-\alpha (6+b-2p)=0$. Therefore the velocity-dependent collisional cascade hasFor $b=p=0$, this recovers the Dohnanyi collisional cascade slope $\alpha =7/2$. The denominator must be nonzero, and the resulting slope and destructive cutoff must satisfy the assumptions used for the collision-rate integral. A size-dependent vertical belt thickness would supply an additional volume factor and change this result.

Let $\mathbf{r}={\mathbf{r}}_2-{\mathbf{r}}_1$ and $r=|\mathbf{r}|$. Newton's second law for the two planetesimals givesSubtracting yields the forced relative two-body problemThe internal gravitational forces cancel from the centre of mass equation but add in the relative equation; external forces affect the internal orbit through their difference in acceleration.

Put $m=m_1+m_2$, $\mu =Gm$, $n=\sqrt{\mu /a^3}$, $v=na$, and $\theta =nt$. Work instantaneously in the inertial frame in which the centre of mass is at rest. The unperturbed circular relative orbit hasThe individual velocities are ${\mathbf{v}}_1=-(m_2/m)\mathbf{v}$ and ${\mathbf{v}}_2=(m_1/m)\mathbf{v}$. Linear gas drag gives ${\mathbf{F}}_i/m_i=-({\mathbf{v}}_i-v_g\hat{\mathbf{y}})/t_i$, where $t_i=t_{si}$ are the aerodynamic stopping times. DefineThen the gas drag on a planetesimal binary isEquivalently, $\mathbf{f}=-\kappa \mathbf{v}+\Delta v_g\hat{\mathbf{y}}$. If the centre of mass moves, $v_g\hat{\mathbf{y}}$ here is its relative gas velocity, held approximately fixed during one short binary orbit.

With $\hat{\mathbf{r}}=(\cos \theta ,\sin \theta )$ and $\hat{\mathbf{\theta }}=(-\sin \theta ,\cos \theta )$, project the relative gas drag to obtainThe specific orbital energy is $\epsilon =v^2/2-\mu /r=-\mu /(2a)$. Differentiation cancels the conservative gravitational work and gives $\dot{\epsilon }=\mathbf{v}\cdot \mathbf{f}$. For the circular orbit,Using $v^2=\mu /a$, the orbit-averaged drag work relation therefore gives the instantaneous leading-order changeThe differential headwind term can change sign during the orbit; the damping term always removes orbital energy.

For weak gas drag, evaluate the work on the unperturbed circular orbit and average over one orbital period. The gas velocity changes negligibly during that period and $⟨\cos \theta ⟩=0$. Hence gas-driven contraction of a planetesimal binary obeysWith constant aerodynamic stopping times,This is the secular, orbit-averaged separation. It assumes the perturbations are small enough for the orbit to remain nearly circular; a large differential headwind can excite eccentricity or disrupt a weakly bound pair instead.

Let $\mathbf{V}$ be the centre of mass velocity and $\mathbf{v}$ the binary relative velocity. Summing the two linear gas drag forces givesThe last term averages to zero over the internal circular orbit. On the nearly circular stellar orbit, $\mathbf{V}\simeq v_K\hat{\mathbf{\theta }}$ and ${\mathbf{v}}_g=(1-\eta )v_K\hat{\mathbf{\theta }}$, so the averaged tangential acceleration is $-\eta {\kappa }_b{v}_K$.

With specific orbital energy $-G{M}_{\star }/(2a_b)$, the orbit-averaged drag work is $-\eta {\kappa }_b{v}_K^2$, where $v_K^2=G{M}_{\star }/a_b$. Thus gas-driven migration of a binary centre of mass givesThis uses the weak-drag, nearly Keplerian approximation implicit in applying the circular-orbit work relation. Strongly gas-coupled orbits require a coupled radial-azimuthal drift solution.

Put $ϵ=m_2/m_1\ll 1$ and $t_{s2}/t_{s1}={ϵ}^k$, with $0<k<1$. The gas-driven contraction of a planetesimal binary and gas-driven migration of a binary centre of mass times obeyThereforeThe low-mass secondary has the shorter aerodynamic stopping time and chiefly controls internal damping, whereas the primary chiefly controls centre of mass drag. The pair can contract substantially, potentially merging, before its stellar orbital radius changes appreciably. This conclusion is restricted to a pair that remains bound and within the adopted linear-drag regime.

The variables $y_j=I_j{e}^{i{\Omega }_j}$ are complex inclinations, despite the question's description as eccentricities. Diagonalize the Laplace-Lagrange inclination matrix: let real eigenvectors ${\mathbf{v}}_k$ satisfy $B{\mathbf{v}}_k={\lambda }_k{\mathbf{v}}_k$, and form $V=({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3)$. The physical matrix is diagonalizable by positive diagonal symmetrization of a matrix, as shown below.

The matrix exponential solution isWrite $c_k=I_k{e}^{i{\gamma }_k}$ with $I_k\ge 0$, and define $I_{jk}=I_k{v}_{jk}$. This givesThe eigenvalues determine nodal precession frequencies, the eigenvectors determine relative inclinations within each secular eigenmode, and the initial complex inclinations determine their amplitudes and phases. The real $I_{jk}$ can have either sign, encoding opposite nodal phases. The zero-frequency mode is a common fixed tilt.

Let $a=B_{12}$, $b=B_{13}$, $c=B_{21}$, $d=B_{23}$, $e=B_{31}$, and $f=B_{32}$. Zero row sums put $B$ in the form used in the spectrum of a three-dimensional matrix with zero row sums. In particular, $B(1,1,1{)}^T=0$. DefineExpanding the characteristic polynomial gives $\det (\lambda I-B)=\lambda ({\lambda }^2+S\lambda +T)$, soFor the physical coefficients, $L_j{B}_{jk}=L_k{B}_{kj}$ with $L_j=M_j\sqrt{G{M}_{\star }{a}_j}$. Thus $W=\mathrm{diag}(L_j)$ gives positive diagonal symmetrization of a matrix, andConsequently the eigenvalues are real and nonpositive, and $B$ has a complete eigenvector basis. This justifies the mode expansion in part (i).

For small $\alpha$, expand the defining integrand of the Laplace coefficient:HenceThe constant term integrates to zero and ${\int }_0^{2\pi }{\cos }^{2}xdx=\pi$, givingThere is no quadratic term: its factors are combinations of $\cos x$ and ${\cos }^{3}x$, which both integrate to zero over a period. The linear approximation is used for the distant planet, not for the closely spaced inner pair.

Keep the inner-pair Laplace coefficient unexpanded because $\alpha =a_1/a_2$ is close to one. Put $\ell =L_1/L_2$ andUsing the small-argument Laplace coefficient for the distant planet and $n_2/n_1={\alpha }^{3/2}$ givesThe reverse-coupling expressions follow from the angular momentum identity $L_j{B}_{jk}=L_k{B}_{kj}$. If desired, $L_1/L_2=(M_1/M_2)\sqrt{\alpha }$, $L_1/L_3=(M_1/M_3)\sqrt{\beta }$, and $L_2/L_3=(M_2/M_3)\sqrt{\beta /\alpha }$ express all ratios in the supplied quantities. The diagonal entries are $B_{11}=-(B_{12}+B_{13})$, $B_{22}=-(B_{21}+B_{23})$, and $B_{33}=-(B_{31}+B_{32})$.

For comparable masses and closely spaced inner orbits, $\alpha \simeq 1$, $\ell \sim 1$, and $L_1/L_3,L_2/L_3=O(\sqrt{\beta })$. ThusThe distant couplings $B_{13},B_{23}$ scale as ${\beta }^3$ relative to a fixed inner architecture, and the reverse couplings have an additional $\sqrt{\beta }$ suppression.

Ignoring the distant planet at first, the inner two-by-two block has eigenvalues zero and $-(B_{12}+B_{21})$. Therefore its fast differential secular eigenmode hasFor the slow mode, strong inner coupling makes $y_1\simeq y_2=y_{\mathrm{in}}$. Add the first two evolution equations weighted by angular momenta and divide by $L_{12}=L_1+L_2$. The rigid inner pair in hierarchical secular inclination dynamics has equationswithThe two-plane matrix has eigenvalues zero and $-c(1+L_{12}/L_3)$. The architecture implies $L_3\gg L_{12}$ for comparable masses, so its leading slow frequency isThe neglected factor $1+L_{12}/L_3$ gives a relative correction of order $\sqrt{\beta }$; imperfect locking of the inner pair supplies higher-order corrections in its weak coupling to the outer planet.

Initially the two inner orbital normals coincide, so their combined angular momentum is $L_{12}=L_1+L_2$ along one direction. Use the outer orbital normal as the vertical axis. The inner pair's vector then has perpendicular component $L_{12}\sin \Delta I$ and parallel component $L_{12}\cos \Delta I$. Adding the outer angular momentum givesThe invariable plane is perpendicular to this total vector. Thus the outer orbit's orbital inclination to that plane equals the angle between its own normal and the total normal. The inclination to the invariable plane of two orbital planes isThis uses the small, prograde inclinations for which the denominator is positive; in general use the corresponding two-argument arctangent to choose the quadrant. If $A\gg 1$, $I_3\simeq \Delta I/(A+1)$ for small $\Delta I$: the distant planet carries most of the angular momentum and lies close to the invariable plane.