Periodic and Quasi-Periodic Orbits in the Collinear Four-Body Problem: A Variational Analysis

Abstract

This paper investigated the periodic and quasi-periodic orbits in the symmetric collinear four-body problem through a variational approach. We analyze the conditions under which homographic solutions minimize the action functional. We compute the minimal value of the action functional for these solutions and show that, for four equal masses organized in a linear configuration, these solutions are the minimizers of the action functional. Additionally, we employ numerical experiments using Poincaré sections to explore the existence and stability of periodic and quasi-periodic solutions within this dynamical system. Our results provide a deeper understanding of the variational principles in celestial mechanics and reveal complex dynamical behaviors, crucial for further studies in multi-body problems.

1. Introduction

In our galaxy, most star systems are multiple star systems. These include triple, quadruple, or larger systems, which can be modeled by the four-body problem. Open star clusters, globular clusters, dwarf galaxies, galaxies, galaxy clusters, and superclusters are all examples of multiple star systems. As we know, the binary stars’ systems provide a basic check on stellar structure and evolution theory, and the statistics of binary and multiple star systems present evidence about star formation mechanisms and environmental effects in the galactic gravitational potential and in clusters. It has long been known that the three or more body problem is non-integrable; however, information on the global stability of a gravitational system can be found by studying the regions of the phase space where real motion is possible. The large number of degrees of freedom in the four or more few-body problems make stability criteria of this form difficult to formulate [1,2,3,4,5].

By way of illustration, in hierarchical n-body systems, the masses involved can be divided into subgroups [6]. The relative motion of the subgroups is dominated by their gravitational interaction with one another. The subgroups themselves may be likewise divisible into further subgroups in a recursive fashion. Analytical studies of the stability of hierarchical systems, of three or more bodies, are challenging because of the greater number of variables involved with increasing numbers of bodies and the limitation of just ten integrals that exist in the gravitational n-body problem. The utilization of symmetries and/or neglecting the masses of some of the bodies compared to others can simplify the dynamical problem and enable global analytical stability conditions to be derived. These symmetric and restricted few-body systems with their analytical stability criteria can then provide useful information on the stability of the general few-body system when near symmetry or the restricted situation. Steves et al. [6] developed a general analytical stability criterion for the Caledonian symmetric five-body problem and investigated how variation of the masses affects the stability of the whole system.

The collinear four-body problem is a special case of the general four-body problem in celestial mechanics, in which all four bodies are restricted to moving in a straight line. Compared to the general four-body problem, this has fewer degrees of freedom and is comparatively easier to handle.

Finding equilibrium locations or configurations in the collinear four-body problem can be accomplished by variational methods that minimize the action function. By studying minor perturbations around these points, these methods aid in the determination of the points at which the net forces acting on the bodies balance. Using variational techniques can be useful in investigating the stability or instability of this particular case of the four-body problem.

In [7], Moulton investigated a collinear problem of straight-line solutions for n positive masses. In his work, he arranged n masses along a straight line and determined the conditions under which these masses, positioned at arbitrary collinear points, would remain collinear under appropriate initial velocities. The general case cannot be handled by the Lagrange method. Finding the total number of real solutions to n simultaneous fractional algebraic equations is how the problem is formulated in [7]. Moulton’s work is significant in part because it solves the algebraic problem of determining the number of solutions to an impressive-looking set of simultaneous equations. In general, Moulton proved that there are $\frac{n!}{2}$ distinct central configurations for the collinear n-body problem, which is an invariant subsystem of the n-body problem.

In the current problem, we are not counting distinct homographic solutions (of which there are 12). We are looking at a symmetric case of this problem and giving a relationship between variational structure and the symmetric homographic solutions.

In [8], Gordon demonstrated that the elliptic Keplerian orbit minimizes the Lagrangian action of the two-body problems with periodic boundary conditions. The minimum value he obtained of the Keplerian action functional was computed as

$$A^K={\int }_0^T\frac{|\dot{q}{|}^2}{2}+\frac{\mu }{\left|q\right|}\ge \frac{3}{2}{\left(\pi \right)}^{2/3}{\left(\mu \right)}^{2/3}{T}^{1/3}.$$

(1)

The authors in [9,10] demonstrated that the Lagrangian and Eulerian elliptical solutions for the planar three-body problem are the variational minimizers of the Lagrangian action functional. Chen [11] discovered a new periodic solution family for the planar four-body problem with equal masses, minimizing solutions over a quarter of the period $[0,T]$ using numerical integration. Furthermore, the authors extended Chen’s solutions to include minimization over the entire period $[0,4T]$ without relying on numerical techniques.

This paper advances the understanding of the collinear four-body problem by demonstrating that homographic solutions for this symmetric configuration minimize the action function without the need for numerical integration. By analyzing the conditions for action minimization in Section 3 and employing numerical experiments with Poincaré sections, we provide a comprehensive exploration of periodic and quasi-periodic solutions. These findings enhance the theoretical framework for studying the stability and instability properties of homographic solutions in the symmetric collinear four-body problem, offering significant implications for both theoretical research and practical applications in celestial mechanics.

2. Symmetric Collinear Equilibrium Configurations

In this section, we are following the discussion in [12] for collinear configuration of a symmetric four-body problem.

Define the two mass ratios as

$${\mu }_1=\frac{m}{M_t},{\mu }_2=\frac{M}{M_t},$$

(2)

where m and M are the two unequal masses, and $M_t$ is the total mass of the system. For the symmetrical arrangement of the two pairs of masses, there are two possible cases. In one case, the pair with masses m lies in the middle of the collinear arrangement, and in the second case, their positions are swapped. Both cases are symmetric about the center of mass C.

Let

$$q_2=k{q}_1,k\in (0,1),q_3=-q_2,q_4=-q_1,M_t=2(m+M).$$

The equations of motion for $q_1$ and $q_2$ are

$${\ddot{q}}_1=-\frac{M_t}{q_1^3}{R}_1{q}_1,{\ddot{q}}_2=-\frac{M_t}{q_1^3}{R}_2{q}_2,$$

where

$$\begin{array}{ccc}\hfill R_1& =& {\mu }_2\left(\frac{1}{{(1-k)}^2}+\frac{1}{{(1+k)}^2}\right)+\frac{{\mu }_1}{4},\hfill \\ \hfill R_2& =& \frac{{\mu }_1}{k}\left(\frac{1}{{(1+k)}^2}-\frac{1}{{(1-k)}^2}\right)+\frac{{\mu }_2}{4k^3}.\hfill \end{array}$$

For the equilibrium solutions to exist, we must have $R_1=R_2$. Therefore,

$$A{\mu }_1-B{\mu }_2=0,$$

(3)

where

$$A=\frac{1}{{(1-k)}^2}-\frac{1}{{(1+k)}^2}+\frac{k}{4},B=\frac{1}{4k^2}-\frac{k}{{(1-k)}^2}-\frac{k}{{(1+k)}^2}.$$

We know from Equation (2) that ${\mu }_1+{\mu }_2=1/2$. Therefore,

$${\mu }_1=\frac{B}{2(A+B)},{\mu }_2=\frac{A}{2(A+B)}.$$

Lemma 1. 

For $k\in (0,1),$ $A+B>0,$ $A>0$, and $B>0$ when $k\in (0,0.41)$

Proof. 

It is trivial to show that $A>0$ and $A+B>0$ when $k\in (0,1)$. To show that $B>0$ for $k\in (0,0.41)$, we need to show that

$$\frac{1}{4k^2}>k\left(\frac{1}{{(1-k)}^2}+\frac{1}{{(1+k)}^2}\right),$$

or

$$\frac{1}{8k^3}>\frac{1+k^2}{{(1-k^2)}^2}.$$

Equivalently,

$$\phi \left(k\right)=8k^5-k^4+8k^3+2k^2-1<0.$$

It is easy to show that ${\phi }^{\prime }\left(k\right)>0$ for k$\in (0,1).$ Therefore, $\phi \left(k\right)$ is an increasing function for k$\in (0,1).$ Also, $\phi \left(0\right)<0$ and $\phi \left(1\right)>0;$ therefore, there exists only one root, $k_0\cong 0.41$. For $k<k_0,$$\phi \left(k\right)<\phi \left(k_0\right)=0$. Hence, $B>0$ for $0<k<k_0.$ The proof is also confirmed numerically, as can be seen in Figure 1. The first graph doesn’t clearly show the function $B\left(k\right)$ because its values are significantly larger than those of $A\left(k\right)$, especially for small values of k. This makes $A\left(k\right)$ appear indistinguishable from $A\left(k\right)+B\left(k\right)$. To address this, we have zoomed in on a smaller range of k so that the behavior of $A\left(k\right)$ and $B\left(k\right)$ is more visible. Note that a logarithmic scale is used on the vertical axis to better distinguish the functions. □

3. Action Minimizing Orbits

This section introduces the analytical discussion of a family of periodic solutions in the collinear four-body problem using variational techniques. For this problem, we let $m_1,m_2,m_3,m_4$ be four point masses moving in ${\mathbb{R}}^2$ in accordance with Newton’s Second Law—that is to say,

$$m_i\ddot{q_i}=\frac{\partial \mathcal{U}}{\partial q_i},$$

where $q_i\in {\mathbb{R}}^2$ denotes the position and $m_i$ the mass of the ith particle [13,14]. With configuration $q=(q_1,q_2,q_3,q_4)$, the force function $\mathcal{U}\left(q\right)$ (negative of the potential energy) is defined as

$$\mathcal{U}\left(q\right)=\sum _{i<j}\frac{m_i{m}_j}{|q_i-q_j|}.$$

The kinetic energy is defined as

$$\mathcal{K}=\sum _{i=1}^4\frac{m_i}{2}{\left|\dot{q_i}\right|}^2,$$

while the Hamiltonian governing the equations of motion is

$$\mathcal{H}(q,p)=\mathcal{K}\left(p\right)-\mathcal{U}\left(q\right),p=\frac{\partial \mathcal{L}}{\partial \dot{q}},$$

where the Lagrangian $\mathcal{L}$ is defined as

$$\mathcal{L}(q,\dot{q})=\sum _{i=1}^4\frac{m_i{\left|\dot{q_i}\right|}^2}{2}+\sum _{i<j}\frac{m_i{m}_j}{|q_i-q_j|},\dot{q_i}=\frac{dq}{dt}.$$

We now define the action functional to be of the form

$$\mathcal{A}\left(q\right)={\int }_0^T\mathcal{L}(q,\dot{q})dt,$$

where the Lagrangian $\mathcal{L}$ is the difference between the kinetic energy and the potential energy, as described above.

We focus on solutions called the homographic solutions of the form

$$q_i\left(t\right)=\varphi \left(t\right)q_{i,0},$$

(4)

where $q_{i,0}\in {\mathbb{R}}^2$, $\varphi \left(t\right):[0,T]\to {\mathbb{C}}^{\ast }$ is smooth, $\varphi \left(0\right)=\varphi \left(T\right)$, and $\mathrm{deg}\left(\varphi \right)\ne 0$. The following theorem proves that these solutions are the variational minimizers of the action function, and the minimum value is also computed.

Theorem 1. 

The Keplerian elliptical solutions are the minimal regular solutions to the symmetric collinear four-body problem with two pairs of equal masses, with the minimum action equal to $\frac{3}{2}{\left(2\pi \right)}^{2/3}{T}^{1/3}{\left(\frac{\xi \left(k\right)}{\eta \left(k\right)}\right)}^{2/3}$.

Proof. 

We aim to restrict the action functional to solutions defined by Equation (4). Let $\dot{q_1}\left(t\right)=\dot{\varphi }\left(t\right)q_{1,0}$ and the total mass $M_t=1$. Observe that $r_1=\left|q_{1,0}\right|\in \mathbb{R}$; then, $|{\dot{q}}_1{\left(t\right)|}^2={\left|\dot{\varphi }\left(t\right)\right|}^2{r}_1^2$.

Using the collinear constraints defined in Section 2, the kinetic energy $\mathcal{K}$ can be calculated as

$$\begin{array}{ccc}\hfill \mathcal{K}& =& \frac{1}{2}\sum _{i=1}^4{m}_i{\left|{\dot{q}}_i\left(t\right)\right|}^2\hfill \\ & =& {\mu }_1|\dot{\varphi }{\left(t\right)|}^2{r}_1^2+{\mu }_2{k}^2{\left|\dot{\varphi }\left(t\right)\right|}^2{r}_1^2\hfill \\ & =& |\dot{\varphi }{\left(t\right)|}^2{r}_1^2\left({\mu }_1+{\mu }_2{k}^2\right)\hfill \\ & =& |{\dot{q}}_1{\left(t\right)|}^2\left({\mu }_1+{\mu }_2{k}^2\right).\hfill \end{array}$$

The potential energy $\mathcal{U}$ for the standard four-body problem is given by

$$\mathcal{U}=\sum _{1\le i<j\le 4}\frac{m_i{m}_j}{|q_i-q_j|}.$$

Using $|q_i-q_j|=|\varphi \left(t\right)\left|\right|q_{i,0}-q_{j,0}|$, and the collinear constraints, we obtain

$$\mathcal{U}=\frac{1}{\left|\varphi \right(t\left)\right|}\left(2{\mu }_1{\mu }_2\left(\frac{1}{r_1(1-k)}+\frac{1}{r_1(1+k)}\right)+\frac{1}{2r_1}\left({\mu }_1^2+\frac{{\mu }_2^2}{k}\right)\right).$$

Multiplying and dividing by $r_1$, the potential energy is further simplified as

$$\mathcal{U}=\frac{1}{|q_1|}\left(2{\mu }_1{\mu }_2\left(\frac{2}{(1-k)(1+k)}\right)+\frac{1}{2}\left({\mu }_1^2+\frac{{\mu }_2^2}{k}\right)\right).$$

The action restricted to this class of homographic solutions can be computed as below:

$$\begin{array}{ccc}\hfill \mathcal{A}& =& {\int }_0^T\left({\mu }_1+{\mu }_2{k}^2\right){\left|\dot{q_1}\right|}^{2}dt+{\int }_0^T\left[2{\mu }_1{\mu }_2\left(\frac{2}{(1-k)(1+k)}\right)+\frac{1}{2}\left({\mu }_1^2+\frac{{\mu }_2^2}{k}\right)\right]\frac{1}{|q_1|}dt\hfill \\ & =& 2\left({\mu }_1+{\mu }_2{k}^2\right){\int }_0^T\frac{|\dot{q_1}{|}^2}{2}dt+\left[{\mu }_1{\mu }_2\left(\frac{4}{(1-k)(1+k)}\right)+\frac{1}{2}\left({\mu }_1^2+\frac{{\mu }_2^2}{k}\right)\right]{\int }_0^T\frac{1}{|q_1|}dt.\hfill \end{array}$$

For simplification, let

$$\eta \left(k\right)=2\left({\mu }_1+{\mu }_2{k}^2\right),$$

(5)

and

$$\xi \left(k\right)={\mu }_1{\mu }_2\left(\frac{4}{(1-k)(1+k)}\right)+\frac{1}{2}\left({\mu }_1^2+\frac{{\mu }_2^2}{k}\right).$$

(6)

Then,

$$\mathcal{A}\left(q\right)=\eta {\int }_0^T\frac{|\dot{q_1}{|}^2}{2}dt+\xi {\int }_0^T\frac{1}{|q_1|}dt.$$

The infimum of $\mathcal{A}\left(q\right)$ is

$$\begin{array}{ccc}\hfill \underset{q}{inf}\mathcal{A}\left(q\right)& =& \underset{k>0}{inf}\underset{q_1}{inf}\left\{\eta {\int }_0^T\frac{|\dot{q_1}{|}^2}{2}dt+\xi {\int }_0^T\frac{1}{|q_1|}dt\right\}\hfill \\ & =& \underset{k>0}{inf}\left\{\eta \underset{q_1}{inf}\left({\int }_0^T\frac{|\dot{q_1}{|}^2}{2}dt+\frac{\xi }{\eta }{\int }_0^T\frac{1}{|q_1|}dt\right)\right\}.\hfill \end{array}$$

We use Gordon’s result [8] to calculate the following infimum

$$\underset{q_1}{inf}\left({\int }_0^T\frac{|\dot{q_1}{|}^2}{2}dt+\frac{\xi }{\eta }{\int }_0^T\frac{1}{|q_1|}dt\right)=\frac{3}{2}{\left(2\pi \right)}^{\frac{2}{3}}{T}^{1/3}{\left(\frac{\eta }{\xi }\right)}^{2/3}.$$

Then,

$$\begin{array}{ccc}\hfill \underset{q}{inf}\mathcal{A}\left(q\right)& =& \underset{k>0}{inf}\left\{\eta \left(k\right)\frac{3}{2}{\left(2\pi \right)}^{2/3}{T}^{1/3}{\left(\frac{\xi \left(k\right)}{\eta \left(k\right)}\right)}^{2/3}\right\}\hfill \\ & =& \underset{k>0}{inf}\left\{\frac{3}{2}{\left(2\pi \right)}^{2/3}{T}^{1/3}{\eta }^{1/3}\left(k\right){\xi }^{2/3}\left(k\right)\right\}.\hfill \end{array}$$

Let

$$\varphi \left(k\right)=\frac{3}{2}{\left(2\pi \right)}^{2/3}{T}^{1/3}{\eta }^{1/3}\left(k\right){\xi }^{2/3}\left(k\right).$$

The function $\varphi \left(k\right)$ attains its infimum at $k_0$ if and only if $\eta \left(k\right){\xi }^2\left(k\right)$ attains its infimum at $k_0$. To show that $\varphi \left(k\right)$ is convex, we need to show that ${\varphi }^{″}\left(k\right)>0.$ For this purpose, we rewrite $\varphi \left(k\right)$ as $\varphi \left(k\right)=\eta \left(k\right){\xi }^2\left(k\right)$ and

$${\varphi }^{″}\left(k\right)={\xi }^2\left(k\right){\eta }^{″}\left(k\right)+4\xi \left(k\right){\eta }^{\prime }\left(k\right){\xi }^{\prime }\left(k\right)+2\eta \left(k\right)\xi \left(k\right){\xi }^{″}\left(k\right)+2\eta \left(k\right){\xi }^{\prime 2}\left(k\right),$$

where ${\eta }^{\prime }\left(k\right),$${\eta }^{″}\left(k\right)$, ${\xi }^{\prime }\left(k\right)$, and ${\xi }^{″}\left(k\right)$ are given in Appendix A. After some subtle simplifications, we obtain ${\varphi }^{″}\left(k\right)$ as below:

$${\varphi }^{″}\left(k\right)=\frac{1}{32}{{(k-1)}^{-8}{(k+1)}^8{\left(k^5-9k^3+k^2-1\right)}^{-7}k\left(k^7-k^5+8k^3+k^2-1\right)P}_{39}\left(k\right),$$

where the polynomial expression $P_{39}\left(k\right)$ was calculated using Mathematica and is expressed below:

$$\begin{array}{ccc}\hfill P_{39}\left(k\right)& =& 120k^{39}+1506k^{38}+1840k^{37}+64665k^{36}+307985k^{35}+96864k^{34}\hfill \\ & & +5491999k^{33}-2687201k^{32}+8272517k^{31}+5933252k^{30}-51540909k^{29}\hfill \\ & & -325466k^{28}-249830675k^{27}+84856200k^{26}-202946269k^{25}-318531830k^{24}\hfill \\ & & +333244705k^{23}-330109312k^{22}+9930530071k^{21}+1592903248k^{20}\hfill \\ & & +11885946259k^{19}-4700048448k^{18}+186685197k^{17}+139086048k^{16}\hfill \\ & & +1764471199k^{15}+4018007916k^{14}-971579575k^{13}-637345046k^{12}\hfill \\ & & -418692385k^{11}+210607336k^{10}+328772337k^9-45096250k^8-2504317k^7\hfill \\ & & -17372690k^6+1890269k^5-45689k^4+57936k^3+4336k^2+561.\hfill \end{array}$$

(7)

We note that the factor ${(k-1)}^{-8}{(k+1)}^8{(k^5-9k^3+k^2-1)}^{-7}k$ is negative when $k\in (0,1)$. The factor $k^7-k^5+8k^3+k^2-1$ has only one positive real root at $k=0.464405$ and is negative when $k<0.464405$ and positive when $k>0.464405.$ The last factor $P_{39},$ which is a polynomial of degree 39, has only 2 real positive roots: $0.238166$ and $0.387316$. The polynomial $P_{39}$ is positive when $k\in (0,0.238)\cup (0.387,1)$ and is negative when $k\in (0.238,0.387),$ see Figure 2. Therefore, ${\varphi }^{″}\left(k\right)$ is positive when $k\in (0,0.238)\cup (0.387,0.4644)$, see Figure 3. This proves that the function $\varphi \left(k\right)$ is convex when $k\in (0,0.238)\cup (0.387,0.4644)$.

For coercivity, we see that $\varphi \left(k\right)$ is continuous for all positive values of k, $\varphi \left(k\right)\to \infty$ as $k\to 0$ and when $k\to 1$, $\varphi \left(k\right)$ tends to ∞, which implies $\varphi \left(k\right)$ is coercive. Hence, $\varphi \left(k\right)$ attains ${inf}_{\left(k\right)>0}\left\{\varphi \left(k\right)\right\}$ at unique positive $\left(k_0\right)$ and satisfies $\varphi \left(k_0\right)=0$. □

4. Numerical Examples

Recent investigations show that the four-star system, particularly two binaries (i.e., two pairs of twin stars that revolve around each other at great distances), is more common in the universe than previously thought [5,15,16,17].

The four-body problem was investigated using the geometrical restriction methods. The stability and dynamical evolution of symmetric quadruple stellar systems, exoplanetary systems of two planets orbiting a binary pair of stars, were analyzed [4,18,19,20,21,22,23,24].

Therefore, a special type of four-body problem with analytical and numerical investigation can contribute to the understanding of the dynamical behavior of quadruple systems (e.g., HD 98800 quadruple system with two pairs of stars orbiting each other [25,26]. For example, HD 98800 is a young and nearby quadruple system composed of two spectroscopic binaries orbiting around each other (AaAb and BaBb), with a gas-rich disk in polar configuration around BaBb. The orbital parameters of BaBb and AB are relatively well constrained [26].

Providing an application for a real planetary system has been a concern of researchers for a long time. Investigations of quadruple stellar systems show that it is not easy to calculate a real star system because each system differs from mathematical models (in our case, a symmetric collinear four-body problem).

By studying the orbits of stars, we can obtain information about the formation processes and evolution of multiple star systems. The orbits of stars preserve information about these formation processes. Thus, by analytically and numerically investigating the collinear four-body problems, we can better understand the dynamic behavior of such quadruple-star systems.

We focus on the symmetric collinear four-body (SC4BP) central configurations presented in Section 2.

Let us consider the case of four point masses, where $m_1=m_4=m$, $m_2=m_3=M$ lie collinear on the x-axis and symmetrical to the y-axis at the origin, as introduced in Section 2 (see Figure 4). Using the position coordinates from Section 2, we obtain the following reduced Hamiltonian of the symmetric collinear four-body problem:

$$H(q_1,q_2,p_1,p_2)=\frac{p_1^2}{4m}+\frac{p_2^2}{4M}-\frac{m^2}{2q_1}-\frac{M^2}{2q_2}-\frac{2Mm}{q_2-q_1}-\frac{2Mm}{q_2+q_1},$$

(8)

where $\mathbf{q}=(q_1,q_2)\in {\mathbb{R}}_{>0}^2$ are generalized coordinates and $\mathbf{p}=(p_1,p_2)\in {\mathbb{R}}^2$ are generalized momenta.

Explicitly, the equations of motion are the following:

$$\begin{array}{ccc}\hfill {\dot{q}}_1& =& \frac{p_1}{2m},\hfill \\ \hfill {\dot{p}}_1& =& -\frac{m^2}{2q_1^2}+\frac{2Mm}{{(q_2-q_1)}^2}-\frac{2Mm}{{(q_2+q_1)}^2},\hfill \\ \hfill {\dot{q}}_2& =& \frac{p_2}{2M},\hfill \\ \hfill {\dot{p}}_2& =& -\frac{M^2}{2q_2^2}-\frac{2Mm}{{(q_2-q_1)}^2}-\frac{2Mm}{{(q_2+q_1)}^2}.\hfill \end{array}$$

(9)

We aim to find stable regions in special cases and analyze their stability. The dynamical parameters of our system are studied with phase portraits and Poincaré maps, which contain all the information necessary to characterize the dynamics of a system.

For the investigation of the reduced Hamiltonian equations of motion, we have tested two cases with different masses: $m<M$ and $m>M$. In these situations, we studied the quasi-periodic orbits implemented using the Poincaré surface of section technique by selecting the phase element $p_1=0$ and $q_1>0$.

The fourth-order Runge–Kutta method is used to compute and plot projections of the Poincaré surface of the Hamiltonian system (8) with optional step-size ($0.05$), number of iterations (5), and a range of time ($t\in [-500,500]$) [27]. Implementing a symbolic computing environment (in our case, Maple, Matlab, and Mathematica) allows us to gain deeper insight into the given problem (Equation (9)).

For the case $m>M$, with tentative different energy levels and initial conditions, we did not obtain any good results for Poincaré sections.

The masses m and M are chosen in the following way: m is arbitrarily set to $0.1$; for increasing values of h—$-0.077$, $9.8$, and 15—M is set to $0.8$, 6, and 7, respectively.

We explored other settings of h, M, and m, but only found illustrative examples in the cases mentioned above. Our numerical experiments focus on specific mass ratios and energy levels. The wider exploration of the parameter space for different mass and energy levels was based on empirical analysis. The obtained results seem very sensitive to initial conditions. We will investigate this matter further to ensure we reach the correct conclusion; in addition, it will be useful to predict the behavior of bodies passing close to collision and in which cases a strong sensitivity to initial conditions appears [28].

A color-coded version of the map is utilized to better visualize the chaotic behavior over small energy level variations.

We discovered some intriguing surfaces of sections in Figure 5 and Figure 6 for the condition $m<M$ and for energy levels $-0.077$, $9.8$, and 15. The three different regions are identified as follows:

• The inner part is the region around $p_2=0$.
• The middle part can be described as follows:

  -

  For $h=-0.77$, $p_2\in [-20,-3]\cup [3,20]$.

  -

  For $h=9.8$, $p_2\in [-700,-200]\cup [200,700]$.

  -

  For $h=15$, $p_2\in [-30,-5]\cup [5,30]$.

• The outer part is dominated by empty space with some stray points.

Note that at energy level $h=-0.077$, the inner part is empty, the outer part is very close to instability, and the middle part of the figure contains some invariant curves, but no stable region.

At $h=9.8$, the inner part is also empty space, the middle part again contains some invariant curves, and the outer part is chaotic, without any stable behavior.

In Figure 6, at $h=15$, it seems that the middle part also surrounds the inner part with small invariant curves while the outer part remains chaotic.

Furthermore, as the distance from the center of the surface to the outside increases, the number of scattered points decreases, suggesting instability.

The invariant curves are just small accumulated points that do not form closed loops, as expected in the case of KAM torus.

All three surfaces of the sections show instability, but the small invariant curves lead us to further research.

Consequently, for proper study of the motion of SC4BP, we need to regularize the equations.

The regularization in this case is necessary because the equations of motion (Equation (9)) show singularities. The continuation of the orbit after close encounters is not feasible since the solution encounters the singularity present in the problem.

In this paper, we focus on the simple binary collisions, which appear in the symmetric collinear four-body problem with collision singularity. In this case, at collisions between two bodies, the corresponding singularities can be removed from the differential equations using classical regularizing transformations (e.g., Levi-Civita or Kustaanheimo–Stiefel regularizations along the orbits or McGehee regularization considering the motion near collision, block regularization) [19,29,30,31,32].

McGehee’s local choice regularization (removing one of the two singularities) is the easiest one from the point of view of the physical meaning of the regularizing variables considered (which are essentially polar coordinates). The system of the symmetric collinear four-body problem has a simple expression, and one can analyze the collision manifold. The ejection or collision orbits in the McGehee variables become asymptotic solutions to equilibrium points. Levi-Civita local regularization results as more suitable from a numerical point of view, since the pass through collision is a regular point, but the expression of the resulting equations of motion is more intricate and appears to be a double covering of the phase space [33].

Moreover, ref. [22] proved that there is an important first integral C in the solution of the collinear four-body problem that ends in simultaneous binary collisions, and this constant C reveals the connection between the two collisions of simultaneous binary collision. Following [14], they use a Levi-Civita canonical transformation to regularize the equations of motion. We applied the McGehee regularization [34], which blows up the collision singularity and regularizes the equations of motion.

We make canonical changes to variables and time. Let us introduce the polar coordinates:

$$\begin{array}{ccc}\hfill q_1& =& rcos\theta ,\hfill \\ \hfill q_2& =& rsin\theta ,\hfill \\ \hfill p_1& =& p_{r}cos\theta -\frac{p_{\theta }}{r}sin\theta ,\hfill \\ \hfill p_2& =& p_{r}sin\theta +\frac{p_{\theta }}{r}cos\theta ,\hfill \end{array}$$

(10)

where $p_r=\frac{dr}{dt}=p_{1}cos\theta +p_{2}sin\theta$, $p_{\theta }=\frac{d\theta }{dt}=-p_{1}rsin\theta +p_{2}rcos\theta$, $r=\sqrt{q_1^2+q_2^2}$, and $\theta =arctan\frac{q_2}{q_1}$. For the case $M=m$, the Hamiltonian in Equation (8) becomes

$$\begin{array}{ccc}\hfill H(r,\theta ,p_r,p_{\theta })& =& \frac{p_r^2}{4m}+\frac{p_{\theta }^2}{4m{r}^2}-\frac{V\left(\theta \right)}{r}\hfill \\ \hfill V\left(\theta \right)& =& \frac{m^2}{2cos\theta }+\frac{m^2}{2sin\theta }+\frac{2m^2}{(sin\theta -cos\theta )}+\frac{2m^2}{(sin\theta +cos\theta )},\hfill \end{array}$$

(11)

and we obtain the following equations of motion:

$$\begin{array}{ccc}\hfill \dot{r}& =& \frac{p_r}{2m},\hfill \\ \hfill \dot{\theta }& =& \frac{p_{\theta }}{2m{r}^2},\hfill \\ \hfill {\dot{p}}_r& =& \frac{p_{\theta }^2}{2m{r}^3}-\frac{m^2}{r^2}·\frac{(sin\theta +cos\theta )}{sin2\theta }-\frac{4m^2}{r^2}·\frac{sin\theta }{{sin}^2\theta -{cos}^2\theta },\hfill \\ \hfill {\dot{p}}_{\theta }& =& \frac{m^2}{2r}·\left(\frac{sin\theta }{{cos}^2\theta }-\frac{cos\theta }{{sin}^2\theta }\right)-\frac{4m^2}{r}·\frac{sin2\theta }{{({sin}^2\theta -{cos}^2\theta )}^2}.\hfill \end{array}$$

(12)

We remark that the total collision is at the manifold $r=0$, but we still have discontinuities when $\theta =k\pi$ or $\theta =(2k+1)\frac{\pi }{2}$, where $k\in \mathbb{N}$. Firstly, we remove the singularity ($r=0$), introducing the following variables: $p_r=\frac{v}{\sqrt{r}},$ $p_{\theta }=\sqrt{r}u$.

The system (12) becomes

$$\begin{array}{ccc}\hfill \dot{r}& =& \frac{v{r}^{-\frac{1}{2}}}{2m},\hfill \\ \hfill \dot{\theta }& =& \frac{u{r}^{-\frac{3}{2}}}{2m},\hfill \\ \hfill \dot{v}& =& \frac{(u^2+\frac{v^2}{2})r^{-\frac{3}{2}}}{2m}-m^2{r}^{-\frac{3}{2}}·\frac{(sin\theta +cos\theta )}{sin2\theta }-4m^2{r}^{-\frac{3}{2}}·\frac{sin\theta }{{sin}^2\theta -{cos}^2\theta },\hfill \\ \hfill \dot{u}& =& -\frac{uv}{4m}{r}^{-\frac{3}{2}}+\frac{m^2}{2}{r}^{-\frac{3}{2}}·\left(\frac{sin\theta }{{cos}^2\theta }-\frac{cos\theta }{{sin}^2\theta }\right)-4m^2{r}^{-\frac{3}{2}}·\frac{sin2\theta }{{({sin}^2\theta -{cos}^2\theta )}^2}.\hfill \end{array}$$

(13)

We mention that $r·H$ is the first integral of the system in these new variables.

Consequently, we specify the following time variables $\frac{dt}{ds}=r^{\frac{3}{2}}$, denoted with ′ $=\frac{d}{ds}$, and the system (13) becomes

$$\begin{array}{ccc}\hfill r^{\prime }& =& \frac{vr}{2m},\hfill \\ \hfill {\theta }^{\prime }& =& \frac{u}{2m},\hfill \\ \hfill v^{\prime }& =& \frac{(u^2+\frac{v^2}{2})}{2m}-m^2·\frac{(sin\theta +cos\theta )}{sin2\theta }-4m^2·\frac{sin\theta }{{sin}^2\theta -{cos}^2\theta },\hfill \\ \hfill u^{\prime }& =& -\frac{uv}{4m}+\frac{m^2}{2}·\left(\frac{sin\theta }{{cos}^2\theta }-\frac{cos\theta }{{sin}^2\theta }\right)-4m^2·\frac{sin2\theta }{{({sin}^2\theta -{cos}^2\theta )}^2}.\hfill \end{array}$$

(14)

Moreover, to regularize the discontinuities, we introduce the following variable:

$$w\left(\theta \right)=\sqrt{\frac{cos\theta (sin\theta -cos\theta )}{V\left(\theta \right)}}·u\left(\theta \right).$$

(15)

We note that w is regular in the Hill region, and the system becomes

$$\begin{array}{ccc}\hfill r^{\prime }& =& \frac{vr}{2m},\hfill \\ \hfill {\theta }^{\prime }& =& \frac{w}{2m}·\sqrt{\frac{V\left(\theta \right)}{cos\theta (sin\theta -cos\theta )}},\hfill \\ \hfill v^{\prime }& =& 2rH-\frac{v^2}{4m}+V\left(\theta \right),\hfill \\ \hfill w^{\prime }& =& -\frac{vw}{4m}·\sqrt{\frac{V\left(\theta \right)}{cos\theta (sin\theta -cos\theta )}}+\frac{m^2}{2}·\left(\frac{sin\theta }{{cos}^2\theta }-\frac{cos\theta }{{sin}^2\theta }\right)-\hfill \\ & & -4m^2·\frac{sin2\theta }{{({sin}^2\theta -{cos}^2\theta )}^2}.\hfill \end{array}$$

(16)

Now, we introduce the time variable s as the second time transformation:

$$\frac{dt}{ds}=\sqrt{\frac{cos\theta (sin\theta -cos\theta )}{V\left(\theta \right)}},$$

(17)

denoted with * $=\frac{d}{ds}$, and we obtain the following new system:

$$\begin{array}{ccc}\hfill r^{\ast }& =& \frac{vr}{2m}·\sqrt{\frac{cos\theta (sin\theta -cos\theta )}{V\left(\theta \right)}},\hfill \\ \hfill {\theta }^{\ast }& =& \frac{w}{2m},\hfill \\ \hfill v^{\ast }& =& [2rH-\frac{v^2}{4m}+V\left(\theta \right)]·\sqrt{\frac{cos\theta (sin\theta -cos\theta )}{V\left(\theta \right)}},\hfill \\ \hfill w^{\ast }& =& -\frac{vw}{4m}+\left[\frac{m^2}{2}·\left(\frac{sin\theta }{{cos}^2\theta }-\frac{cos\theta }{{sin}^2\theta }\right)-4m^2·\frac{cos\theta }{{({sin}^2\theta -{cos}^2\theta )}^2}\right]·\hfill \\ & & ·\sqrt{\frac{cos\theta (sin\theta -cos\theta )}{V\left(\theta \right)}}.\hfill \end{array}$$

(18)

In the new variable $(r^{\ast },{\theta }^{\ast },v^{\ast },w^{\ast })$, we can write the following part using Equation (15) and the first integral of the system:

$$\frac{w^{\ast 2}}{2cos{\theta }^{\ast }(sin{\theta }^{\ast }-cos{\theta }^{\ast })}+\frac{v^{\ast 2}}{2V\left({\theta }^{\ast }\right)}=1+\frac{rH}{V\left({\theta }^{\ast }\right)},$$

(19)

and to study the periodic orbits of Equation (18) for a fixed energy level H and mass m, we introduce the solution of $\frac{dV}{d\theta }=0$, namely, the critical value ${\theta }_{crit}$.

For negative energy, the left part of relation (19) becomes less than or equal to 1, which depicts in $(w,v)$ variables the intersection between the total collision and the section $\theta ={\theta }_{crit}$.

For this critical section, where we took $h=-0.077$ and $m=1$, we plot the Poincaré surface of sections (Figure 7), which gives an inside look at how the given system acts.

We can observe a quasi-periodic region as a vent in the center of Figure 7, where there are twisted lines on a chaotic backdrop.

The obtained results appear very sensitive to initial conditions. To the best of our knowledge, a traditional way to detect chaotic behavior of orbits in dynamical systems is the calculation of the maximum Lyapunov characteristic exponent (LCE). The development of fast and efficient methods of detecting chaotic behavior in gravitational dynamical systems, such as the method of the relative Lyapunov indicator (RLI), can enable detailed pictures of the chaotic behavior of the studied system [35]. In the future, we will analyze the four-body problem from this point of view.

These results show how much the stability of the collinear four-body problem depends on the initial conditions. The quasi-periodic region could denote a stable manifold, and the outer part is an unstable manifold with unstable orbits.

5. Conclusions and Future Work

In this study, we investigated periodic and quasi-periodic orbits in the symmetric collinear four-body problem using variational methods. The symmetry is used here for the purposes of allowing simplified numerical calculations. Our primary findings show that homographic solutions with equal masses minimize the action functional. We also computed the minimum value of the action specifically for these solutions. Additionally, through numerical experiments using Poincaré sections, we explored the existence and stability of periodic and quasi-periodic solutions within the broader dynamics of the symmetric collinear four-body problem.

These results have significant implications for understanding celestial mechanics and multi-body system stability. By confirming that homographic solutions minimize the action, we lay the groundwork for future stability analyses in the collinear four-body problem. The numerical techniques employed, particularly Poincaré sections, offer valuable insights into complex dynamical behaviors. Studying the stability and its relation to chaotic and quasi-periodic regions of motion in four-body systems is, therefore, fundamental to understanding the evolution of quadruple stellar and exoplanetary systems.

Future research could extend these methods to non-symmetric configurations and explore the effects of mass variations among the bodies. Investigating the application of these techniques to three-dimensional problems could further enhance our understanding of celestial dynamics. Advanced numerical methods and computational power could refine solutions and contribute to astrodynamics and space mission planning. In the case of non-symmetric configurations, we can work with collinear configurations $q=(q_1,q_2,q_3,q_4)\in {\mathbb{R}}^4$, and the action functional is determined from the four-body problem action functional.