Stability of the Rhomboidal Symmetric-Mass Orbit

We study the rhomboidal symmetric-mass 4-body problem in both a two-degree-of-freedom and a four-degree-of-freedom setting. Under suitable changes of variables in both settings, isolated binary collisions at the origin are regularizable. Linear stability analysis is performed in both settings. For the two-degree-of-freedom setting, linear stability is established for a wide interval of mass ratios. A Poincar\'{e} section analysis is also performed, showing stability. In the four-degree-of-freedom setting, linear stability fails except for a very small interval of mass ratios.


Introduction
In the Principia Mathematica, published in 1687, Newton outlined many governing principles of the motion of physical objects. Combining the laws F = ma and the law of universal gravitation gave a relation that could describe the motions of bodies in space. The resulting equations helped to explain many of the behaviors that astronomers of the time were aware of (most notably Kepler). Mathematically, the study of determining the motion of n point masses in space is known as the Newtonian n-body problem. Notationally, if {q 1 , q 2 , ..., q n } represent the positions of the bodies in R k (k = 1, 2, or 3) with masses {m 1 , m 2 , ..., m n } respectively, then their motion is governed by the system of differential equations where the dot represents the derivative with respect to time. Despite hundreds of years of study and the relatively recent development of computer ODE solvers, many open questions about the n-body problem remain.
One aspect of the n-body problem that has been getting much attention of late are orbits involving collision singularities. A collision singularity occurs when q i = q j for some i = j. In the equations governing motion, this results in a zero denominator in one or more terms in the sum. Under certain conditions, these collisions can be regularized and continued past collision.
Schubart [15] was one of the first to study periodic orbits with collisions. He was able to find a collinear three-body equal-mass orbit where the central body alternated between collisions with the outer two. This was further extended to the case of arbitrary masses numerically by Hénon [6] in 1977. Analytic existence of the equal-outer-mass orbit was established independently by Venturelli [21] and Mockel [11], both in 2008. Shibiyama [16] recently demonstrated the existence of the arbitrary-mass version. The study of linear stability of Schubart's orbit was performed by Hietarinta and Mikkola [7] in 1993.
Sweatman found a Schubart-like collinear four-body symmetric orbit in 2002 [19], and later studied its linear stability [20]. This orbit features simultaneous binary collisions between two outer pairs of bodies followed by an interior collision between the two central bodies. Analytic existence of this orbit was given by Ouyang and Yan in [12].
Planar orbits with singularities have also been studied. A planar four-body orbit featuring simultaneous binary collisions was described in [13]. The orbit was shown to be linearly stable in [4]. It was later shown that this orbit could be extended to symmetric masses in [2] (see also [3]), and linear stability for this extension was shown for an interval of certain mass ratios in [1].
Analytic existence of large families of orbits with singularities was recently proven by Shibayama in [16] and Martinez in [9]. Each orbit can be reduced to two position and two momentum variables (the so-called two-degree-of-freedom problem). One orbit of note in this family is the rhomboidal four-body orbit, which features two pairs of bodies on the xand y-axes. The pairs collide at the origin in an alternating fashion. This orbit was shown to exist analytically in multiple independent papers (by Yan in [23] and Martinez in [9] for equal masses, [16] for symmetric masses). Additionally, Yan showed that for equal masses, the orbit is linearly stable.
In a separate study of the rhomboidal four-body problem with unequal masses, Waldvogel [22] notes that "sufficiently simple systems may bear the chance of permitting theoretical advances," and identifies the rhomboidal configuration as one such system. Indeed, much of the analysis performed in this paper is far simpler than that of [1].
The remainder of this paper is divided into two principal sections. Section 2 concerns the orbit in the two-degree-of-freedom (2DF) setting. In Section 2.1 we give the notation and mathematical description of the orbit. Section 2.2 outlines some basic theory for linear stability. Section 2.3 describes some preliminary numerical calculations that are needed to study the orbit. Section 2.4 gives a Poincaré section analysis of the orbit in the 2DF setting. In Section 3 we further the study of the orbit in the four-degree-of-freedom (4DF) setting. Section 3.1 sets up the additional mathematical notation for the new setting. Section 3.2 describes the symmetries of the orbit, which will be needed for the stability calculations. Section 3.3 reviews how symmetries of a periodic orbit can be used to simplify the linear stability calculation. Section 3.4 shows how the stability calculation can be reduced to the calculation of three entries of a particular matrix. Section 3.5 gives the results and implications of the remaining calculation.
2 The Rhomboidal Two-Degree-of-Freedom Symmetric-Mass Problem

The Periodic Orbit
We consider the planar Newtonian 4-body problem with bodies located at and masses 1, m, 1, m respectively for some m ∈ (0, 1]. The bodies travel along the x and y axes, forming the vertices of a rhombus at all times away from collision. Binary collisions occur between the bodies with equal masses at the origin. For the periodic orbit, the non-colliding bodies have zero momentum at collision time. Following collision, the colliding bodies eject along the appropriate coordinate axis, and the remaining two bodies travel toward collision at the origin, where a similar (zero momentum of non-colliding bodies) behavior occurs. (See Figure 1.) Analytic existence of this orbit was shown in [23] in the m = 1 case, and for m ∈ (0, 1] in [16]. The Hamiltonian for the system is given by where w 1 = 2ẋ 1 , w 2 = 2mẋ 2 . We can continue the orbits past collision via a regularization under which binary collision corresponds to an elastic bounce. To regularize these collisions, we use a Levi-Civita-type change of coordinates. Using the canonical transformations Q 2 i = x i , P i = 2Q i w i for i = 1, 2, with a change of time satisfying dt/ds = x 1 x 2 , the regularized Hamiltonian in extended phase space is given by This yields the equations of motion where denotes the derivative with respect to the new time variable s.
At the time of collision of the two bodies on the x-axis, we have Q 1 = 0 and P 2 = 0. At this time, setting Γ = 0 yields Hence, P 1 = ±8 3/2 , with the sign being the same as the sign on Q 2 . Similarly, at the time of collision of the two bodies on the y-axis, Q 2 = P 1 = 0. The condition Γ = 0 then gives and so P 2 = ±(8m) 3/2 .

Let
A standard proof shows that if γ(s) is a T -periodic solution to (3) -(6), both −Sγ(T /2 − s) and Sγ(T − s) are solutions as well. Existence and uniqueness of solutions then imply that for all s. Hence the symmetry group for the rhomboidal four-body orbit is isomorphic to the Klein four group, with S and −S as generators.

Linear Stability
Note that Γ is a smooth function defined on R 4 \{Q 1 = Q 2 = 0}. Suppose γ(s) is a T -periodic solution of the system z = JDΓ(z), where = d/ds, and I and O are the 2 × 2 identity and zero matrices, respectively. If X(s) is the fundamental matrix solution of the linearized equations then the monodromy matrix is given by X(T ) and satisfies X(s + T ) = X(s)X(T ) for all s. Eigenvalues of the monodromy matrix are also the characteristic multipliers of γ, and therefore determine the linear stability of γ. In particular, γ is spectrally stable if all of its characteristic multipliers lie on the unit circle, and γ is linearly stable if it is spectrally stable and semisimple apart from trivial eigenvalues.
Linear stability is typically established by numerical integration. Some elegant techniques for simplifying the numerical work were presented by Roberts in [14], and will be presented in Section 3.3.

Numerical Determination of Initial Conditions
In order to determine the initial conditions for the rhomboidal orbit, we model each of Q 1 , Q 2 , P 1 , and P 2 by truncated trigonometric polynomials: The choice of trigonometric polynomials is natural for modeling periodic behavior. A similar technique was carried out by Simó in [18]. The time shifts and choice of odd-only multiples of s correspond to symmetries of the orbit. In particular, for these polynomials, the timereversing symmetries shown earlier are built-in, and the non-colliding bodies have zero net momentum at collision time. For a fixed n, we numerically minimize the value of where the minimization is taken over the space of coefficients We combine this with a root-finding technique to find the appropriate value of E for a 2π-periodic orbit, as in [2] (see also [3]). Once these trigonometric polynomials are determined, we can extract the initial conditions for the orbit by evaluatingQ 1 ,Q 2 ,P 1 ,P 2 at any fixed time s ∈ [0, 2π].
We obtain the initial conditions for the 2π-periodic orbit for m = 1 by rescaling the conditions given in [23]. It is easy to check that if γ(s) = (Q 1 (s), Q 2 (s), P 1 (s), P 2 (s), E) is a solution for (3) -(6), then the solution γ ε (s) = (εQ 1 (εs), εQ 2 (εs), P 1 (εs), P 2 (εs), E/ε 2 ) also satisfies (3) - (6). Moreover, this rescaling does not change the linear stability of the orbit. Given the initial conditions in [23] and rescaling with ε ≈ 1.55, the orbit is roughly 2π-periodic, which we verify by integration using the standard fourth-order Runge-Kutta-Fehlberg algorithm. A standard curve-fitting technique can then be used to give the coeffi- (12). After that, a gradual "step-down" technique (as in [2]) can be used to find the initial conditions for other values of m ∈ (0, 1]. If we assume the initial conditions occur at the time the two bodies with mass 1 collide, we have Q 1 = P 2 = 0. Also, as before, we know that P 1 = −2 √ 2, so we need only to find the values of Q 2 and E. The results of the numerical calculation are shown in Figures 2 -4. As a result of having determined the initial conditions for the orbit, we can perform a numerical integration to determine the linear stability. An elegant decomposition will yield the following stability result for the 2DF setting as a corollary to the stability in the 4DF setting. As such, we postpone the proof of the following until Section 3.5.
Theorem 1 There exists some positive number ε such that the 2DF symmetric-mass periodic orbit of the regularized planar rhomboidal four-body problem is linearly stable for m ∈ (.01 + ε, 1].

Poincaré Section Analysis
To numerically analyze nonlinear stability, we find a suitable Poincaré section for the orbit. This was done in the m = 1 case in [22]. Our more general Poincaré section is based on techniques presented in [7] and [19]. For any value of m, we seek a number α such that  is maintained throughout the orbit, with x 1 and x 2 as defined in (14) earlier. In other words, the value of α corresponds to the ratio of x 1 and x 2 in a homographic orbit where the trajectories of the bodies correspond to total collapse (or ejection from total collapse). We find this value of α by solving the standard equations of motion (1) with the substitutions x 1 = αx 2 andẍ 1 = αẍ 2 . Doing so, we find that the required value of α for a given mass m is a root of the 12th-degree polynomial Notice that if the ratio x 1 /x 2 is constant throughout the orbit, then the ratio x 2 /x 1 is also constant throughout. It can be verified by numerical integration that the roots of 13 corresponding to the ratio x 1 /x 2 lie in the interval [0, 1]. This will be preferred for ease of numerical calculation. The value of α as a function of m is plotted in Figure 5. For fixed E = −1, we define a Poincaré section Σ to be the two-dimensional surface given by x 1 = αx 2 in the phase space defined by the variables x 1 , x 2 ,ẋ 1 , andẋ 2 . (Note thatẋ 1 andẋ 2 are simply linear re-scalings of w 1 and w 2 .) Restricting to E = −1, we find a bound on the possible values of x 1 . Specifically, ifẋ 1 =ẋ 2 = 0 on Σ, the condition E = −1 requires that For a set of initial conditions on Σ, the requirement E = −1 necessarily implies that x 1 ≤ r max , and if either ofẋ 1 orẋ 2 are non-zero, then the strict inequality x 1 < r max holds.
We define coordinates (r, θ) on Σ by Under this change of coordinates, the homographic orbit corresponds to the line θ = π/4. For a 9 × 15 grid of equally spaced initial conditions in (r, θ) we numerically integrate (3) - (6) for the corresponding initial conditions and record the first 200 intersections of the orbit with Σ. (Integration was preemptively terminated if any of Q i , P i exceeded 1000 in absolute value.) The results of this are shown in Figures 6 -10. The observed concentric rings numerically match the predicted result of Moser's Invariant Curve Theorem in [17], and show that the rhomboidal symmetric-mass orbit is nonlinearly stable for m ∈ (.01 + ε, 1] for the same ε as in Theorem 1.

Description and Existence
We now consider the planar Newtonian 4-body problem with bodies located at and masses 1, m, 1, m respectively for some m ∈ (0, 1]. (It is important to note here that x 2 does not correspond to x 2 from the previous section.) For the periodic orbit, the bodies still travel along the x and y axes, forming the vertices of a rhombus at all times away from collision, with the same behaviors of the 2DF orbit (such as zero momentum of non-colliding bodies at collision time) still holding.
The Hamiltonian for this system is given by where the w i are the conjugate momenta defined by The angular momentum for the system is given by We can regularize the system under a change of spatial variables and a re-scaling of time. Define Then F induces the canonical change of variables (x i , w i ) ↔ (Q i , P i ) given by Each of the P i is linear in w i . Solving the resulting system of equations yields We can regularize the collisions at the origin by multiplying by a change of time satisfying . At the time of collision between the two bodies of mass 1, we have Q 1 = Q 2 = 0. The condition and so at collision the momenta P 1 and P 2 are both finite and satisfy P 2 1 + P 2 2 = 8. Similarly, when Q 3 = Q 4 = 0, we get so the momenta P 3 and P 4 are both finite and satisfy P 2 3 + P 2 4 = 8m 3 .
Let A denote the set where This corresponds to the regularized coordinates Then, when A holds, the four bodies and their respective momenta lie on the xand y-axes, as in the two-degree of freedom (2DF) problem. We also have which are the same coordinate transformations used in [23] and in our work in Section 2. Furthermore, we haveQ so A is invariant, and corresponds to the 2DF rhomboidal configuration. Hence, the 2DF problem embeds nicely into the 4DF problem, and initial conditions from the 2DF problem can also be used to study the 4DF orbit. This result, combined with the existence of the 2DF orbit from [16], and [9], demonstrates the analytic existence of the 4DF rhomboidal orbit.

Symmetries of the Rhomboidal Four-Degree-of-Freedom Orbit
Let and define the block matrix where 0 represents the 2 × 2 identity matrix. Then we have

Stability Calculations with Symmetry
Continuing from the brief introduction given in 2.2, if Y is the fundamental matrix solution of ξ = JD 2 Γ(γ(s))ξ, ξ(0) = Y 0 for some invertible matrix Y 0 , then by definition of X(s), Y (s) = X(s)Y 0 , implying X(T ) = Y (T )Y −1 0 . Then we have and so X(T ) and Y −1 0 Y (T ) are similar, and stability can be determined by the eigenvalues of either. For our purposes, the latter will be preferred.
The following can be found in [14]:
Then the fundamental matrix solution X(s) satisfies X(−s + T /N) = SX(s)S T (X(T /N)).
Note that the matrix S given in (16) satisfies all the required hypotheses.  Similar results for time-preserving symmetries are also presented in [14], but are not needed for this orbit. Using these results may allow the computation of the eigenvalues (hence stability) to be accomplished using only a fraction of the orbit. Applying Corollary 2 with N = 2, S as defined in (16), and noting that S T = S yields Similarly, if N = 1, since S 2 = I, we get

This yields
Hence, in order to analyze the stability of the orbit, we need only compute the eigenvalues of Y along a quarter of the orbit.
Again, from [14]: Lemma 2 For a symplectic matrix W , suppose there is a matrix K such that Then W is stable if and only if all of the eigenvalues of K are real and have absolute value less than or equal to 1.
We now show that there is an appropriate choice of Y 0 for which W has the required form, further reducing the stability calculations for the orbit. If we let Also, since Λ 2 = D 2 = I, we know immediately that Thus, Similarly, we find that Thus, we have for some 4 × 4 matrix K.
Remark: The given matrix Y 0 in (18) is not unique. Different choices of Y 0 will give different properties of the monodromy matrix. It is also worth noting that Y 0 is independent of the value of m for this orbit, which is not always true (see [1].) We can give formulas for the entries of K in terms of W . Since B is symplectic, we have J = B T JB, and hence B −1 = −JB T J.
Using W = ΛD for D = −B −1 SB and the relation −SJ = JS, we find Directly computing ΛJ and using the block form of B, we find that Define col i (−SJB) to be the ith column of the matrix −SJB. Then we have col i (−SJB) = −SGc i where c i is the ith column of B. Using the above two formulas, this implies that the (i, j) entry of W is given by −c T i SJC j . Equation (19) shows that the (i, j) entry of K is the (i + 4, j + 4) entry of W . Hence, Remark: Computing the entries of K this way will allow us to bypass computing W −1 . This is preferred as a numerical method as W may be very poorly conditioned.
From the form of W , this implies that Here, the zeros denote entries for which the mixed partials evaluate to zero identically, and the entries denoted a are entries for which the mixed partials evaluate to zero assuming the conditions given by (15)  Proof We verify that Mη has the proper form. Note that Hence Remarks: (i) In terms of the 4DF Rhomboidal orbit, this form very nicely decomposes phase space into a direct sum of the subspaces A = {Q 2 = Q 3 = P 2 = P 3 = 0} and A ⊥ = {Q 1 = Q 4 = P 1 = P 4 = 0}. This decomposition is due in part to the coordinate transformation we chose. The choice of notation for A ⊥ is appropriate in that A ⊥ and A are orthogonal complements in R 8 . The two subspaces are also skew-orthogonal: if a 1 ∈ A and a 2 ∈ A ⊥ , then a T 1 Ja 2 = 0.
(ii) Matrices of the form M and M 2 are similar to the diamond product discussed in [8]. Specifically, Σ −1 MΣ = A 1 3A 2 for some matrices A 1 and A 2 , where M ∈ M 2 and Σ is the permutation matrix corresponding to σ = (1 2 3 4 5 6 7 8). Furthermore, one of these two matrices corresponds to the 2DF setting.
(iii) The particular choice of Y 0 given in (18) is again important for this argument.
Assuming the initial condition ξ(0) = Y 0 , then ξ(s) ∈ M 2 for all s. Hence, W ∈ M 2 . It is easily shown that if M ∈ M 2 is invertible, then M −1 ∈ M 2 . By Equation (19), it follows that K ∈ M , and so we know that by the above result about one of the eigenvectors of W . Hence, the remaining three eigenvalues of K are simply those of the central 2 × 2 matrix together with e.
Remark: Owing to the decomposition mentioned earlier, the position of e in the matrix K indicates that it should be an eigenvalue corresponding to the behavior of the orbit in the subspace A, along with the trivial eigenvalue −1 from the (1, 1) position. Hence, computing the linear stability of the 4DF orbit automatically gives the stability of the 2DF orbit.
The coordinate changes that we used to regularize the system did not "factor out" the angular momentum (as the polar symplectic transformation does), and so we expect that K will have another trivial eigenvalue ±1 corresponding to this integral. To demonstrate this, let γ(s) again represent the periodic orbit with period T as above, and definev(s) = ∇A(γ(s)). As in [10], p. 134, Lemma 7, we considerv(s) as a left eigenvector of a particular matrix related to the monodromy matrix. We find that Since If φ(s, z) represents the solution to the system of linearized differential equations with initial condition z, we know that A(φ(s, z)) = A(z).
This additional eigenvector and eigenvalue gives us more information about the structure of W , and hence, of K. We know that We readily computev(0)Y 0 = 1 2 (0, −ζ 4 , √ 8, 0, 0, 0, 0, 0). From this, we know that Since K ∈ M , this requires that the additional −1 eigenvalue comes from the central 2 × 2 block in K. Furthermore, this imposes some relations on the entries a, b, c, d in (23). In particular, Remarks: (i) Since K is real-valued, this result, along with other results about the form of K, force all of the eigenvalues of K to be real numbers. (ii) This analysis is an improvement over work done in [1], in which the −1 eigenvalue corresponding to angular momentum showed up numerically but could not be factored out a priori. This is most likely due to the relative simplicity of the rhomboidal orbit.

Results
We numerically obtain the matrix W by a numerical integration of the linearized systems (17) and the initial conditions computed in Section 2.3. The values of a, d, and e in the matrix K, as given in (23), are readily computed using (20) and the computed value of ζ 4 from Section 2.3. The resulting eigenvalue calculations are represented in Figures 11 and 12.
We note first that if we restrict to the subspace A, then the eigenvalue of K given by the (4, 4) entry corresponds to linear stability of the 2DF orbit. This value stays in the interval [−1, 1] for m ∈ (0.01 + ε, 1], giving the linear stability result claimed earlier. Based on the additional results of the numerical calculations, we conclude that the 4DF rhomboidal orbit is linearly unstable, hence unstable, for all m except for a small interval about m = 0.4. Additionally, there are three values of m for which we establish only spectral stability, due to repeated eigenvalues on the unit circle. Roberts' argument (see [14]) demonstrates that each of the computed eigenvalues of K in [−1, 1] correspond to the real part of a Fig. 11 A plot of the nontrivial eigenvalue of the central 2 × 2 submatrix of K as a function of m. This eigenvalue crosses the y-axis for some value of m ≈ 0.4. d Fig. 12 A plot of the nontrivial eigenvalues of K as functions of m. The thicker line represents the (4, 4) entry of K. The thinner line is the same curve as plotted in Figure 11 with the value at m = 0.4 emphasized.
square root of an eigenvalue on the complex unit circle. Accordingly, the value of m = m 1 where the two curves in Figure 12 cross is a point with duplicated eigenvalues, hence only spectral stability. Similarly, the value m = m 0 where the curve in Figure 11 crosses the x-axis gives spectral stability, as (±i) 2 = 1. A third point occurs where cos(2α(m)) = cos(2β(m)), where α(m) and β(m) are the the two curves plotted in Figure 12.
In order to obtain a more precise intervals of mass values for which the orbit is linearly stable (excluding the above-mentioned m i ), the initial conditions for mass values m = 0.39, 0.391, ..., 0.409, 0.41 were obtained using the same trigonometric polynomial approximation/optimization as used in Section 2.3. The same linear stability calculations demonstrate that the 4DF rhomboidal orbit is linearly stable for m contained in some subinterval of (0.395, 0.401). In other words, the orbit was found to be linearly unstable for m = 0.395 and m = 0.401, but linearly stable for all computed values in between, with the exclusion of the three critical mass values.
In [5], Bounemoura shows that in an n-dimensional Hamiltonian system, orbits beginning close to a linearly stable invariant torus generically remain "close" to the invariant torus for a super-exponential amount of time, eventually drifting away. The theory can also be ap-plied to other cases, such as elliptic fixed points of maps. This analysis leads us to believe that, even though we have linear stability for some open interval containing m = 0.4, the orbit is likely still be unstable. Numerical perturbations off of the invariant subspace A give evidence that this is the case.