Radial actions of the Kepler problem and the isotropic oscillator

A particle of mass m > 0 in a central potential U with angular momentum ℓ moves radially in the effective potential U_ℓ(r) = ℓ² / (2 m r²) + U(r). At energy E its radial momentum is p_r = ±√(2 m (E - U_ℓ(r))), and the radial motion oscillates between two turning points 0 < r₁ < r₂, the zeros of 2 m (E - U_ℓ). Between them one has the classical quadratures:

• the radial action I_r = (1 / π) ∫_{r₁}^{r₂} √(2 m (E - U_ℓ(r))) dr, the area enclosed by the radial phase curve divided by 2π (TauCeti.radialAction);
• the radial period T_r = 2 ∫_{r₁}^{r₂} m dr / √(2 m (E - U_ℓ(r))), twice the time from pericentre to apocentre, since dr/dt = p_r / m (TauCeti.radialPeriod);
• the apsidal angle Δθ = ∫_{r₁}^{r₂} ℓ dr / (r² √(2 m (E - U_ℓ(r)))), the polar angle swept from pericentre to apocentre (TauCeti.apsidalAngle).

All three are integrals of TauCeti.absRadialMomentum, the function r ↦ √(2 m (E - U(r)) - ℓ² / r²) = √(2 m (E - U_ℓ(r))).

This file evaluates them in closed form for the two potentials singled out by Bertrand's theorem, whose bounded orbits are all closed. Each turning-point hypothesis is the equation 2 m (E - U(rᵢ)) = ℓ² / rᵢ².

• Kepler, U(r) = -κ / r: I_r = κ √(m / (-2 E)) - |ℓ|, hence E = -m κ² / (2 (I_r + |ℓ|)²); Kepler's third law T_r = 2π √(m / κ) a^(3/2) for the semi-major axis a = (r₁ + r₂) / 2; and Δθ = π sign ℓ, so the radial and angular frequencies agree up to the sign of ℓ.
• Isotropic oscillator, U(r) = m ν² r² / 2: I_r = (E / ν - |ℓ|) / 2, T_r = π / ν and Δθ = (π / 2) sign ℓ, so the radial frequency is twice the magnitude of the angular one.

In both cases the turning-point equations force 2 m (E - U_ℓ(r)) = c (r - r₁) (r₂ - r) / r² (in the variable r² for the oscillator), so every quadrature is one of the integrals of TauCeti.Analysis.SpecialFunctions.Integrals.SqrtQuadratic.

Main declarations

• TauCeti.absRadialMomentum, TauCeti.radialAction, TauCeti.radialPeriod, TauCeti.apsidalAngle, with TauCeti.absRadialMomentum_nonneg and TauCeti.absRadialMomentum_eq_zero_iff (vanishing at the turning points).
• TauCeti.kepler_kappa_eq, TauCeti.kepler_sq_eq_neg_two_mul_mul_mul (Vieta's formulas for the Kepler turning points), TauCeti.kepler_energy_neg, TauCeti.kepler_kappa_pos and TauCeti.kepler_sqrt_mul_sqrt_eq_abs.
• TauCeti.kepler_radialAction, TauCeti.kepler_energy_eq_neg_div_sq_radialAction_add_abs, TauCeti.kepler_radialPeriod, TauCeti.kepler_apsidalAngle.
• TauCeti.isotropicOscillator_energy_eq, TauCeti.isotropicOscillator_sq_eq_sq_mul (Vieta's formulas for the oscillator turning points) and TauCeti.isotropicOscillator_abs_eq_mul.
• TauCeti.isotropicOscillator_radialAction, TauCeti.isotropicOscillator_radialPeriod, TauCeti.isotropicOscillator_apsidalAngle.

References

• V. I. Arnold, Mathematical Methods of Classical Mechanics, 2nd ed., Graduate Texts in Mathematics 60, Springer, 1989, Sections 8 and 50.
• H. Goldstein, C. Poole, J. Safko, Classical Mechanics, 3rd ed., Addison-Wesley, 2002, Chapter 3 and Section 10.8.

noncomputable def TauCeti.absRadialMomentum (m : ℝ) (U : ℝ → ℝ) (ℓ E r : ℝ) : ℝ

The absolute radial momentum √(2 m (E - U(r)) - ℓ² / r²) at radius r of a particle of mass m with angular momentum ℓ and energy E in the central potential U. Where the radicand is nonnegative this is |p_r|, written as a function of r at fixed ℓ and E rather than as a phase-space coordinate; it vanishes at the turning points.

noncomputable def TauCeti.radialAction (m : ℝ) (U : ℝ → ℝ) (ℓ E r₁ r₂ : ℝ) : ℝ

The radial action (1 / π) ∫_{r₁}^{r₂} |p_r| dr between the turning points r₁ and r₂, the area enclosed by the radial phase curve divided by 2π.

noncomputable def TauCeti.radialPeriod (m : ℝ) (U : ℝ → ℝ) (ℓ E r₁ r₂ : ℝ) : ℝ

The radial period 2 ∫_{r₁}^{r₂} m dr / |p_r|: the radial velocity is p_r / m, so this is the time to travel from r₁ to r₂ and back.

noncomputable def TauCeti.apsidalAngle (m : ℝ) (U : ℝ → ℝ) (ℓ E r₁ r₂ : ℝ) : ℝ

The apsidal angle ∫_{r₁}^{r₂} ℓ dr / (r² |p_r|), the polar angle swept while the radius moves from r₁ to r₂: along the orbit dθ / dr = ℓ / (r² p_r).

theorem TauCeti.absRadialMomentum_def (m : ℝ) (U : ℝ → ℝ) (ℓ E r : ℝ) : absRadialMomentum m U ℓ E r = √(2 * m * (E - U r) - ℓ ^ 2 / r ^ 2)

Unfolds TauCeti.absRadialMomentum.

theorem TauCeti.radialAction_def (m : ℝ) (U : ℝ → ℝ) (ℓ E r₁ r₂ : ℝ) : radialAction m U ℓ E r₁ r₂ = Real.pi⁻¹ * ∫ (r : ℝ) in r₁..r₂, absRadialMomentum m U ℓ E r

Unfolds TauCeti.radialAction.

theorem TauCeti.radialPeriod_def (m : ℝ) (U : ℝ → ℝ) (ℓ E r₁ r₂ : ℝ) : radialPeriod m U ℓ E r₁ r₂ = 2 * ∫ (r : ℝ) in r₁..r₂, m / absRadialMomentum m U ℓ E r

Unfolds TauCeti.radialPeriod.

theorem TauCeti.apsidalAngle_def (m : ℝ) (U : ℝ → ℝ) (ℓ E r₁ r₂ : ℝ) : apsidalAngle m U ℓ E r₁ r₂ = ∫ (r : ℝ) in r₁..r₂, ℓ / (r ^ 2 * absRadialMomentum m U ℓ E r)

Unfolds TauCeti.apsidalAngle.

theorem TauCeti.absRadialMomentum_nonneg (m : ℝ) (U : ℝ → ℝ) (ℓ E r : ℝ) : 0 ≤ absRadialMomentum m U ℓ E r

The absolute radial momentum is nonnegative.

theorem TauCeti.absRadialMomentum_eq_zero_iff {m : ℝ} {U : ℝ → ℝ} {ℓ E r : ℝ} : absRadialMomentum m U ℓ E r = 0 ↔ 2 * m * (E - U r) ≤ ℓ ^ 2 / r ^ 2

The absolute radial momentum vanishes exactly where 2 m (E - U(r)) ≤ ℓ² / r²; in particular it vanishes at a turning point, where 2 m (E - U(r)) = ℓ² / r².

theorem TauCeti.absRadialMomentum_pos_iff {m : ℝ} {U : ℝ → ℝ} {ℓ E r : ℝ} : 0 < absRadialMomentum m U ℓ E r ↔ ℓ ^ 2 / r ^ 2 < 2 * m * (E - U r)

The absolute radial momentum is positive exactly where ℓ² / r² < 2 m (E - U(r)).

The Kepler problem

theorem TauCeti.kepler_kappa_eq {m κ ℓ E r₁ r₂ : ℝ} (hm : m ≠ 0) (h₁ : r₁ ≠ 0) (h₂ : r₂ ≠ 0) (hne : r₁ ≠ r₂) (hr₁ : 2 * m * (E + κ / r₁) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E + κ / r₂) = ℓ ^ 2 / r₂ ^ 2) : κ = -E * (r₁ + r₂)

Vieta's formula for the Kepler sum of turning points. If r₁ ≠ r₂ are nonzero turning points of a Kepler orbit with m ≠ 0, then κ = -E (r₁ + r₂); equivalently, for E ≠ 0 the semi-major axis (r₁ + r₂) / 2 is κ / (-2 E).

theorem TauCeti.kepler_sq_eq_neg_two_mul_mul_mul {m κ ℓ E r₁ r₂ : ℝ} (hm : m ≠ 0) (h₁ : r₁ ≠ 0) (h₂ : r₂ ≠ 0) (hne : r₁ ≠ r₂) (hr₁ : 2 * m * (E + κ / r₁) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E + κ / r₂) = ℓ ^ 2 / r₂ ^ 2) : ℓ ^ 2 = -2 * m * E * (r₁ * r₂)

Vieta's formula for the Kepler product of turning points. If r₁ ≠ r₂ are nonzero turning points of a Kepler orbit with m ≠ 0, then ℓ² = -2 m E r₁ r₂.

theorem TauCeti.kepler_energy_neg {m κ ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hℓ : ℓ ≠ 0) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E + κ / r₁) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E + κ / r₂) = ℓ ^ 2 / r₂ ^ 2) : E < 0

Two distinct positive turning points of a Kepler orbit with m > 0 and ℓ ≠ 0 force negative energy.

theorem TauCeti.kepler_kappa_pos {m κ ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hℓ : ℓ ≠ 0) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E + κ / r₁) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E + κ / r₂) = ℓ ^ 2 / r₂ ^ 2) : 0 < κ

Two distinct positive turning points of a Kepler orbit with m > 0 and ℓ ≠ 0 force an attractive potential: κ > 0.

theorem TauCeti.kepler_sqrt_mul_sqrt_eq_abs {m κ ℓ E r₁ r₂ : ℝ} (hm : m ≠ 0) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E + κ / r₁) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E + κ / r₂) = ℓ ^ 2 / r₂ ^ 2) : √(-2 * m * E) * √(r₁ * r₂) = |ℓ|

If 0 < r₁ < r₂ are turning points of a Kepler orbit with m ≠ 0, then √(-2 m E) √(r₁ r₂) = |ℓ|.

theorem TauCeti.kepler_radialAction {m κ ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hℓ : ℓ ≠ 0) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E + κ / r₁) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E + κ / r₂) = ℓ ^ 2 / r₂ ^ 2) : radialAction m (fun (r : ℝ) => -κ / r) ℓ E r₁ r₂ = κ * √(m / (-2 * E)) - |ℓ|

The Kepler radial action. Between two distinct turning points 0 < r₁ < r₂ of a Kepler orbit of mass m > 0 and angular momentum ℓ ≠ 0 in the potential -κ / r, the radial action (1 / π) ∫_{r₁}^{r₂} √(2 m (E + κ / r) - ℓ² / r²) dr is κ √(m / (-2 E)) - |ℓ|.

theorem TauCeti.kepler_energy_eq_neg_div_sq_radialAction_add_abs {m κ ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hℓ : ℓ ≠ 0) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E + κ / r₁) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E + κ / r₂) = ℓ ^ 2 / r₂ ^ 2) : E = -(m * κ ^ 2) / (2 * (radialAction m (fun (r : ℝ) => -κ / r) ℓ E r₁ r₂ + |ℓ|) ^ 2)

The Kepler energy in terms of the actions. Between two distinct turning points 0 < r₁ < r₂ of a Kepler orbit with ℓ ≠ 0, E = -m κ² / (2 (I_r + |ℓ|)²), where I_r is the radial action and |ℓ| the angular action.

theorem TauCeti.kepler_radialPeriod {m κ ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hℓ : ℓ ≠ 0) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E + κ / r₁) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E + κ / r₂) = ℓ ^ 2 / r₂ ^ 2) : radialPeriod m (fun (r : ℝ) => -κ / r) ℓ E r₁ r₂ = 2 * Real.pi * √(m / κ) * ((r₁ + r₂) / 2) ^ (3 / 2)

Kepler's third law. Between two distinct turning points 0 < r₁ < r₂ of a Kepler orbit with ℓ ≠ 0, the radial period is 2π √(m / κ) a^(3/2) for the semi-major axis a = (r₁ + r₂) / 2.

theorem TauCeti.kepler_apsidalAngle {m κ ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hℓ : ℓ ≠ 0) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E + κ / r₁) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E + κ / r₂) = ℓ ^ 2 / r₂ ^ 2) : apsidalAngle m (fun (r : ℝ) => -κ / r) ℓ E r₁ r₂ = Real.pi * ↑(SignType.sign ℓ)

The Kepler apsidal angle. Between two distinct turning points 0 < r₁ < r₂ of a Kepler orbit with ℓ ≠ 0, the polar angle swept from pericentre to apocentre is π sign ℓ: the orbit closes after one radial period, with no precession of its apsides.

The isotropic oscillator

theorem TauCeti.isotropicOscillator_energy_eq {m ν ℓ E r₁ r₂ : ℝ} (hm : m ≠ 0) (h₁ : r₁ ≠ 0) (h₂ : r₂ ≠ 0) (hne : r₁ ^ 2 ≠ r₂ ^ 2) (hr₁ : 2 * m * (E - m * ν ^ 2 * r₁ ^ 2 / 2) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E - m * ν ^ 2 * r₂ ^ 2 / 2) = ℓ ^ 2 / r₂ ^ 2) : 2 * E = m * ν ^ 2 * (r₁ ^ 2 + r₂ ^ 2)

Vieta's formula for the isotropic oscillator, sum of squared turning points. If r₁, r₂ are nonzero turning points of an isotropic oscillator with m ≠ 0 and r₁² ≠ r₂², then 2 E = m ν² (r₁² + r₂²).

theorem TauCeti.isotropicOscillator_sq_eq_sq_mul {m ν ℓ E r₁ r₂ : ℝ} (hm : m ≠ 0) (h₁ : r₁ ≠ 0) (h₂ : r₂ ≠ 0) (hne : r₁ ^ 2 ≠ r₂ ^ 2) (hr₁ : 2 * m * (E - m * ν ^ 2 * r₁ ^ 2 / 2) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E - m * ν ^ 2 * r₂ ^ 2 / 2) = ℓ ^ 2 / r₂ ^ 2) : ℓ ^ 2 = (m * ν * (r₁ * r₂)) ^ 2

Vieta's formula for the isotropic oscillator, product of squared turning points. If r₁, r₂ are nonzero turning points of an isotropic oscillator with m ≠ 0 and r₁² ≠ r₂², then ℓ² = (m ν r₁ r₂)².

theorem TauCeti.isotropicOscillator_abs_eq_mul {m ν ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hν : 0 ≤ ν) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E - m * ν ^ 2 * r₁ ^ 2 / 2) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E - m * ν ^ 2 * r₂ ^ 2 / 2) = ℓ ^ 2 / r₂ ^ 2) : |ℓ| = m * ν * (r₁ * r₂)

If 0 < r₁ < r₂ are turning points of an isotropic oscillator with m > 0 and ν ≥ 0, then |ℓ| = m ν r₁ r₂.

theorem TauCeti.isotropicOscillator_radialAction {m ν ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hν : 0 < ν) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E - m * ν ^ 2 * r₁ ^ 2 / 2) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E - m * ν ^ 2 * r₂ ^ 2 / 2) = ℓ ^ 2 / r₂ ^ 2) : radialAction m (fun (r : ℝ) => m * ν ^ 2 * r ^ 2 / 2) ℓ E r₁ r₂ = (E / ν - |ℓ|) / 2

The radial action of the isotropic oscillator. Between two distinct turning points 0 < r₁ < r₂ of a particle of mass m > 0 in the potential m ν² r² / 2 with ν > 0, the radial action (1 / π) ∫_{r₁}^{r₂} √(2 m (E - m ν² r² / 2) - ℓ² / r²) dr is (E / ν - |ℓ|) / 2.

theorem TauCeti.isotropicOscillator_radialPeriod {m ν ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hν : 0 < ν) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E - m * ν ^ 2 * r₁ ^ 2 / 2) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E - m * ν ^ 2 * r₂ ^ 2 / 2) = ℓ ^ 2 / r₂ ^ 2) : radialPeriod m (fun (r : ℝ) => m * ν ^ 2 * r ^ 2 / 2) ℓ E r₁ r₂ = Real.pi / ν

The radial period of the isotropic oscillator. Between two distinct turning points 0 < r₁ < r₂ in the potential m ν² r² / 2 with ν > 0, the radial period is π / ν, half the period 2π / ν of the oscillator.

theorem TauCeti.isotropicOscillator_apsidalAngle {m ν ℓ E r₁ r₂ : ℝ} (hm : 0 < m) (hν : 0 < ν) (h₁ : 0 < r₁) (hlt : r₁ < r₂) (hr₁ : 2 * m * (E - m * ν ^ 2 * r₁ ^ 2 / 2) = ℓ ^ 2 / r₁ ^ 2) (hr₂ : 2 * m * (E - m * ν ^ 2 * r₂ ^ 2 / 2) = ℓ ^ 2 / r₂ ^ 2) : apsidalAngle m (fun (r : ℝ) => m * ν ^ 2 * r ^ 2 / 2) ℓ E r₁ r₂ = Real.pi / 2 * ↑(SignType.sign ℓ)

The apsidal angle of the isotropic oscillator. Between two distinct turning points 0 < r₁ < r₂ in the potential m ν² r² / 2 with ν > 0, the polar angle swept from pericentre to apocentre is (π / 2) sign ℓ: the orbit closes after two radial periods.