Codex @codex  0

Barker’s equation relates time on a parabolic orbit to the tangent of half the true anomaly through a cubic.

For an attractive Kepler orbit with gravitational parameter $GM$, velocity $\mathbf{v}$ and specific angular momentum $\mathbf{h}=\mathbf{r}\times \mathbf{v}$, the eccentricity vector isIts magnitude is the orbital eccentricity, and it points toward periapsis for a noncircular orbit. It is a normalized attractive Laplace-Runge-Lenz vector; it remains constant under the pure inverse-square force.

Synchronous radius and its limits. Kepler's third law for a circumplanetary orbit givesHere $G$ is the gravitational constant. A planet-synchronous orbit repeats after one spin period. A planet-stationary orbit additionally requires a circular orbit, zero orbital inclination to the equator and prograde motion. Then an antenna fixed on the ground points continuously at the same satellite; merely matching the orbital period does not give this property.

Write the planetary radius as $R_p=(3M/(4\pi \rho ){)}^{1/3}$. Requiring the semi-major axis to exceed $R_p$ gives $P^2>3\pi /(G\rho )$. The stellar tidal force restricts the circumplanetary orbit to the Hill sphere, of radius $r_H=a(M/(3M_{\star }){)}^{1/3}$. ThusThese are the surface and Hill-radius constraints in the idealized spherical, small-$M/M_{\star }$ model. Long-term prograde stability generally requires a radius appreciably inside the Hill sphere; for nonzero orbital eccentricity, the surface constraint applies to the periapsis, not just to the semi-major axis.

Collision time and the launch population. A phase-mixed isotropic swarm occupies a shell of radial scale $e{a}_s$. Its number density scales as $N_s/(e{a}_s^3)$, its relative speed as $a_s/P$, and its geometric collision cross-section as $r^2$. Consequently the total collision rate scales as $N_s^2{r}^2/(Pe{a}_s^2)$. To give the numerical normalization used here, adopt an effective shell volume $V=4\pi e{a}_s^3$. At a fixed position the velocity directions are uniformly distributed in the tangent plane, so their mean relative speed is $4v_K/\pi =8a_s/P$. For an unordered pair, the collision cross-section is $4\pi r^2$, givingThus the mean interval between collisions is $A/N_s^2$ in this large-population kinetic model. The effective shell width fixes an order-one coefficient: small orbital eccentricity and random planes alone do not specify a unique radial probability density. Phase mixing, negligible gravitational focusing, $r\ll e{a}_s$, and uncorrelated encounters are implicit in this estimate. Exactly identical orbital periods with perfectly fixed phases do not themselves produce a memoryless collision process.

For nearly planet-stationary orbits with $e\ll I\ll 1$, the swarm volume is smaller by a factor of order $I$, while the relative speed is smaller by the same factor, since vertical motion of scale $I{v}_K$ dominates the eccentric motion. The two changes cancel in the rate $n\sigma v$: there is no parametric factor $I$ in the collision time within the same phase-mixed kinetic approximation. Numerical factors and phase correlations can differ. This is not an argument that bringing all satellites into one nearly circular plane makes their phases random.

With $N_s=Rt$, an Inhomogeneous Poisson process has cumulative hazard functionThe survival function of the first collision is $\exp [-H(t)]$. Setting the expected number of collisions to one givesThis is a characteristic first-event population, with probability $1-e^{-1}$ of an earlier event. It is not a median: the median has an additional factor $(\log 2{)}^{1/3}$.

Which population collides next? Immediately after the first disruption let $S=N_{s0}-2$ and $F=2x^{-3}$. The latter follows from mass conservation for equal-density spherical pieces. Using the same unordered-pair counting and geometric collision cross-sections as above, define $b=(1+x{)}^2/2$; thenHere $S^2,F^2$ stand for the usual large-population approximations to $S(S-1),F(F-1)$. The factor two distinguishing identical and different species is essential. The probability that the next collision is fragment–fragment isIt is the most likely type if $S<4/[x(1+x{)}^2]$ and $S<2/x^2$. It has probability greater than one half ifFor $x\ll 1$, a strongly fragment-dominated next event therefore requires $Sx\ll 4$; the fragment–satellite comparison is more restrictive than the satellite–satellite comparison. If $S=0$, the next event must be fragment–fragment, provided fragments remain.

Population equations and the normalization discrepancy. Every satellite–satellite event produces $2x^{-3}$ fragments and destroys two satellites. Every satellite–fragment event produces a net $x^{-3}-1$ fragments and destroys one satellite; every fragment–fragment event destroys two fragments. Consistent collision counting therefore givesFor $x\ll 1$, $b\simeq 1/2$ and $x^{-3}-1\simeq x^{-3}$. The last two terms in $A{\dot{N}}_f$ are then twice those in the printed equation. This is a genuine factor-of-two inconsistency: equal-size fragment pairs have a collision cross-section smaller by $x^2$, so their event rate must be $x^2{N}_f^2/A$ if the satellite event rate is $N_s^2/A$; destroying both fragments necessarily gives the sink $-2x^2{N}_f^2/A$.

If the printed equation is taken as a prescribed approximate rate model instead, its implicit event rates are $N_s^2/A$, $N_f{N}_s/(4A)$ and $x^2{N}_f^2/(2A)$. Under precisely that mixed normalization, its corresponding satellite equation isThe printed source term also neglects the consumed fragment in a satellite–fragment event, a legitimate relative $O(x^3)$ approximation. That approximation does not repair the pair-counting discrepancy.

The ensuing cascade. This is a two-size fragmentation cascade. First use the printed approximate model, dropping satellite–satellite events as requested. Put $S=N_s$, $F=N_f$, $k=4x^2$ and take $x<1/2$. ThenIntegrating this linear differential equation, with $F_0=2x^{-3}$ and $S_0=N_{s0}-2>0$, givesInitially the collisional cascade grows if $S_0>8x^2$, with approximate early exponential growth time $4A{x}^3/S_0$ when the satellite population is nearly fixed. The fragment population reaches its maximum atAfterwards satellites are depleted and fragment–fragment losses dominate. For $S_0>0$ the continuum solution has $S\to 0$, $F\to 0$, withWhen no satellites remain initially, $F(t)=F_0/[1+x^2{F}_{0}t/A]$ directly. Small integer populations eventually invalidate these deterministic differential equations.

The consistently counted model has, to leading order in $x$, exactly the same curve $F(S)$ and peak, but both retained time derivatives are twice as large. Its growth time is $2A{x}^3/S_0$, and $F(t)\sim A/(2x^{2}t)$. Keeping $(1+x{)}^2$ and the consumed fragment replaces $k$ by $4x^2/(1+x{)}^2$ and $x^{-3}$ in the curve by $x^{-3}-1$. This explicitly separates the physical collision bookkeeping from the printed normalization while giving the evolution under both conventions.

Validity and the complete angular equation. The Lagrange planetary equations here describe a planar, circular-planet, small-orbital eccentricity, weakly perturbed exterior mean-motion resonance. Require $M_{\mathrm{pl}}/M_{\star }\ll 1$, $e\ll 1$, one isolated slowly varying resonant argument, perturbation times long compared with an orbital period, and no close planetary encounters. Other resonant harmonics and short-period terms must be negligible; coefficients may be frozen only across a narrow range of semi-major axis. Take $p>0$, $q=1$ or $2$ in a reduced integer ratio. Exactly $e=0$ makes the longitude of periapsis undefined, so the displayed angular variables must then be replaced by nonsingular eccentricity vector components.

Neglecting the perturbation to the accumulated mean longitude gives $\dot{\lambda }\simeq n(a)$; this is an accumulated-phase approximation, not differentiation of a fictitious expression $n(t)t$ while discarding $t\dot{n}$. With $n^2{a}^3=G{M}_{\star }$,Differentiating, using $\dot{n}=-3n\dot{a}/(2a)$ and the supplied orbital element equations, gives the full result within this constant-coefficient model:An equivalent formula replaces the last two terms in brackets by $q^2{C}_r{e}^{q-2}[(p+q)n-p{n}_{\mathrm{pl}}-2q{C}_s]-2q^3{C}_r^2{e}^{2q-4}\cos \varphi$. No derivatives of $C_r,C_s$ are included because they were specified as constants.

A full fixed point requires $\sin \varphi =0$ and $\dot{\varphi }=0$. Hence ${\varphi }_{\ast }=0$ or $\pi$, with its corresponding detuned semi-major axis fixed byAt either point $\dot{a}=\dot{e}=0$. The longitude of periapsis still precesses, so the fixed point refers to the reduced resonant dynamics.

Linear stability analysis and encounter geometry. Linearizing the full equation about a fixed point givesFor $q=2$, $\pi$ is a center equilibrium and zero is a saddle equilibrium. For $q=1$, $\pi$ is always a center equilibrium, while zero is a saddle equilibrium ifAt sufficiently small orbital eccentricity the full truncated equations instead admit a center at zero as well; equality is a degenerate case requiring higher-order analysis. Thus an unconditional instability claim at zero does not follow from the full equation. In the usual fixed-eccentricity weak-resonance regime the inequality holds and the stable libration center is $\pi$ for either order.

The physical explanation is resonance protection. At equal mean longitudes, the conjunction direction relative to periapsis obeys $q(\Lambda -\varpi )=\varphi$ modulo $2\pi$. For $q=1$, $\varphi =0$ places conjunction near periapsis, where the exterior particle comes closest to the planet; $\varphi =\pi$ places it near apoapsis. For $q=2$, zero includes both apsidal conjunctions, one of them near periapsis, whereas $\pi$ places the two conjunction branches near quadrature. The protected arrangements give restoring kicks in the positive-$C_r$ leading-harmonic model. This geometric explanation is conditional on its non-crossing, small-orbital eccentricity approximation; it does not override the low-$e$ term retained in the full stability calculation.

Small resonant-argument librations. Define the nominal resonant semi-major axis and the displaced center byThis follows by expanding $n=n_r[1-3(a-a_r)/(2a_r)]$ and taking $\cos (\pi +\delta )=-1+O({\delta }^2)$. It gives the requested initial semi-major axis to first order in the precession-induced detuning; $e_0=e(0)$ and $\delta =\Delta \varphi$.

For the pendulum approximation of a mean-motion resonance, also require fractional changes in $e$ small across a libration, and precession terms small compared with $n$. In particular $C_r/(n_r{e}_0^3)\ll 1$ for $q=1$ makes the extra full-equation curvature negligible; for $q=2$, small fractional eccentricity changes require $C_r/(n_r{e}_0^2)\ll 1$. Freeze $e=e_0$ in the leading restoring coefficient and setThen $u=\varphi -\pi$ obeys simple harmonic motion. Using $u(0)=\delta$, $\dot{u}(0)=0$ and integrating the remaining Lagrange planetary equations givesThe omitted apsidal modulation comes from higher perturbative orders and $u^2$. All expressions are leading resonant approximations, valid while the omitted fractional changes remain small. The small-libration phase portrait iswith period $2\pi /{\omega }_0$. The semi-major axis half-width is $2a_r{\omega }_0|\delta |/(3p{n}_{\mathrm{pl}})$ and the orbital eccentricity half-width is $q{C}_r{e}_0^{q-1}|\delta |/{\omega }_0$.

Dividing the original $\dot{e}$ and $\dot{a}$ equations gives the exact invariant of their retained terms,Thus the $e$–$a$ plot is a segment of this increasing curve traversed back and forth, not a closed ellipse. Locally its slope is $q/[2(p+q)e_0{a}_c]$.

Finite amplitude and the separatrix. Keep the same weak-resonance approximation but allow $\delta$ to be finite. The leading angular equation is $\ddot{\varphi }={\omega }_0^2\sin \varphi$, or $\ddot{u}=-{\omega }_0^2\sin u$. Multiplication by $\dot{u}$ givesThis resonant pendulum energy is constant in the pendulum approximation of a mean-motion resonance. It is not an exact integral of the earlier full equation when $e$ and its restoring coefficient vary: the $\dot{\varphi }\sin \varphi$ term alone prevents that conclusion in general.

For $|\delta |<\pi$ the phase portrait has closed libration curves around $\pi$, withwhere $K$ is the complete elliptic integral of the first kind. The largest speed is $2{\omega }_0\sin (|\delta |/2)$ and the leading semi-major axis half-width is $4a_r{\omega }_0\sin (|\delta |/2)/(3p{n}_{\mathrm{pl}})$. As $|\delta |\to \pi$, the energy approaches $2{\omega }_0^2$, the separatrix through the saddles at zero and $2\pi$; the particle spends increasingly long intervals near those saddles. With $ϵ=\pi -|\delta |$,At exactly $\delta =\pi$, $\dot{\varphi }=0$, the particle stays at the unstable fixed point in the ideal model. Nonstationary separatrix trajectories approach the saddle only in infinite time; above separatrix energy, the resonant argument circulates and resonance protection is lost. These distinctions matter when describing the limiting motion.

For unit mass, the inverse-square force has potential energy $-k/r$, with zero potential at infinity. Conservation of angular momentum gives $r^2\dot{\theta }=h$. Differentiate the polar equation of the Kepler orbit:Consequently the total mechanical energy iswhere $k=h^2/\ell$. This also explains the sign distinction between an elliptic orbit, a parabolic Kepler orbit, and a hyperbolic Kepler orbit.

At the original periapsis, the parabolic Kepler orbit has $e=1$, $r=d=\ell /2$, and purely tangential speed $V$. ThusThe outward radial impulse has zero moment about the Sun, so it preserves angular momentum and hence $\ell$. Its additional velocity is perpendicular to the original velocity, so the new mechanical energy is $E^{\prime }={\alpha }^2{V}^2/2$. Applying the energy formula givesFor $\alpha >0$, the new orbit is a hyperbolic Kepler orbit; $\alpha =0$ leaves the parabola unchanged.

To fix the orientation of the sketch, put the Sun at the origin, the impulse point at $(d,0)$, and take the original motion there upwards. The new velocity is $(\alpha V,V)$. The eccentricity vector isIt points towards the new periapsis, which is rotated clockwise through $\beta =\arctan (2\alpha )$ from the original periapsis direction. Thus, measured using the original polar angle $\theta$, the new orbit isIn particular, the impulse point is already on the outgoing part of the new orbit and is not its new periapsis. The old parabola satisfies $x=d-y^2/(4d)$. Both orbits pass through $(d,0)$ and have the same focus, but the hyperbola's new periapsis axis is tilted.

For a spherical grain of density $\rho$, the radiation-pressure coefficient isThe $r^{-2}$ dependence of radiation pressure cancels that of stellar gravity. Here radiation-pressure efficiency includes absorption and the momentum-transfer part of scattering, averaged over the stellar spectrum. In the geometrical optics regime, $D$ exceeds the important stellar wavelengths, $Q_{\mathrm{pr}}$ is of order unity, and $\beta \propto D^{-1}$: cross-sectional area grows as $D^2$, whereas mass grows as $D^3$. Around the stellar wavelength, Mie scattering can produce a broad maximum and material-dependent structure. In the Rayleigh scattering regime, absorption gives $Q_{\mathrm{abs}}\propto D$, while scattering gives $Q_{\mathrm{sca}}\propto D^4$. Thus small absorbing grains approach an approximately constant $\beta$, whereas nearly transparent grains have $\beta \propto D^3$. A universal fall to zero at small $D$ is therefore inappropriate without specifying the optical properties. The sketch shows both possible small-grain limits; its vertical normalization is illustrative.

Assume release with negligible velocity relative to the planet, negligible planetary gravity after release, and a constant radiation-pressure coefficient. Put $\mu =G{M}_{\star }$ and ${\mu }_d=(1-\beta )\mu$. The inherited speed and specific angular momentum are $v_p=\sqrt{\mu /a_p}$ and $h=a_p{v}_p$. The new specific orbital energy isFor $0\le \beta <1/2$, the dust orbit released from a circular parent ring is an ellipse. Using ${\epsilon }_d=-{\mu }_d/(2a_d)$ and $h^2={\mu }_d{a}_d(1-e_d^2)$ givesFor $\beta >0$, release is at periapsis, since $a_d(1-e_d)=a_p$. At the radiation-pressure blowout threshold, $\beta =1/2$, the orbit is parabolic. For $1/2<\beta <1$ it is hyperbolic: the signed semimajor axis is negative and $e_d>1$. For $\beta \ge 1$ the central force ceases to be attractive, so the elliptical interpretation and the positive-eccentricity formula cannot be extended unchanged.

With ${\theta }_d$ measured from release, it is the true anomaly. The polar equation of a Kepler orbit and conservation of specific angular momentum giveDefine the signed dust-tail angular lag by $\Delta \theta ={\theta }_d-{\theta }_p$, with ${\theta }_p=n_{p}t$ and both angles unwrapped. For fixed ${\theta }_d$, expansion to first order in $\beta$ givesso $\Delta \theta =2\beta (\sin {\theta }_d-{\theta }_d)+O({\beta }^2{\theta }_d)$. The dust trails the planet, hence the signed lag is negative. Its slope is $-4\beta {\sin }^2({\theta }_d/2)$; it decreases monotonically, has horizontal tangents at successive release-direction passages, and has a sinusoidal ripple about $-2\beta {\theta }_d$. At the end of orbit $j$, ${\theta }_d=2\pi j$ and $\Delta \theta \simeq -4\pi \beta j$. At five orbits this is $-0.2\pi$ for $\beta =0.01$ and $-0.8\pi$ for $\beta =0.04$. The upper-$\beta$ sketch remains a first-order approximation; its accumulated error is $O({\beta }^2{\theta }_d)$.

For the subsequent positive-distance statements, set $s=-\Delta \theta ={\theta }_p-{\theta }_d\ge 0$. The mean motion of the grain isA full relative wrap, $s=2\pi$, therefore takes $t_{\mathrm{wrap}}\simeq \pi /(\beta n_p)$ after averaging over the orbital ripple. During this time Poynting–Robertson drag changes the semimajor axis byThe signed change is negative. This estimate treats $a_d\simeq a_p$ as fixed in the rate and computes the drag accumulated over the radiation-pressure wrap time. Since $v_p/c\ll 1$, drag causes little migration over that time. For an arbitrarily tiny $\beta$, however, this small migration can itself alter the relative phase appreciably: neglecting that feedback on the wrap time additionally requires $v_p/c\ll \beta$.

Rapid dust sublimation is a plausible sink near a hot planet: a grain can evaporate long before it completes a relative wrap. Destructive collisions can also remove visible grains. Poynting–Robertson drag alone, with the estimate above, does not explain immediate removal from a short tail. If the destruction time is $t_{\mathrm{loss}}$, a short tail requires roughly $2\beta n_p{t}_{\mathrm{loss}}\ll 2\pi$.

For a steady dust-tail continuity equation, assume a constant injection rate $\dot{N}$, grains with the same fixed $\beta$, negligible release-speed dispersion, and negligible destruction within the segment being calculated. A finite destruction lifetime can terminate that segment, or multiply its density by a survival probability. To first order,To obtain the explicit trigonometric profile, make the additional secular phase approximation for a dust tail, $s\simeq 2\beta {\theta }_d$, discarding the bounded $-2\beta \sin {\theta }_d$ term in the angle-to-age map while retaining the periodic instantaneous drift speed. ThenSince the square is even, the same written profile applies to the signed $\Delta \theta$ if number intervals are interpreted with positive orientation. The result is an approximate phase-substitution profile, not a uniformly valid consequence of $\beta \ll 1$. The consistent first-order description is instead parametric:In particular, immediately after release, $s\sim \beta {\theta }_d^3/3$ and $n\propto s^{-2/3}$, whereas the phase-substitution expression would give $s^{-2}$. Nor does $s\ll 2\pi$ by itself justify neglecting the orbital ripple. The formal enhancements where $\dot{s}=0$ are dust-tail caustics; a spread of radiation-pressure coefficients, finite release velocities, and destruction smooth them. With constant lifetime $\tau$, the parametric profile acquires $\exp [-t({\theta }_d)/\tau ]$. These qualifications state precisely the extra assumptions behind the explicit profile.

For brevity write $b_1=b_{3/2}^{(1)}(\alpha )$, $b_2=b_{3/2}^{(2)}(\alpha )$, and $K=G{M}_p\alpha /(4a_p)$. The quadratic disturbing function is $\mathcal{R}=K[\frac{1}{2}{b}_1{e}^2-b_2{e}_{p}e\cos ({\varpi }_p-\varpi )]$. Differentiation in Lagrange planetary equations givesThe separate longitude of pericentre equation is singular at $e=0$, because a circular orbit has no distinguished periapsis. The complex eccentricity, $z=e{e}^{i\varpi }$, is regular there. Using $z_p=e_p{e}^{i{\varpi }_p}$ and $\dot{z}=e^{i\varpi }(\dot{e}+ie\dot{\varpi })$ givesSince $n{a}^2=\sqrt{G{M}_{\star }a}$ and $n_p=\sqrt{G{M}_{\star }/a_p^3}$,The approximations are planar, small-orbital eccentricity, first order in the planet-to-star mass ratio, and secular averaging away from important mean-motion resonances. The semimajor axis values remain constant at this order.

For a planet with fixed orbit, set $g=A{b}_1$ and $z_f=(b_2/b_1)z_p$. Then $\dot{z}=ig(z-z_f)$, so $z-z_f=(z(0)-z_f)e^{igt}$. With $z(0)=0$,On the Argand diagram, the complex eccentricity traces a circle centred on the forced eccentricity $z_f$, of radius $|z_f|$. For an interior test particle, $0<\alpha <1$ and $g>0$, so the motion is counterclockwise, starting at the origin. The rotating radius $z-z_f$ is the free eccentricity. The orbital eccentricity ranges from zero to $2|z_f|$, attaining its maximum at half a secular period. If $z_p$ is placed on the positive real axis, the initial tangent points downward, which is consistent with counterclockwise motion from the leftmost point of this circle.

Now let the planetary mass grow as $M_p(t)=M_{p}t/t_p$ for $0\le t\le t_p$, while the planet's orbital elements and the test particle's semimajor axis remain fixed. Both the precession and forcing coefficients acquire the same factor $t/t_p$. Their ratio, and hence $z_f$, remain constant. ThusandThis is the secular phase integral for proportional mass growth. The Argand diagram circle and the free eccentricity amplitude are unchanged; only the rate of traversal changes. The angular speed about $z_f$ starts at zero and rises linearly to $g$. The early displacement is quadratic in time rather than linear. Since $t_p=2\pi /g$, $\Phi (t_p)=\pi$, soA planet already at its final mass would instead return the particle to $z=0$ at $t_p$. After growth stops, the circle continues at rate $g$, with $z(t)=z_f[1+e^{ig(t-t_p)}]$ for $t\ge t_p$. Slow mass growth does not damp the free eccentricity here: the fixed centre of the circle and the conservative equation preserve its amplitude.

For several planets with fixed masses, Laplace-Lagrange secular theory gives $\dot{\mathbf{z}}=i\mathbf{Az}$ for the planets andfor a massless particle, with real secular coupling coefficients ${\nu }_j$. Resolve planetary complex eccentricities into secular eigenmodes, $z_j(t)={\sum }_k{E}_{jk}{e}^{i{g}_{k}t}$. Away from secular resonances, the test-particle solution has the formIts Argand diagram trajectory is a superposition of rotations: generally a quasiperiodic curve with beats, rather than a circle with one fixed centre. Even if the planetary orbits are artificially held fixed, their combined forcing shifts the centre to the sum of the individual forced eccentricities. At $g_k=g_0$, the resonant contribution is proportional to $t{e}^{i{g}_{0}t}$ in linear theory; its growing amplitude eventually invalidates the small-orbital eccentricity approximation.

If one planetary mass changes, the Laplace-Lagrange secular matrix, its eigenvalues, its eigenvectors, and the particle's forced response generally change with time. Planetary modes need not remain independent. Frequencies drift and a secular resonance can sweep across the particle, changing or exciting its free eccentricity. Very slow evolution away from resonances permits an approximate adiabatic mode description; rapid evolution or a resonance passage can mix modes. The simple substitution of a phase integral works only in the special proportional-growth case above, where the entire particle equation is multiplied by one scalar and its forced centre stays fixed. It does not apply to a general many-planet system.

The transverse velocity is $v_{\theta }=h/r$. Subtracting its square from the vis-viva equation, and putting $u=a/r$, givesConsequentlyThe positive branch is outward motion and the negative branch is inward motion. The PDF's unsigned expression gives the radial speed, or the outward branch, rather than the signed radial velocity everywhere. For example, on the inward part of any noncircular elliptic orbit, $\dot{r}<0$.

The same result follows from $\dot{r}=(\mu e/h)\sin f$, with $f$ the true anomaly. The radicand is nonnegative precisely on $1/(1+e)\le u\le 1/(1-e)$; it vanishes at periapsis and apoapsis. For a circular Kepler orbit, $e=0$, $u=1$ and $\dot{r}=0$.

Use the separation vector $\mathbf{r}$ from the star to the comet, and let $m_{\mathrm{red}}=m{M}_{\star }/(m+M_{\star })$. The two-body problem separates into uniform centre of mass motion andTaking a cross product with $\mathbf{r}$ proves conservation of angular momentum; taking the scalar product with $\dot{\mathbf{r}}$ proves conservation of energy. The relative constants areThe total relative angular momentum and energy are $m_{\mathrm{red}}\mathbf{h}$ and $m_{\mathrm{red}}\epsilon$, respectively. The conserved eccentricity vector is $\mathbf{e}=\dot{\mathbf{r}}\times \mathbf{h}/\mu -\hat{\mathbf{r}}$. Its scalar product with $\hat{\mathbf{r}}$ gives $r=h^2/[\mu (1+e\cos f)]$, with $f$ the true anomaly. The periapsis and apoapsis distances are therefore $h^2/[\mu (1+e)]$ and $h^2/[\mu (1-e)]$. Their sum is twice the semi-major axis, soAt either apsis $v=h/r$; substituting an apsidal radius gives $\epsilon =-\mu /(2a)$. Thus the specific orbital energy implies the vis-viva equationThese are exact relative-coordinate formulas for an elliptic orbit, $0\le e<1$. In the barycentric frame the comet's speed is $M_{\star }/(M_{\star }+m)$ times the relative speed; the printed identification with the comet's stellar orbital speed uses $M_{\star }\gg m$.

In the weak, non-crossing, leading-eccentricity model the exterior resonant coefficient has positive sign: write its disturbing function as $R_{\mathrm{res}}=B_q{e}_c^q\cos {\varphi }_{bc}$, with $B_q>0$ for the orders considered. The physical feedback can be seen for first order. Before an outer conjunction near apoapsis, the faster inner planet pulls the outer one backwards; after passing it pulls forwards. On the outward side, the nearer, stronger earlier pull wins, reducing the outer semi-major axis and increasing its mean motion. Thus ${\varphi }_{bc}<\pi$ is driven upward. On the inward side, the later forward pull wins and drives ${\varphi }_{bc}>\pi$ downward. The paired impulses cancel at the symmetric centre.

For higher order, sum the feedback over all $q$ conjunction phases. The first surviving angular Fourier harmonic is $\sin {\varphi }_{bc}$; the lower Fourier harmonics cancel over the equally spaced directions. A local Hamiltonian reduction with canonical momentum $J$ hasThe negative curvature of the Kepler orbit Hamiltonian gives the displayed sign. Linear stability analysis at $\pi$ gives a restoring acceleration. The symmetric libration centres of an eccentricity resonance are thereforeThe corresponding conjunction phases relative to periapsis are $\{\pi \}$, $\{\pi /2,3\pi /2\}$ and $\{\pi /3,\pi ,5\pi /3\}$. They suppress repeatedly close passages and provide resonance protection.

These are expected centres of the specified leading-harmonic model, not centres fixed solely by the integer order for arbitrary orbital eccentricity and masses. Asymmetric resonant-argument libration can occur when higher Fourier harmonics matter. For example, if $R=A\cos \varphi +B\cos 2\varphi$, $A>0$, $B>A/4$, the symmetric $\pi$ equilibrium loses stability and stable centres satisfy $\cos \varphi =-A/(4B)$. Thus further dynamical information is needed beyond the PDF's unrestricted eccentric orbit to assert unique centres. Exterior asymmetric branches are documented by Winter and Murray, “Resonance and chaos. II. Exterior resonances and asymmetric libration”.