The symmetry of the Kepler problem, the inverse Ligon-Schaaf mapping and the Birkhoff conjecture

The Ligon-Schaaf regularization (LS mapping) was introduced in 1976 and has been used several times. However, we are not aware of any direct usage of the inverse mapping, perhaps since it appears at first sight to be quite complex, involves the use of a transcendental equation (referred to as the generalized Kepler equation) that cannot be solved in closed form, and lacks smoothness near the collision point. Here, we provide some insight into the significance of this equation, along with a very simple derivation and confirmation of the inverse LS mapping. Then we use numerical methods to investigate three applications: 1) solutions of the Kepler function, 2) calculation of orbits including time-of-flight data based on the Delaunay Hamiltonian, and 3) numerical evidence for the Birkhoff conjecture for the circular restricted 3-body problem.

Proof. This is a simple substitution. ∎ Definition S2. For negative energy and ≠ 0, the inverse LS mapping is defined by (i) Mapping Proof.
(i) Derivation of Φ −1 = −1 ∘ −1 by substitution. This is the Cushman and Bates [1] notation for the derivation of = 1 ∘ 2 in the main paper. Note: This is exercise 12 in Cushman and Bates [1], which has one typographical error with <r,r> instead of <s,s>.

Identities
Proof of Proposition 3.
We'll begin with equation (30). Substitute from mapping F1, equations (7) and (10) Substitute the definition of H, equation (1) Collect similar terms In contrast with the original paper/thesis, we will confirm the following three identities by using mapping F2 and the previous three identities. This is much simpler than using the full mapping F.
Substitute from mapping F2, equations (13) and (14) Expand and collect similar terms Then we continue with equation (34) • Substitute from mapping F2, equations (13) through (16) Expand and collect similar terms Use previous identities and trigonometric identity Then we continue with equation (35) Substitute from mapping F2, equations (15) and (16) Expand and collect similar terms Use the previous identity of 2 = − 1 2 Since this identity is so important, we will also calculate it using the full forward mapping. We begin with the right-hand side of equation (36) = 0 ( ) − 0 ( ) Substitute from mapping F, equations (3) and (5) Expand and collect similar terms The next identity, equation (37), is important, because it appears as a denominator in equation (29) of the inverse LS mapping.
Substitute from mapping G1, equations (24) and (26) Use previous identities Resolve square root, choosing a minus sign.
Then, we continue with the identity for angular momentum (38).
Continue with the Hamiltonian of the rotating Kepler problem (41). This is just a combination of the previous identities for the Kepler Hamiltonian and the angular momentum and requires no further calculation. This completes the proof of Proposition 3. ∎
We begin with the definition of q from G1, substitute the values or r and s from G2, expand the expression, collect similar terms and use a trigonometric identity, use the previous identity of 2 = − 1 2 Now we take the definition of p from G1.
and substitute the values or r and s from G2.
Now we can solve (47) for q: Then we can solve (48) for p: With this, we can solve (50) for q: substitute q, • and p from above and (49), and simplify.
The result is now mapping G1.
(iii) and (iv) This is analogous to the case of negative energy. This completes the proof of Proposition 5. ∎
Then we substitute from mapping F1, equations (47) and (48) Substitute the definition of H, equation (1) Collect similar terms Then we continue with equation (69) Then we continue with equation (70) Substitute from mapping F1, equations (49) and (50) Substitute the definition of H, equation (1) Collect similar terms In contrast with the original paper/thesis, we will confirm the following three identities by using mapping F2 and the previous three identities. This is much simpler than using the full mapping F.
Substitute from mapping F2, equations (53) and (54) Expand and collect similar terms Then we continue with equation (72) Expand and collect similar terms Then we continue with equation (73) (55) and (56) Expand and collect similar terms Substitute from mapping F1, equation (49) = − 0 Now we begin with the right-hand side of equation (74) = 0 ℎ( ) + 0 √−η 2 ℎ( ) Substitute from mapping F2, equations (43) and (45) Expand and collect similar terms Since this identity is so important, we will also calculate it using the full forward mapping. We begin with the right-hand side of equation (74) = 0 ℎ( ) + 0 √−η 2 ℎ( ) Substitute from mapping F, equations (43) and (45) Expand and collect similar terms The next identity, equation (75), is important, because it appears as a denominator in equation (67) of the inverse LS mapping.
Substitute from mapping G1, equations (63) and (64) Use previous identities Resolve square root, choosing a minus sign.
Then, we continue with the identity for angular momentum (76).
Use previous identity for = √2 Substitute interim result for q from equation (47) = 1 Substitute from mapping G1, equations (63) and (64) Expand and collect similar terms Substitute from mapping G2, equations (59) through (62) Collect similar terms Use the previous identity of 2 = − 1 2 and trigonometric identity = ( 0 − 0 ) The next identity, equation (78), is just a restatement of the previous identity for 2 .
Continue with the Hamiltonian of the rotating Kepler problem (79). This is just a combination of the previous identities for the Kepler Hamiltonian and the angular momentum and requires no further calculation. This completes the proof of Proposition 6. ∎

Applications Investigate Kepler function
Our numerical solution is comprised of the two programs that were used to create S2

Calculate Orbits and Time of Flight
We provide software that accepts the coordinates p and q in phase space and calculates coordinates on the cotangent bundle of the sphere (after the LS mapping), the anomalies, Kepler elements, and Delaunay variables. It also calculates orbits, including time-of-flight data, by solving the Delaunay Hamiltonian on the cotangent bundle of the sphere and using the inverse LS mapping to convert the result to normal phase space. This last mapping implicitly solves the Kepler equation. The button show orbit was used to create the images in S1 Fig. The MATLAB App Designer program LSOrbit.mlapp accepts values for q and p and calculates basic values (including ξ and η), anomalies, Kepler elements, and Delaunay variables (see details below). The button show orbit creates a time vector of "num. points" and, for all points, converts them using the forward mapping, uses the Delaunay flow to calculate the orbit based on time t, and then converts back to phase space using the inverse mapping. Then it displays the orbit and the Delaunay variables for all points of the orbit. The button testDV creates "num. points" points based on q, q/2, q/4 etc., and calculates and displays the Delaunay variables for those points.
The calculations are based on the following formulas: LS forward mapping F: equation (1) Angular momentum G, Cordani [2] pg. 22.

Birkhoff conjecture for the circular restricted 3-body problem
Proof of Proposition 8. For the Kepler Hamiltonian, these expressions can be calculated in closed form. The Hamiltonian: The gradient of the Hamiltonian: The Hessian of the Hamiltonian: 2 ) The tangential component of the gradient: The curvature on the tangential hyperspace: We'll calculate that last expression. * ( ) * = 8 (‖ ‖ 2 + ‖ ‖ 2 ) 9 ( The curvature tangential to the hyperspace: In order to calculate the last expression, we'll begin by introducing some new notation, which consists of giving a name to each of the 2 terms of the gradient and the 3 terms of the Hessian. That results in a total of 12 terms for the product * ( ) * . Then we calculate each of these 12 terms singly, simplify and combine them, and finally calculate the determinant.
As a preliminary step, we observe where is the semimajor axis. Then we conclude By examining the ellipse, we can see that the semimajor axis < (| |) so, since | | < 1 for all points in the bounded component of the Hill region, < 1.

Then we recall
where is the eccentricity and 0 < < 1 for an ellipse, and we conclude altogether: Now, we make some observations and calculate some special cases for 5 . First, we state the main constraint We used numerical optimization, specifically MATLAB fmincon, to calculate the maximum value of 5 assuming that 0 ≤ ≤ 1, −1 ≤ ≤ 1, and ≤ − 3 2 . The result is (approximately) max ( 5 ) = 0 = 1 = −1.
The conditions show us the limits for the resulting optimum, so we will define them as constants. First, we simultaneously solve the following equations for X: 2 ( 6 + 12 3 − 1) −X ≤ L ≤ X L ≤ −  For Sun-Mars, the grid used did not find any points of negative curvature; however, our other calculations did find a few. For Sun-Earth, the grid used found a few points (8 to be exact) with negative curvature and relative distance to the Sun of less than .01.
For Sun-Jupiter, S5 Fig shows the distribution of energy and curvature of the energy, first over a wide range of points, then over a much smaller range of points centered around one point that has negative curvature. we can see that the energy (HC for circular restricted Hamiltonian, red points) is very close to the energy of the Lagrangian point L1 when the third body is close to L1 (q near 1). When the body is close to the heavy primary (q near 0), the potential energy can achieve large negative values, but can also approach the energy of L1, since it can be balanced by a large kinetic energy. It appears as though this large variance of the energy also makes a negative curvature possible.
In S5  Description of Software S1_Software.zip MATLAB software for calculations and plots. ".m" is a MATLAB script and can be opened in a standard text editor. Most of the code is fairly easy to read without experience in MATLAB,