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

1. Introduction

Let $m_1,m_2,m_3,m_4$ be four point masses moving in ${\mathcal{R}}^2$ in accordance with Newton’s Second Law, that is to say,

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

where $q_i\in {\mathcal{R}}^2$ denotes the position, and $m_i$ the mass of the ith particle [1,2,3,4]. 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}{\displaystyle \frac{m_i{m}_j}{|q_i-q_j|}}.$$

The Kinetic energy is defined as

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

while the Hamiltonian governing the equations of motion is

$$\mathcal{H}(q,\dot{q})=\mathcal{K}\left(\dot{q}\right)-\mathcal{U}\left(q\right).$$

For this problem, we 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 defined as

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

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

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

(1)

The authors in [6,7] have demonstrated that the Lagrangian and Eulerian elliptical solutions for the planar three-body problem are the variational minimizers of the Lagrangian action functional. In [8], it was shown that the homographic solutions to the rhombus four-body problem are the variational minimizers of the action functional restricted to a rhombus loop space. Additionally, the authors in [9,10] identified the minimizers of the classical action functional restricted to the homographic solutions of the kite and trapezoidal four-body problems as the Keplerian elliptical solutions. 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 in [12] 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 [16], for collinear configuration of a symmetric four-body problem.

Define the two mass ratios as

$${\mu }_1={\displaystyle \frac{m}{M_t}},{\mu }_2={\displaystyle \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 lie in the middle of the collinear arrangement and in the second case their position is 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=-{\displaystyle \frac{M_t}{q_1^3}}{R}_1{q}_1,{\ddot{q}}_2=-{\displaystyle \frac{M_t}{q_1^3}}{R}_2{q}_2,$$

where

$$\begin{array}{ccc}\hfill R_1& =& {\mu }_2\left({\displaystyle \frac{1}{{(1-k)}^2}}+{\displaystyle \frac{1}{{(1+k)}^2}}\right)+{\displaystyle \frac{{\mu }_1}{4}},\hfill \\ \hfill R_2& =& {\displaystyle \frac{{\mu }_1}{k}}\left({\displaystyle \frac{1}{{(1+k)}^2}}-{\displaystyle \frac{1}{{(1-k)}^2}}\right)+{\displaystyle \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={\displaystyle \frac{1}{{(1-k)}^2}}-{\displaystyle \frac{1}{{(1+k)}^2}}+{\displaystyle \frac{k}{4}},B={\displaystyle \frac{1}{4k^2}}-{\displaystyle \frac{k}{{(1-k)}^2}}-{\displaystyle \frac{k}{{(1+k)}^2}}.$$

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

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

Lemma: 

 

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

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

or

$${\displaystyle \frac{1}{8k^3}}>{\displaystyle \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.$ □

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. We focus on solutions called the homographic solutions of the form $q_i\left(t\right)=\varphi \left(t\right)q_{i,0}$ where $\varphi :[0,T]\to {\mathbb{C}}^{*}$ 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: 

 

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 ${\displaystyle \frac{3}{2}}{\left(2\pi \right)}^{2/3}{T}^{1/3}{\left({\displaystyle \frac{\xi \left(k\right)}{\eta \left(k\right)}}\right)}^{2/3}$

Proof. 

We aim to restrict the action functional to homographic solutions of the above-mentioned form.

Let $|{\dot{q}}_1{\left(t\right)|}^2={\left|\dot{\varphi }\left(t\right)\right|}^2{r}_1^2$, and the total mass $M_t=1$, then the kinetic energy $\mathcal{K}$ is equal to

$$\begin{array}{ccc}\hfill \mathcal{K}& =& {\displaystyle \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 is given by

$$\mathcal{U}=\sum _{1\le i<j\le 4}{\displaystyle \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}|$, we get

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

Multiplying and dividing by $r_1$ we obtain

$$\mathcal{U}={\displaystyle \frac{1}{|q_1|}}\left(2{\mu }_1{\mu }_2\left({\displaystyle \frac{2}{(1-k)(1+k)}}\right)+{\displaystyle \frac{1}{2}}\left({\mu }_1^2+{\displaystyle \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({\displaystyle \frac{2}{(1-k)(1+k)}}\right)+{\displaystyle \frac{1}{2}}\left({\mu }_1^2+{\displaystyle \frac{{\mu }_2^2}{k}}\right)\right]{\displaystyle \frac{1}{|q_1|}}dt\hfill \\ & =& 2\left({\mu }_1+{\mu }_2{k}^2\right){\int }_0^T{\displaystyle \frac{|\dot{q_1}{|}^2}{2}}dt+\left[{\mu }_1{\mu }_2\left({\displaystyle \frac{4}{(1-k)(1+k)}}\right)+{\displaystyle \frac{1}{2}}\left({\mu }_1^2+{\displaystyle \frac{{\mu }_2^2}{k}}\right)\right]{\int }_0^T{\displaystyle \frac{1}{|q_1|}}dt.\hfill \end{array}$$

Let

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

(4)

and

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

(5)

Then

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

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

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

Then

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

Let

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

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

$${\varphi }^{\prime \prime }\left(k\right)=\xi {\left(k\right)}^2{\eta }^{\prime \prime }\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 }^{\prime \prime }\left(k\right)+2\eta \left(k\right){\xi }^{\prime 2}\left(k\right),$$

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

$${\varphi }^{\prime \prime }\left(k\right)={\displaystyle \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

$$\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}$$

(6)

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 (1). Therefore ${\varphi }^{\prime \prime }\left(k\right)$ is positive when $k\in (0,0.238)\cup (0.387,0.4644),$ see Figure (2).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 \infty$, $\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$. □

Figure 1. 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)$

Figure 1. 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)$

Figure 2. The function $\varphi \left(k\right)$ is convex when $k\in (0,0.238)\cup (0.387,0.4644)$

Figure 2. The function $\varphi \left(k\right)$ is convex when $k\in (0,0.238)\cup (0.387,0.4644)$

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. 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.

Figure 3. Configuration of the symmetric collinear four-body (SC4BP) problem.

Figure 3. Configuration of the symmetric collinear four-body (SC4BP) problem.

Let us consider the case of four point masses, where $m_1=m_4=m$, $m_2=m_3=M$ lying collinearly on the x-axis and symmetrical to the y-axis at the origin, as introduced in Section 2 (see Figure 3). 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)={\displaystyle \frac{p_1^2}{4m}}+{\displaystyle \frac{p_2^2}{4M}}-{\displaystyle \frac{m^2}{2q_1}}-{\displaystyle \frac{M^2}{2q_2}}-{\displaystyle \frac{2Mm}{q_2-q_1}}-{\displaystyle \frac{2Mm}{q_2+q_1}},$$

(7)

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& =& {\displaystyle \frac{p_1}{2m}},\hfill \\ \hfill {\dot{p}}_1& =& -{\displaystyle \frac{m^2}{2q_1^2}}+{\displaystyle \frac{2Mm}{{(q_2-q_1)}^2}}-{\displaystyle \frac{2Mm}{{(q_2+q_1)}^2}},\hfill \\ \hfill {\dot{q}}_2& =& {\displaystyle \frac{p_2}{2M}},\hfill \\ \hfill {\dot{p}}_2& =& -{\displaystyle \frac{M^2}{2q_2^2}}-{\displaystyle \frac{2Mm}{{(q_2-q_1)}^2}}-{\displaystyle \frac{2Mm}{{(q_2+q_1)}^2}}.\hfill \end{array}$$

(8)

We aim to find some stable regions in special cases and analyze their stability over time. 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 quasiperiodic orbits implemented using the Poincaré surface of section technique by selecting the phase element $p_1=0$ [13,14].

In Figure 4 and Figure 5, with $m<M$, and for energy levels $-0.07762,9.8,15$ we found some interesting surface of sections.

At the energy level $H=-0.07762$, 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$ in the inner part is also empty space, and the outer part is chaotic, without any stable behaviors.

In Figure 5, at $H=15$, it seems that the chaotic outer part surrounds the inner part with small invariant curves.

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

Consequently, for proper study of the motion of SC4BP, we need to regularize the equations [15]. We make canonical changes of 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 -{\displaystyle \frac{p_{\theta }}{r}}sin\theta ,\hfill \\ \hfill p_2& =& p_{r}sin\theta +{\displaystyle \frac{p_{\theta }}{r}}cos\theta ,\hfill \end{array}$$

(9)

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

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

(10)

and we obtain the following equations of motion:

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

(11)

We remark the total collision is at the manifold $r=0$, but still have discontinuities when $\theta =k\pi$ or $\theta =(2k+1){\displaystyle \frac{\pi }{2}}$, where $k\in \mathbf{N}$.

Firstly, we remove the singularity ($r=0$), introducing the following variables: $p_r={\displaystyle \frac{v}{\sqrt{r}}},$ $p_{\theta }=\sqrt{r}u$.

The system (11) becomes

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

(12)

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

Consequently, we specify the following time variables ${\displaystyle \frac{dt}{ds}}=r^{{\displaystyle \frac{3}{2}}}$, and we note with ${}^{\prime }={\displaystyle \frac{d}{ds}}$, and the system (12) becomes

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

(13)

Furthermore, to regularize the discontinuities we introduce the variable:

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

(14)

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

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

(15)

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

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

(16)

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

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

(17)

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

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

(18)

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

For negative energy the left part of relation (18) 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.07762$ and $m=1$, we plot the Poincaré surface of sections (Figure 6), which gives an inside look at how the given system acts.

We can observe a quasiperiodic region as a vent in the center of Figure 6, and there were twisted lines on a chaotic backdrop. These results show how much the stability of the collinear four-body problem depends on the initial conditions. The quasiperiodic region could denote a stable manifold, and the outer part is an unstable manifold with unstable orbits.

5. Conclusions

In this study, we have investigated the periodic and quasi-periodic orbits in the symmetric collinear four-body problem using variational methods. Our primary findings demonstrate that the homographic solutions for this problem with equal masses indeed minimize the action functional. We have also computed the minimum value of the action functional specifically for these solutions. Furthermore, through numerical experiments employing Poincaré sections, we have 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 the understanding of celestial mechanics and the stability of multi-body systems. By confirming that homographic solutions minimize the action, we provide a foundation for future stability analyses of such solutions in the collinear four-body problem. The numerical techniques used in this study, particularly the Poincaré sections, offer valuable insights into the complex dynamical behaviors and can be extended to study more general configurations and other types of symmetry in multi-body problems.

Future research could extend these methods to non-symmetric configurations and explore the implications of mass variations among the bodies. Additionally, investigating the application of these variational techniques to three-dimensional multi-body problems could provide a more comprehensive understanding of celestial dynamics. The use of advanced numerical methods and computational power can further refine the solutions and contribute to the broader field of astrodynamics and space mission planning.

In summary, our findings contribute to the theoretical understanding of action-minimizing orbits in the collinear four-body problem, opening avenues for further research in both theoretical and applied celestial mechanics.