Stability of Broucke's Isosceles Orbit

TL;DR

This paper analyzes the linear stability of Broucke's isosceles orbit with specific mass configurations, extending previous results and showing stability conditions through analytical and numerical methods.

Contribution

It extends Yan's results to a new mass configuration, introduces a regularization method for binary collisions, and applies an analytical stability analysis using Roberts' method.

Findings

Four-degrees-of-freedom orbit is unstable except for specific mass intervals.

Two-degrees-of-freedom subset exhibits wider stability interval.

Stability analysis reduces to computing three entries of a 4x4 matrix.

Abstract

We extend the result of Yan to Broucke's isosceles orbit with masses $m_1$, $m_1$, and $m_2$ with $2m_1 + m_2 = 3$. Under suitable changes of variables, isolated binary collisions between the two mass $m_1$ particles are regularizable. We analytically extend a method of Roberts to perform linear stability analysis in this setting. Linear stability is reduced to computing three entries of a $4 \times 4$ matrix related to the monodromy matrix. Additionally, it is shown that the four-degrees-of-freedom setting has a two-degrees-of-freedom invariant set, and linear stability results in the subset comes `for free' from the calculation in the full space. The final numerical analysis shows that the four-degrees-of-freedom orbit is linearly unstable except for the interval $0.555 < m_1 < 0.730$, whereas the two-degrees-of-freedom orbit is linearly stable for a much wider interval.

Figures

Equations

$$m_i{\ddot{q}}_i=\sum _{i\ne j}\frac{m_i{m}_j(q_j-q_i)}{|q_i-q_j{|}^3},$$

$$A_0=q_1{p}_2-q_2{p}_1+q_3{p}_4-q_4{p}_3+q_5{p}_6-q_6{p}_5.$$

$$q_5=-\frac{m_1(q_1+q_3)}{m_2},q_6=-\frac{m_1(q_2+q_4)}{m_2},$$

$$p_5=-(p_1+p_3),p_6=-(p_2+p_4).$$

$$K_0=\frac{1}{2}\left(\frac{p_1^2}{m_1}+\frac{p_2^2}{m_1}+\frac{p_3^2}{m_1}+\frac{p_4^2}{m_1}+\frac{(p_1+p_3{)}^2}{m_2}+\frac{(p_2+p_4{)}^2}{m_2}\right)$$

$$U_0$$

$$+\frac{m_1{m}_2}{{\left({\left(q_1+\frac{m_1(q_1+q_3)}{m_2}\right)}^2+{\left(q_2+\frac{m_1(q_2+q_4)}{m_2}\right)}^2\right)}^{\frac{1}{2}}}$$

$$+\frac{m_1{m}_2}{{\left({\left(q_3+\frac{m_1(q_1+q_3)}{m_2}\right)}^2+{\left(q_4+\frac{m_1(q_2+q_4)}{m_2}\right)}^2\right)}^{\frac{1}{2}}}.$$

$$F_1=\left(\frac{u_1+u_3}{2}\right)p_1+\left(\frac{u_2+u_4}{2}\right)p_2+\left(\frac{u_3-u_1}{2}\right)p_3+\left(\frac{u_4-u_2}{2}\right)p_4.$$

$$q_1$$

$$q_2$$

$$q_3$$

$$q_4$$

$$p_1=v_1+v_3,p_2=v_2+v_4,p_3=v_3-v_1,p_4=v_4-v_2.$$

$$K_1=\frac{1}{m_1}(v_1^2+v_2^2+v_3^2+v_4^2)+\frac{2}{m_2}(v_3^2+v_4^2)$$

$$U_1$$

$$+\frac{m_1{m}_2}{{\left({\left(\frac{1}{2}{u}_1+\mu u_3\right)}^2+{\left(\frac{1}{2}{u}_2+\mu u_4\right)}^2\right)}^{\frac{1}{2}}}$$

$$+\frac{m_1{m}_2}{{\left({\left(\mu u_3-\frac{1}{2}{u}_1\right)}^2+{\left(\mu u_4-\frac{1}{2}{u}_2\right)}^2\right)}^{\frac{1}{2}}}.$$

$$A_1=u_1{v}_2-u_2{v}_1+2\mu (u_3{v}_4-u_4{v}_3).$$

$$F_2=v_1(Q_1^2-Q_2^2)+2v_2{Q}_1{Q}_2+v_3{Q}_3+v_4{Q}_4.$$

$$u_1$$

$$u_2$$

$$u_3$$

$$u_4$$

$$v_1=\frac{Q_1{P}_1-Q_2{P}_2}{2(Q_1^2+Q_2^2)},v_2=\frac{Q_2{P}_1+Q_1{P}_2}{2(Q_1^2+Q_2^2)}.$$

$$K_2=\frac{P_1^2+P_2^2}{4m_1(Q_1^2+Q_2^2)}+\left(P_3^2+P_4^2\right)\left(\frac{1}{m_1}+\frac{2}{m_2}\right)$$

$$U_2$$

$$+\frac{m_1{m}_2}{\sqrt{\frac{1}{4}(Q_1^2+Q_2^2{)}^2+\mu \mathbf{Q}+{\mu }^2(Q_3^2+Q_4^2)}}$$

$$+\frac{m_1{m}_2}{\sqrt{\frac{1}{4}(Q_1^2+Q_2^2{)}^2-\mu \mathbf{Q}+{\mu }^2(Q_3^2+Q_4^2)}},$$

$$A_2=\frac{1}{2}\left(Q_1{P}_2-Q_2{P}_1\right)+2\mu (Q_3{P}_4-Q_4{P}_3).$$

$$\Gamma$$

$$-\frac{m_1{m}_2\left(Q_1^2+Q_2^2\right)}{\sqrt{\frac{1}{4}(Q_1^2+Q_2^2{)}^2+\mu \mathbf{Q}+{\mu }^2(Q_3^2+Q_4^2)}}$$

$$-\frac{m_1{m}_2\left(Q_1^2+Q_2^2\right)}{\sqrt{\frac{1}{4}(Q_1^2+Q_2^2{)}^2-\mu \mathbf{Q}+{\mu }^2(Q_3^2+Q_4^2)}}$$

$$-m_1^2-E(Q_1^2+Q_2^2).$$

$$q_1=-q_3\text{, }{q}_2=q_4\text{, }{p}_1=-p_3\text{, and }{p}_2=p_4.$$

$$Q_2=Q_3=P_2=P_3=0$$

$$\dot{Q}_{2}\big{|}_{\mathcal{A}}=\dot{Q}_{3}\big{|}_{\mathcal{A}}=\dot{P}_{2}\big{|}_{\mathcal{A}}=\dot{P}_{3}\big{|}_{\mathcal{A}}=0,$$

$$J=\left[\begin{array}{cc}O& I\\ -I& O\end{array}\right],$$

Full text

1. Introduction

Mathematically, the study of determining the motion of $n$ point masses in space whose motion is governed by Newton’s gravitational law is known as the Newtonian $n$-body problem. Notationally, if $\{q_{1},q_{2},...,q_{n}\}$ represent the positions of the bodies in $\mathbb{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

[TABLE]

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\neq j$. In the equations governing the motion, this results in a zero denominator in one or more terms in the sum. Under certain conditions, these collisions can be regularized and the solutions can be continued past collision.

In this paper, we will study Broucke’s isosceles triangle orbit, originally presented in [6]. This is among the earliest-introduced orbits featuring collisions in the planar setting. The orbit features regularizable collisions between two of the bodies, while the third oscillates along the vertical axis.

The chief aim of this paper is to give stability results, including linear stability, for the orbit being considered in both a four-degrees-of-freedom (4DF) and a two-degrees-of-freedom (2DF) setting at the same time. With our choice of coordinates in the 4DF setting, the linearized phase space of the regularized equations has an elegant decomposition into two invariant subspaces. The linear stability analysis in the 2DF setting corresponds to one of these subspaces. Additionally, we are able to use analytical techniques to reduce the linear stability analysis of the numerical calculation of three entries of a $4\times 4$ matrix $K$ related to the monodromy matrix. This type of stability calculation based on an invariant-subspace decomposition was done earlier in [1]. A suitable coordinate transformation reveals this same structure in the Broucke’s setting.

Schubart [18] was one of the first to study periodic orbits with regularizable 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 [7] in 1977. Analytic existence of the equal-outer-mass orbit was established independently by Venturelli [22] and Mockel [13], both in 2008. Shibiyama [19] recently demonstrated the existence of the arbitrary-mass version. The study of linear stability of Schubart’s orbit was performed by Hietarinta and Mikkola [8] in 1993.

Sweatman found a Schubart-like collinear four-body symmetric orbit in 2002 [20], and later studied its linear stability [21]. 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 [15].

Apart from the Broucke’s orbit, other planar orbits with regularizable collisions have also been studied. Among the first of these is the rhomboidal four-body orbit, which features two pairs of bodies: one pair on the $x$-axis and the second on the $y$-axis. 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 [10] for equal masses, [19] for symmetric masses). Additionally, Yan showed that for equal masses, the orbit is linearly stable. This was followed up by work in [1], in which linear stability was shown for a wide range of mass ratios.

Other planar orbits with singularities have also been studied. A planar four-body orbit featuring simultaneous binary collisions was described in [16]. The orbit was shown to be linearly stable in [5]. It was later shown that this orbit could be numerically extended to symmetric masses in [3] (see also [4]), and linear stability for this extension was shown for an interval of certain mass ratios in [2].

More generally, analytic existence of large families of orbits with two degrees of freedom and regularizable singularities was recently proven by Shibayama in [19] and Martinez in [10].

The remainder of the paper is as follows: In Section 2, we describe the orbit and some of its properties. In 2.1, we give formal notation describing the orbit. We also perform the coordinate transformations that regularize the collisions. Section 2.2 gives an analytic proof of the existence of the orbit in the regularized 4DF setting.

Section 3 deals with the linear stability of the orbit. In 3.1, we review some of the basic properties of linear stability. Next, 3.2 describes the symmetries of the orbit, which are needed to perform the analysis in 3.3. Section 3.4 describes all of the remaining linear stability analysis that can be done before any numerical work, including the decomposition into two invariant sets mentioned earlier.

Finally, in Section 4, we present the results obtained from the numerical calculations of stability.

2. Broucke’s Isosceles Orbit

2.1. Setting and Regularization

We denote the position of the three bodies with the variables $\{q_{i}\}$. The two bodies at $(q_{1},q_{2})$ and $(q_{3},q_{4})$ are given mass $m_{1}$, and the third body has mass $m_{2}$, where both $m_{1}$ and $m_{2}$ are non-negative real numbers. Also, $p_{i}=m_{1}\dot{q}_{i}$ for $i=1,2,3,4$, and $p_{i}=m_{2}\dot{q}_{i}$ for $i=5,6$. (See Figure 1.)

The orbit maintains the shape of an isosceles triangle for all time while alternating between collisions of the two bodies of mass $m_{1}$, as shown below. Beginning at collision on the positive $y$-axis (shown by the lightest-colored circles in Figure 2), the two bodies of mass $m_{1}$ move outward in a symmetric fashion as the mass $m_{2}$ body moves upward along the negative $y$-axis. This eventually leads to a collinear configuration along the $x$-axis (darkest circles). The motion is then continued in reverse with the vertical coordinates negated, leading to a collision of the two mass $m_{1}$ bodies on the negative $y$-axis. After this collision, the bodies move back along the same paths leading back to a collision of the two mass $m_{1}$ bodies on the positive $y$-axis. Throughout, the mass $m_{2}$ body oscillates along the $y$-axis. Tracing out the full path the bodies take yields a shape resembling the letter $\Phi$.

For simplicity, we will assume that $2m_{1}+m_{2}=3$ throughout. Thus, the equal-mass case occurs when $m_{1}=m_{2}=1$.

The angular momentum is given by

[TABLE]

Assuming the center of mass is located at the origin and net momentum is zero yields

[TABLE]

[TABLE]

The Hamiltonian function in this setting, which we will denote $H_{0}$, is given by $H_{0}=K_{0}-U_{0}$, where

[TABLE]

and

[TABLE]

We wish to regularize the collision between the two particles of mass $m_{1}$. We define a first transformation of variables $(q_{i},p_{i})\leftrightarrow(u_{i},v_{i})$ by means of the generating function

[TABLE]

Then, with $q_{i}=\partial F_{1}/\partial p_{i}$ and $v_{i}=\partial F_{1}/\partial u_{i}$, we have

[TABLE]

Solving for each $p_{i}$ gives

[TABLE]

Substituting, the new Hamiltonian $H_{1}$ is given by $K_{1}-U_{1}$, where

[TABLE]

and setting $\mu=\frac{1}{2}+\frac{m_{1}}{m_{2}}$, we have

[TABLE]

The angular momentum after this transformation is

[TABLE]

Remark: This transformation “codes in” some of the symmetry present in the orbit. In particular, when $u_{2}=0$ we have that $q_{2}-q_{4}=0$, which places the two bodies of mass $m_{1}$ on the same horizontal line. Similarly, if $u_{3}=0$ then $q_{1}+q_{3}=0$, which places those same bodies symmetrically across the $y$-axis. Similar results are present for $v_{2}0$ and $v_{3}=0$. As these behaviors are expected along the orbit, these are natural targets for the coordinate transform. This also gives a similar setting to what is present in the rhomboidal problem discussed in [1], wherein “symmetry-preserving zeros” facilitated analysis of the orbit.

The second transformation $(u_{i},v_{i})\leftrightarrow(Q_{i},P_{i})$ is generated by

[TABLE]

With $u_{i}=\partial F_{2}/\partial v_{i}$ and $P_{i}=\partial F_{2}/\partial Q_{i}$, we have

[TABLE]

Solving for $v_{1}$ and $v_{2}$ gives

[TABLE]

Making the substitution gives the new Hamiltonian $H_{2}=K_{2}-U_{2}$, with

[TABLE]

and (after some simplification)

[TABLE]

where $\mathbf{Q}=(Q_{1}^{2}-Q_{2}^{2})Q_{3}+2Q_{1}Q_{2}Q_{4}$.

The angular momentum is now

[TABLE]

As a final step, we create the regularized Hamiltonian $\Gamma=\frac{dt}{ds}(H_{2}-E)$, where $\frac{dt}{ds}=Q_{1}^{2}+Q_{2}^{2}$. Then

[TABLE]

It is easy to verify that when $Q_{1}^{2}+Q_{2}^{2}=0$, the condition $\Gamma=0$ forces $P_{1}^{2}+P_{2}^{2}=4m_{1}^{3}$. Hence, the binary collisions between the two equal-mass particles has been regularized. Attempting to regularize the other binary collisions is not needed for this orbit, and is in general a difficult problem (see [14]). Also, regularization of the triple collision where all $Q_{i}=0$ (corresponding to total collapse of the system at the origin) is generally not possible (see [11]).

2.2. Analytic Existence of the Orbit

In [24], existence of Broucke’s orbit was established in a reduced, two-degree-of-freedom setting with $m_{1}=1$. We demonstrate the existence of the orbit in the larger, four-degree-of-freedom setting established here.

Theorem 2.1**.**

Broucke’s isosceles orbit described in Section 2.1 analytically exists for the Hamiltonian system given by $\Gamma$.

Proof.

Let $\mathcal{A}$ denote the set where

[TABLE]

This corresponds to the setting used in [24], wherein a symmetry across the $y$-axis is forced for all time. In that paper, it was assumed that $m_{2}$ (corresponding to the body running along the $y$-axis) was arbitrary, and $m_{1}$ (corresponding to the mass of the colliding pair) was fixed to be 1. With a suitable re-scaling of mass ratios, the proof in [24] shows that the orbit exists in this setting.

In the set $\mathcal{A}$, we have $u_{2}=u_{3}=v_{2}=v_{3}=0$ in the first coordinate transformation. In turn, $u_{3}=v_{3}=0$ forces $Q_{3}=P_{3}=0$. The condition $u_{2}=0$ requires that either $Q_{1}=0$ or $Q_{2}=0$. Additionally, since $v_{2}=0$, we have that $P_{2}=2v_{1}Q_{2}$ and $P_{1}=2v_{1}Q_{1}$. In keeping with our analysis from [1], we will make the choice to set $Q_{2}=0$. Then, when $\mathcal{A}$ holds, we have

[TABLE]

Furthermore, we have

[TABLE]

so $\mathcal{A}$ is invariant. Since analytic existence of the 2DF was given in [24], analytic existence in the 4DF setting follows from the invariance of $\mathcal{A}$ following a rescaling of the masses. ∎

As a further consequence of the invariant set $\mathcal{A}$, initial conditions for the periodic orbit in the 2DF setting automatically give the initial conditions for the 4DF setting. This is very useful numerically as it reduces the number of calculations required to find initial conditions.

3. Linear Stability Analysis

The material in this section is very similar to that of [1], which in turn draws on material in [17]. We repeat it here for completeness.

3.1. Linear Stability Theory

We will begin this section by briefly reviewing the basic definitions for linear stability. Let $\Gamma$ be a smooth function defined on an open subset of $\mathbb{R}^{8}$. Let $\gamma(s)$ be a $T$-periodic solution of the system $z^{\prime}=JD\Gamma(z)$, where ${}^{\prime}=d/ds$,

[TABLE]

and $I$ and $O$ are the $4\times 4$ identity and zero matrices, respectively. If $X(s)$ is the fundamental matrix solution of the linearized equations

[TABLE]

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 $\gamma$, and therefore determine the linear stability of $\gamma$. In particular, $\gamma$ is spectrally stable if all of its characteristic multipliers lie on the unit circle, and $\gamma$ is linearly stable if it is spectrally stable and $X(T)$ is semisimple apart from trivial eigenvalues.

Additionally, if $Y$ is the fundamental matrix solution of

[TABLE]

for some invertible matrix $Y_{0}$, then by definition of $X(s)$, $Y(s)=X(s)Y_{0}$, implying $X(T)=Y(T)Y_{0}^{-1}$. Then we have

[TABLE]

and so $X(T)$ and $Y_{0}^{-1}Y(T)$ are similar, and stability can be determined by the eigenvalues of either.

Linear stability is typically established by numerical integration. Some elegant techniques for simplifying the numerical work by using the symmetries of an orbit were presented by Roberts in [17]. The relevant theory will be reviewed in Section 3.3.

3.2. Symmetries of the Orbit

In order to utilize the techniques developed by Roberts, it is necessary to first identify the symmetries of the periodic orbit.

Lemma 3.1**.**

The regularized Hamiltonian $\Gamma$ has symmetry group isomorphic to the Klein four group.

Proof.

Let

[TABLE]

and define the block matrix

[TABLE]

where [math] represents the $2\times 2$ identity matrix. Then we have

[TABLE]

Hence, $S$ and $-S$ generate a group isomorphic to the Klein four group. For fixed values of $m$ and $E$, we have

[TABLE]

so $\pm S$ generate a Klein-four symmetry group for $\Gamma$ as well.

∎

Symmetries for $\Gamma$ also help to determine symmetries for the periodic orbit.

Theorem 3.2**.**

Let $\gamma$ be a solution to the Hamiltonian system for some $m_{1}$, $m_{2}$ and $E<0$ such that

[TABLE]

and

[TABLE]

for some values of $\zeta_{1}$, $\zeta_{4}$, and $\zeta_{8}$. (In other words, $\gamma(0)$ corresponds to collision of the two bodies of mass $m_{1}$ somewhere along the $y$-axis away from the origin, and $\gamma(s_{0})$ corresponds to a collinear arrangement of the three bodies along the $x$-axis.) Then $\gamma$ extends to a $T=4s_{0}$-periodic orbit, wherein $S$ and $-S$ are symmetries of the orbit.

Remark: Although the $\Gamma=0$ condition forces a relationship between $\zeta_{1}$ and $\zeta_{8}$, the most relevant detail is that both quantities are non-zero.

Proof.

Note that if $\gamma(s)$ is a $T$-periodic solution to the regularized equations of motion resulting from the Hamiltonian $\Gamma$, a standard calculation shows that both $-S\gamma(T/2-s)$ and $S\gamma(T-s)$ are solutions as well. Existence and uniqueness of solutions then imply that

[TABLE]

for all $s$. Hence the symmetry group for Broucke’s isosceles orbit is isomorphic to the Klein four group, with $S$ and $-S$ as generators. ∎

3.3. Stability Reduction using Symmetries

The following can be found in [17]:

Lemma 3.3**.**

Suppose $\gamma(s)$ is a $T$-periodic solution of a Hamiltonian system with Hamiltonian $\Gamma$ and a time-reversing symmetry $S$ such that:

1. (i)

For some $n\in\mathbb{N}$, $\gamma(-s+T/N)=S(\gamma(s))$ for all $s$; 2. (ii)

$\Gamma(Sz)=\Gamma(z)$; 3. (iii)

$SJ=-JS$; 4. (iv)

$S$* is orthogonal.*

Then the fundamental matrix solution $X(s)$ satisfies

[TABLE]

Note that the matrix $S$ given in (3.3) satisfies all the required hypotheses.

Corollary 3.4**.**

Under the same hypotheses,

[TABLE]

Corollary 3.5**.**

If $Y(s)$ is the fundamental matrix solution to (3.2), then

[TABLE]

and

[TABLE]

Similar results for time-preserving symmetries are also presented in [17], 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 3.5 with $N=2$, $S$ as defined in (3.3), and noting that $S^{T}=S$ yields

[TABLE]

Similarly, if $N=1$, since $S^{2}=I$, we get

[TABLE]

This yields

[TABLE]

with $W=Y_{0}^{-1}SY_{0}Y(T/4)^{-1}SY(T/4)$. Hence, in order to analyze the stability of the orbit, we need only compute the entries of $Y$ along a quarter of the orbit.

Again, from [17]:

Lemma 3.6**.**

For a symplectic matrix $W$, suppose there is a matrix $K$ such that

[TABLE]

Then the eigenvalues of $W$ lie on the complex unit circle if and only if all of the eigenvalues of $K$ are real and lie in the interval $[-1,1]$

A matrix whose eigenvalues all lie on the complex unit circle is often referred to as stable.

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

[TABLE]

then setting

[TABLE]

yields $-Y_{0}^{-1}SY_{0}=\Lambda$. (The lines here are provided for ease in reading. Much of our later analysis will involve breaking $8\times 8$ matrices down into $4\times 4$ blocks.) Furthermore, it is easy to check that $Y_{0}$ is both orthogonal and symplectic. If we set $D=-B^{-1}SB$ for $B=Y(T/4)$, we then have

[TABLE]

Also, since $\Lambda^{2}=D^{2}=I$, we know immediately that

[TABLE]

Since $B=Y(T/4)$ is symplectic, setting

[TABLE]

gives

[TABLE]

where $L_{1}$, $L_{2}$, and $K$ are $4\times 4$ matrices satisfying $L_{1}=B_{4}^{T}S_{1}B_{2}+B_{2}^{T}S_{1}B_{4}$, $L_{2}=B_{3}^{T}S_{1}B_{1}-B_{1}^{T}S_{1}B_{3}$, and $K=-B_{2}^{T}S_{1}B_{2}-B_{1}^{T}S_{1}B_{4}$. Thus,

[TABLE]

Similarly, we find that

[TABLE]

Thus, we have

[TABLE]

Remark: The given matrix $Y_{0}$ in (3.4) is not unique. Different choices of $Y_{0}$ are possible, but our particular choice is helpful for much of our later analysis. It is also worth noting that our choice of $Y_{0}$ is independent of the values of $m_{1}$ and $m_{2}$ for this orbit, which is not always true (see [2].)

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

[TABLE]

Using $W=\Lambda D$ for $D=-B^{-1}SB$ and the relation $-SJ=JS$, we find

[TABLE]

Directly computing $\Lambda J$ and using the block form of $B$, we find that

[TABLE]

Define $\text{col}_{i}(-SJB)$ to be the $i$th column of the matrix $-SJB$. Then we have $\text{col}_{i}(-SJB)=-SGc_{i}$ where $c_{i}$ is the $i$th column of $B$. Using the above two formulas, this implies that the $(i,j)$ entry of $W$ is given by $-c_{i}^{T}SJC_{j}$. Equation (3.5) shows that the $(i,j)$ entry of $K$ is the $(i+4,j+4)$ entry of $W$. Hence,

[TABLE]

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.

3.4. Entries of $K$ from Invariant Quantities

Before any numerical work is done, we can determine many of the values of entries of $K$ by using properties of the orbit. We first introduce some notation to simplify the analysis. Let $\mathcal{M}$ denote the set of matrices of the form

[TABLE]

where all of the listed $m_{ij}\in\mathbb{R}$ ($m_{ij}=0$ is permitted). Further, let $\mathcal{M}_{2}$ denote the set of $8\times 8$ matrices whose $4\times 4$ blocks are in $\mathcal{M}$. That is to say, $\mathcal{M}_{2}$ consists of matrices of the form

[TABLE]

where each of the $M_{i}\in\mathcal{M}$. It is easy to check that $\mathcal{M}$ forms a ring with the standard definitions of matrix addition and multiplication. Furthermore, each element of $\mathcal{M}$ with nonzero determinant has its inverse in $\mathcal{M}$ as well. The same ring structure exists for $\mathcal{M}_{2}$.

The following two lemmas will help establish an important theorem about the form of $K$:

Lemma 3.7**.**

If $M\in\mathcal{M}_{2}$, then the system of differential equations given by

[TABLE]

and initial conditions

[TABLE]

has solutions of the form

[TABLE]

Proof.

We verify that $M\eta$ has the proper form. Note that

[TABLE]

Hence, the zeros in the 2nd, 3rd, 6th, and 7th are preserved under multiplication by $M$. So

[TABLE]

is a solution of $\eta^{\prime}=M\eta$. Existence and uniqueness of solutions implies that $\eta(s)$ is the only solution of the system. ∎

Lemma 3.8**.**

If $M\in\mathcal{M}_{2}$, then the system of differential equations given by

[TABLE]

and initial conditions

[TABLE]

has solutions of the form

[TABLE]

Theorem 3.9**.**

If $Y_{0}$ is given by (3.4), then $K\in\mathcal{M}$.

Proof.

Using a computer algebra system, we find that the matrix $JD^{2}\Gamma$ is of the form

[TABLE]

Here, the zeros denote entries for which the mixed partials evaluate to zero identically, and the entries denoted $\omega$ are entries for which the mixed partials evaluate to zero assuming the conditions in $\mathcal{A}$ which hold along the periodic orbit $\gamma(s)$. Under such conditions, we have $JD^{2}\Gamma\in\mathcal{M}_{2}$. If $\xi(0)\in\mathcal{M}_{2}$, then each of the columns of $\xi(0)$ has the same form as in either Lemma 3.7 or 3.8. Hence, the solution to the system of linearized equations given by (3.2) satisfies $\xi(s)\in\mathcal{M}_{2}$ for all $s$. Since $Y_{0}\in\mathcal{M}_{2}$, $\xi(s)\in\mathcal{M}_{2}$ for all $s$. Lastly, because $M\in\mathcal{M}_{2}$ implies $M^{-1}\in\mathcal{M}_{2}$ when $M^{-1}$ exists, Equation (3.5), gives $K\in\mathcal{M}$ as claimed. ∎

Remarks:

1. (i)

The structure of $\mathcal{M}_{2}$ very nicely decomposes phase space into a direct sum of $\mathcal{A}=\{Q_{2}=Q_{3}=P_{2}=P_{3}=0\}$ and $\mathcal{A}^{\perp}=\{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 $\mathcal{A}^{\perp}$ is appropriate in that $\mathcal{A}^{\perp}$ and $\mathcal{A}$ are orthogonal complements in $\mathbb{R}^{8}$. The two are also skew-orthogonal: if $a_{1}\in\mathcal{A}$ and $a_{2}\in\mathcal{A}^{\perp}$, then $a_{1}^{T}Ja_{2}=0$. 2. (ii)

Matrices of the form $\mathcal{M}$ and $\mathcal{M}_{2}$ are similar to the diamond product discussed in [9]. Specifically, $\Sigma^{-1}M\Sigma=\mathcal{A}_{1}\Diamond\mathcal{A}_{2}$ for some matrices $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$, where $M\in\mathcal{M}_{2}$ and $\Sigma$ is a (not necessarily unique) permutation matrix. Furthermore, one of the $\mathcal{A}_{i}$ (depending on choice of $\Sigma$) corresponds to the 2DF setting. 3. (iii)

The particular choice of $Y_{0}$ given in (3.4) is important for this argument.

In light of Theorem 3.9, we need only to find eight of the entries of the matrix $K$. We can, in fact, reduce this number further by using invariant properties of the orbit $\gamma(s)$. As is well-known, invariant quantities of the $n$-body problem are center of mass, net momentum, angular momentum, and energy (the Hamiltonian itself). Each of these correspond to trivial eigenvalues of the monodromy matrix. The center of mass and net momentum were “factored out” by our choice of coordinates at the beginning. The remaining two invariant quantities will be used to reduce the number of entries of $K$ needed to find its eigenvalues.

Theorem 3.10**.**

The matrix $K$ has a right eigenvector $[1,0,0,0]^{T}$, corresponding to eigenvalue $-1$.

Proof.

Let $v=Y_{0}^{-1}\gamma^{\ {}^{\prime}}(0)/||\gamma^{\ {}^{\prime}}(0)||$ or, equivalently, $Y_{0}^{T}\gamma^{\ {}^{\prime}}(0)/||\gamma^{\ {}^{\prime}}(0)||$. By Corollary 3.5, since $Y_{0}$ is orthogonal and $S$ is symmetric, we have

[TABLE]

Since $\gamma^{\ {}^{\prime}}(s)$ is a solution of $\dot{\xi}=JD^{2}\Gamma(\gamma(s))\xi$ and $\gamma^{\ {}^{\prime}}(0)=Y(0)Y_{0}^{-1}\gamma^{\ {}^{\prime}}(0)=Y(0)v$, we also know that $\gamma^{\ {}^{\prime}}(s)=Y(s)Y_{0}^{-1}\gamma^{\ {}^{\prime}}(0)=Y(s)v$. This implies

[TABLE]

By the symmetry $\gamma(s)=-S\gamma(T/2-s)$, we also have $\gamma^{\ {}^{\prime}}(s)=S\gamma^{\ {}^{\prime}}(T/2-s)$. Setting $s=0$ in this setting tells us that $\gamma^{\ {}^{\prime}}(0)=S\gamma^{\ {}^{\prime}}(T/2)$. Since

[TABLE]

for some real number $\omega$, we have $-S\gamma^{\ {}^{\prime}}(0)=\gamma^{\ {}^{\prime}}(0)$. Thus

[TABLE]

Combining (3.7) and (3.8) gives $Wv=-v$, and so $-1$ is an eigenvalue of $W$ with eigenvector $v$. By definition, we have that

[TABLE]

From the form of $W$, this implies that

[TABLE]

So $[1,0,0,0]^{T}$ is an eigenvector of $K$ with eigenvalue $-1$. Consequently, the first column of $K$ must be $[-1,0,0,0]^{T}$.

∎

Combining the results of the previous two theorems gives

[TABLE]

Remark: Owing to the decomposition of linearized phase space into two invariant subspaces and the ordering of the coordinates, the position of $e$ in the matrix $K$ indicates that it should be an eigenvalue corresponding to the behavior of the orbit in $\mathcal{A}$. This eigenvalue, along with the trivial eigenvalue $-1$ from the $(1,1)$ position, completely classify linear stability in the 2DF setting. Hence, computing the linear stability of the 4DF orbit in the chosen coordinates automatically gives the stability of the 2DF orbit. (The results will be discussed further in Section 4.2.)

We can make use of the final invariant quantity, angular momentum, to further simplify our calculations. This is an extension of Roberts’ method from [17], in which coordinate transformations “factor out” the angular momentum before linearization is performed. In the following theorem, we are able to show that this invariant quantity can be used to simplify linear stability calculations after linearization.

Theorem 3.11**.**

The matrix $K$ has a left eigenvector $\nabla A(\gamma(0))Y_{0}$ with eigenvalue $-1$.

Proof.

This proof is based on ideas given in [12], p. 134, Lemma 7. Define $\hat{v}(s)=\nabla A(\gamma(s))$, where $A$ represents the regularized angular momentum given in (2.1). Then

[TABLE]

Since

[TABLE]

we know that

[TABLE]

Let $\phi(s,z)$ be the general solution to the system of regularized differential equations with initial condition $z$. Then

[TABLE]

Differentiating with respect to $z$ gives

[TABLE]

or, equivalently

[TABLE]

where $X(s)$ is the fundamental matrix solution. Setting $s=T/2$ and substituting $X(T/2)=Y_{0}(Y_{0}^{-1}Y(T/2))Y_{0}^{-1}$ gives

[TABLE]

and so

[TABLE]

By the symmetry of the orbit, $\gamma(T/2)=-\gamma(0)$, which gives

[TABLE]

and therefore

[TABLE]

Hence $\hat{v}(0)Y_{0}$ is a left eigenvector for $W$ with eigenvalue $-1$. ∎

Remarks:

1. (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. 2. (ii)

This analysis is an improvement over work done in [2], in which the $-1$ eigenvalue corresponding to angular momentum showed up numerically but could not be factored out a priori.

To find the eigenvector, we readily compute $\hat{v}(0)Y_{0}=(0,-2\mu\zeta_{4},m_{1}^{3/2},0,0,0,0,0)$. From this, we know that

[TABLE]

Since $K\in\mathcal{M}$, this requires that the additional $-1$ eigenvalue comes from the central $2\times 2$ block in $K$. Furthermore, this imposes some relations on the entries $a,b,c,d$ in (3.9). In particular,

[TABLE]

[TABLE]

Since Theorems 3.10 and 3.11 account for two of the eigenvalues of $K$, we have shown that we only need to find three entries of the matrix $K$ in order to determine the linear stability of the orbit. Specifically, finding $e$ gives the linear stability in the 2DF setting. Combining this result with one element of $\{a,b\}$ and one element of $\{c,d\}$ allows us to determine the second eigenvalue corresponding to linear stability in the 4DF setting.

4. Results

4.1. Initial Conditions of the Orbit

Recall that the initial conditions for the orbit in the regularized setting determined by $\Gamma$ are

[TABLE]

For $m_{1}=0$, the corresponding orbit consists of a body of mass 3 at the origin and two massless particles on the $x$-axis. If $m_{1}=1.5$, the corresponding orbit contains a massless particle at the origin and two bodies of mass 1.5 on the $x$-axis. Note that any vertical motion of the positive-mass particle(s) in either case would cause the center of mass to drift away from the origin, so the orbit must remain collinear for all time. Hence, we expect that at as $m_{1}\to 0^{+}$ or $m_{1}\to 1.5^{-}$, the value of $Q_{4}(0)$ should approach 0.

Setting $E=-1$, and using numerical methods, we find the value of $Q_{4}(0)$. The results of that calculation are shown in Figure 3. Numerical difficulties occurred as $m_{1}\to 1.5^{-}$ (equivalently $m_{2}\to 0^{+}$), likely due to the prevalence of $m_{1}/m_{2}$ terms occurring in the equations of motion. (In fact, numerical methods were not able to find initial conditions for values of $m_{1}\geq 1.47$.)

4.2. Stability in Two Settings

As discussed in Section 2.2, in order to find the initial conditions for the 4DF orbit, we need only to find the initial conditions for the 2DF orbit. This is done numerically for values of $m_{1}=0.005,0.010,0.015,...,1.485,1.490,1.495$.

We numerically obtain the matrix $W$ (hence $K$) by a numerical integration of the linearized systems using the initial conditions discussed previously. The values of $a$, $d$, and $e$ in the matrix $K$, as given in (3.9), are readily computed using (3.6). Knowing these, we are able to determine the eigenvalues of $K$. The results are represented in Figures 4 through 10.

Figures 4 and 5 are used as verification of the numerical method, comparing the proven value $k_{11}=-1$ against the numerically computed value.

In the 2DF setting, the numerical result is that $|e|<1$ for all values of $m_{1}$ computed between $0.125$ through $1.4$ inclusive. However, for values of $m_{1}>1.295$, the curve $y=e(m_{1})$ is no longer monotonically decreasing, which corresponds to the sharp bend in Figure 3. (This sharp bend and the particular numerical difficulties will be further discussed in an upcoming paper.)

Following these calculations, we obtain the following:

Theorem 4.1**.**

Broucke’s orbit in the 2DF setting is linearly stable for $m_{1}$ lying in an open interval containing $[0.125,1.295]$.

For $m_{1}\in[0.01,1.25]$, the computed values of $e$ all lie within the interval $[0.998,1.002]$, with the computed values of $k_{11}$ in the interval $[0.905,1.095]$. Allowing for minor numerical error, we propose the following:

Conjecture 4.2**.**

Broucke’s orbit in the 2DF setting is linearly stable for all $m_{1}$ in an open interval containing $[0.01,0.125]$.

To compute stability in the 4DF setting, we find the second eigenvalue of the central $2\times 2$ block of $K$. From Theorem 3.11, we know that one eigenvalue of this block is $-1$. Hence, using the property that the trace of a matrix is equal to the sum of the eigenvalues, we know the second eigenvalue is given by $a+d+1$, with $a$ and $d$ as in Equation (3.9). This value is plotted in Figures 8 and 9.

Remark: Computing the second eigenvalue could also be done calculating the negative of the determinant of the central $2\times 2$ block. Both were done during the research. Both yielded the same stability results.

These results give the following:

Theorem 4.3**.**

Broucke’s orbit in the 4DF setting is spectrally stable for $m_{1}$ lying in an open interval containing $[0.595,0.715]$.

We also note that there are possible values of $m_{1}$ for which we establish only spectral stability, due to repeated eigenvalues on the unit circle. For simplicity of explanation, let $f(m)$ denote the lighter curve in Figure 10, and $g(m)$ denote the darker curve. Also let $I$ be the interval specified by Theorem 4.3. Roberts’ argument (see [17]) demonstrates that each of the computed eigenvalues of $K$ in $[-1,1]$ correspond to the real part of a square root of an eigenvalue on the complex unit circle. Accordingly, the value of $m_{1}=\mathbf{m_{1}}$ where $f(\mathbf{m_{1}})=g(\mathbf{m_{1}})$ is a point with duplicated eigenvalues, hence only spectral stability. As noted above, $\mathbf{m_{1}}\approx 0.71$. Similarly, a value of $\mathbf{m_{2}}$ satisfying $g(\mathbf{m_{2}})=\pm 1$ would give a pair of eigenvalues of $1$ in $W^{2}$. (If this curve is asymptotic to the line $y=1$ as conjectured, this may not exist. Certainly no such point exists for $m_{1}\in I$) For $\mathbf{m_{3}}$ with $g(\mathbf{m_{3}})=0$, we get $(\pm i)^{2}=-1$ eigenvalues of $W^{2}$. However, this also does not occur for $m_{1}\in I$. Finally, there is a fourth possibility when $\cos(2f(\mathbf{m_{4}}))=\cos(2g(\mathbf{m_{4}}))$, which arises by equating the real parts of $(e^{i\pi\theta_{1}})^{2}=(e^{i\pi\theta_{2}})^{2}$ when $\theta_{1}\neq\theta_{2}$. Within $m_{1}\in I$, this occurs for some value between 0.595 and 0.6, and again at the value of $\mathbf{m_{1}}$ where the two curves meet. Owing to these results, we conclude with the following:

Theorem 4.4**.**

Apart from the two values of $m_{1}$ listed above, Broucke’s orbit is linearly stable in the same interval given by Theorem 4.3.