## Abstract

Thousands of solutions have been found to the three-body problem; however, the stability of the discovered solutions is rarely discussed. This paper tests the stability of eight solutions to the three-body problem found by Suvakov and Dmitrasinovic in 2013. The solutions are tested for stability in the case of equal masses () and the addition a fourth body to the system of a limited mass (between 0.001 and 0.01). It was found that many solutions exhibit some kind of gravitational stability after the addition of a fourth body. The results showed that each solution had a unique range of masses at which it retained its gravitational binding. Other systems showed complete chaos after introducing even a small fourth body. Factors of stability such as the position of the fourth body in addition to its initial velocity were tested, adding to the precision of the tests. All four-body systems studied exhibited significantly chaotic configurations despite retaining relative gravitational binding in some cases. The results can help us further understand these complex systems and expand on to the n-body problem and the production of weak gravitational waves due to the chaos exhibited in such systems.

# Introduction

The three-body problem, referring to the gravitational interaction of three bodies in space, has been puzzling since the days of Newton. Written at the end of the volume of Opticks, Newton expressed doubts about the stability of a system of mutual gravitational interaction among more than two bodies, stating that their Keplerian motion will be disturbed^{1}. Different attempts have been made to solve the problem in the same way as that of two bodies, but no solutions were found, and the realization of this failure led to Poincareâ€™s statement that the system of equations of motion of the three-body problem is non-integrable^{2}. Despite the chaos that exists in the system, Euler and Lagrange were able to find instances at which the system exhibited stability over periods of time. The solutions found were called periodicals, and scientists have since been finding additional periodic solutions to the system.

The solutions to the three-body problem are categorized in families according to the topographical methods of Montgomery. The first discovered family was Lagrange-Euler family, dating back to the 18th century^{3}. The second discovered family was the The Broucke-Hadjidemetriou-Henon family (BHH family), discovered by Broucke (1975), Hadjidemetriou (1975), and Henon (1976)^{4}. The third family, the figure-8 family, was discovered by Moore in 1993 and rediscovered by Montgomery in the 2000s^{5}. A breakthrough occurred when Suvakov and DmitraÅ¡inovi? discovered 13 new families in 2013; hundreds of solutions have been discovered afterwards^{6}.

The three-body problem has applications in different fields like astronomy, physics, and chemistry. It was used to study the effect of the sun on the moon and the trojan system. Using the idea of the three-body problem, the literature was able to expand it to the n-body problem, a system collecting nâ€“bodies gravitationally, allowing the study of the solar system and clusters of stars containing millions of stars. In 2012, Laskar used the idea of the n-body problem to study the stability of the solar system^{1}. In chemistry, collision theory can be used to describe the reaction of three elements^{7}. In this paper, the stability of zero angular momentum systems will be scrutinized for stability after the introduction of a fourth body considering three factors: the mass of the fourth body, its initial position, and its initial velocity.

Methodology

**Notation and Conditions**

I will be using Newtonâ€™s notation for differentiation (i.e, ). In addition, a geometrized unit system is used for the sake of simplicity in computations (G = 1; , where G is the universal gravitational constant and is the mass of the body i in the plane). The masses, distances, and time in the system are all relative, unitless quantities.

The equations used in this paper are all described for two-dimensional space (the x-y plane). Since trio-systems always exist in a two-dimensional plane, then can be used to describe the system; this was proved in Montgomeryâ€™s paper about the topology of the three-body systems^{8}. The four-bodies case will, however, require working in the three-dimensional space (x-y-z). Dealing with the system in a 3-dimensional space is power-consuming, so to decrease dependency on computational power, the systems will all be tested in the two-dimensional plane; however, the methods developed in the paper can be used to test the stability in a three-dimensional notion.

## Methodology

**Notation and Conditions**

I will be using Newton’s notation for differentiation (i.e ). In addition, a geometrized unit system is used for the sake of simplicity in computations ( , where G is the universal gravitational constant and is the mass of the body in the plane). The masses, distances, and time in the system are all relative, unitless quantities.

The equations used in this paper are all described for two-dimensional space (the x-y plane). Since trio-systems always exist in a two-dimensional plane, then can be used to describe the system; this was proved in Montgomeryâ€™s paper about the topology of the three-body systems^{8}. The four-bodies case will, however, require working in the three-dimensional space (x-y-z). Dealing with the system in a 3-dimensional space is power-consuming, so to decrease dependency on computational power, the systems will all be tested in the two-dimensional plane; however, the methods developed in the paper can be used to test the stability in a three-dimensional notion.

**Equations of Motion**

Objects are attracted to each other gravitationally due to their masses according to the equations of motion proposed by Newton, so the gravitational acceleration of each body follows Equation (1):

where is the unit direction of the position vector between the specified body and every other body surrounding it, is the magnitude of the positional vector between the sought object and each surrounding object i. is the second derivative of the position vector, is the gravitational constant, and is the mass of each body. The differential equations of motion that describe, the interactions among bodies can be determined by positional vectors and passes as follows in. Equations (2) and(3):

(1)

The above equation describes the motion of body 1 . Similarly, we can describe the game two equations for both bodies 2 and 3 in the system ( and ), hence we have gix equations of motion to describe the system fully. When the system contains a fourth body (as in the case tested in the paper), we must include a fourth factor to both and equations of motions of the body. The terms are expressed in Equations (4) and (5).

To study the evolution of the system, we have to solve for all the objects’ positions ( ) and their derivative with respect to time ( .As a result, we are working with a 12 -vector denoted as follows in Equation (6):

Note that each component has two dimensions and , which makes the 12 -vector

The aforementioned equations can guarantee conservation in angular momentum, linear momentum, and mechanical energy of the system. The angular momentum can be calculated using Equation (7):

where and are the position and momentum vectors respectively. In addition, canbedefined as the mass multiplied the velocity of the body ( ). Since we are working in a twodimensional plane, the magnitude of the angular momentum will be as shown in Equation (8).

Where and are the and co-ordinates of the momentum of a body . In addition, as the system in hand is a geometries unit system , the first derivative of the position vector can be used instead of the linear momentum in the above equation to get the angular momentum.

The system also guarantees a conservation in the mechanical energy (the sum of the kinetic energy and the potential energy of the system). Equations (9) and (10) areused to calculate these two quantities.

The mechanical energy of the system will be the summation of the magnitude of both kinetic energy and potential energy as shown in Equation (11).

**Return Proximity Function**

All the solutions to the three-body problem that waill be discossed in this paper can be checked for periodicity using the retuon proximity flontion. The fanction is defined as the absolute minimum of the distance from the initial condition by as shown in Equation (12).

To find a solution with some period , the zeros of the function must be found at the chosen^{9}. The algorithm used to solve the return proximity function basically depends on the choosing some period Â and an initial condition Â and solving the equations of the motion to find Â from . After that, the distance between Â and the initial condition is measured using linear interpolation. Initially, the distance to the initial condition will increase from . When the distance starts decreasing, one must check whether it is minimal or not. If the distance is minimal and is smaller than some arbitrary tolerance (error), we consider this initial condition to be a candidate for a periodic solution^{10}. Every configuration of the three bodies (the shape of the triangle formed by them, independent of size), can be represented by a point in a unit sphere as the one shown in **fig. 1**, known as the shape sphere^{10}.

The north and the south poles of the shape sphere corresponds to equilateral triangles that differ only in orientation of the bodies. In fact, a similar property holds true for the whole northern and southern hemispheres; two triangles of the same shape but with opposite orientations are represented by two points that are symmetrical relative to the equatorial plane. The equator corresponds to degenerate triangles, where bodies are in syzygies â€“ instants when all three bodies are collinear. The geometry of the shape sphere can be understood by the Jacobi coordinates. Basically, Jacobi coordinates are relative coordinate vectors introduced by Carl Jacobi (even though Lagrange have reached them earlier), and they are represented by Equation (13).

The perk of these coordinate vectors is that they are mass-weighted vectors, so they work efficiently for the case of unequal masses in the three-body problem. To use the shape sphere for representing the solutions to the three-body problem, Equations (2) and (3) have to be solved then using and , so we can graphically represent the solutions.

Although the return proximity function and the shape sphere are not used directly in this paper, itâ€™s crucial to define them to understand the topologies of the solutions tested Montgomery, R. (2015).

**The ODE Solver**

As mentioned earlier in the introduction, the six differential equations associated with describing the equations of motion cannot be solved unless at specific initial conditions. For that, the code used in the simulation will use an ordinary differential equations solver, namely solveivp^{12}. The solver will be solving the ODEs numerically and is a part of the renowned python library SciPy. The method used in this paper to solve the ODEs will be the Rung-Kutta-Fehlberg Method for its precision in solving the equations while maintaining a non-draining memory usage^{13}.

**2. ****Selection of Families**

Through decades of work, hundreds of solutions to the three-body problem have been discovered; however, many of them showed instability after a short period of time, hence they are improbable to exist in real life. Contrary to that, some families with zero angular momentum in the system showed considerable stability over long periods of time, which implies that these systems are more likely to exist in real life. Consequently, this paper will focus mainly on testing the stability of zero-angular momentum families of solutions and categorizing them, all taken from the Three-body Gallery created by Suvakov and DmitraÅ¡inovi?. The tested solutions are the following: Figure-8; Butterfly-I; Butterfly-II; Bumblebee; Dragonfly; Goggles; Moth-I; and Moth-II^{14}.

**3. ****The Code**

In this section, significant snippets of the code will be explained.

**Solution**

This snippet represents laying out the time of simulation and calling out the function. NumPy library was used to lay down the line space of the time frame used. Here, N represents the number of points and T represents the time span of each step. For instance, **T = 0.001** means that each step is 1 millisecond. The function **ThreeBody** was solved using **solveivp**, which is a feature found in SciPy library. **solveivp** solves an initial value problem for a system of ordinary differential equations. This function numerically integrates a system of ordinary differential equations given an initial value:

where ( t ) is a 1-dimensional independent variable (time), and y(t) and f(t, y) are ( n )-dimensional vector-valued functions. These elements determine the differential equations. The goal is to find some y(t) , which approximately satisfies the given initial conditions (i.e., y(t_{0}) = y_{0} ).

**Stability Simulation**

In this piece of code, a main function was defined: **VectorStability**. The function takes one argument ( , , , , ) to determine what factor will be tested. After that the function iterates over a certain range to change either mass of the fourth body, initial positions, or the initial velocity.

# Results and Discussion

While testing the stability of the solutions, it was found that the system under the influence of a fourth body doesnâ€™t show linear stability; instead, islands of stability were found to be stable among instability regions, which caused some issues. To solve this problem, only a specific interval was taken into account. The range chosen for the mass was **(0.001, 0.01) **with an interval of **0.0001**. This can be justified by the fact that the introduced body will always have a very small mass in comparison with the masses of the existing three-body system; this can be found in the case of planets orbiting trinary star systems. In addition, initial positions and velocities varied between **0.1** and **1** with an interval of **0.01**.

After testing the systems for change in initial positions and velocities, it was found that the systems turn chaotic without forming any topologically clear configuration in addition to a loose gravitational binding. An example of one of the systems is shown in **fig. 2**, and the changes in x and y coordinates for the system are shown in **fig. 3**and **fig. 4**. For that, only the stability due to change in masses will be discussed in detail.

**Figure-8 Solution**

The figure-8 solution, shown in **fig. 5**, was found to be sensitive to any changes in the system; however, the addition of a fourth mass did not disrupt the gravitational binding of the system. The figure-8 solution with the addition of a fourth mass (0.0079) is shown in **fig. 6**. To have a better understanding of the disturbed system, the x and y coordinates of the four bodies were graphed, and theyâ€™re shown in **fig. 7** and **fig. 8 **respectively. Interestingly, the solution exhibited stability points at which the overall topology of the system was not changed significantly, and these points are resembled in the following list: [0.00163, 0.00170, 0.00260, 0.00830]. The solution at 0.00163 is shown in **fig. 9**.

**Butterfly-I**

The butterfly-I solution, which is shown in fig. 10, was found to be sensitive to any changes in the system. The addition of any fourth body disturbed the system significantly. The disturbed system after the addition of a fourth body with mass 0.00163 is shown in fig. 11. The x and y coordinates of the unstable system are shown in fig. 12 and fig. 13. We can see the drastic change in the y-coordinates that caused the chaos in the system.

**Butterfly-II**

The butterfly-II is shown in **fig. 14**. It was found to have a strong gravitational binding over

the range [0.001:0.0019]. Over the specified range, the system failed to exist. The only effect that took place was a compression in the coordinates in the system. This can be traced back to the increase in the gravitational force in the system due to the introduction of a fourth body. The system with the addition of a fourth body with a mass of 0.0019 is shown in fig. 15. To have the full picture of the system, the x and y coordinates of the system were graphed and are shown in fig. 16 and fig. 17.

**Bumblebee**

This solution showed a very strange behavior under the influence of the fourth body. The system, in the stable cases, showed a considerable deformation from the original configuration shown in fig. 18. However, among the gravitationally bound masses, certain points showed a gravitationally loose configuration with zero stability. The stability was in the range [0.001:0.0031]. fig. 19 is an example of gravitationally stable masses with its x and y coordinates shown in fig. 20 and fig. 21. The chaotic configuration is shown in fig. 22 and its x and y

coordinates are shown in fig. 23 and fig. 24. The points at which the system exhibited chaos are [0.0021, 0.0023, 0.0025, 0.0028, 0.003]. This solution can be marked as stable over the defined range except for the unstable points.

**Dragonfly**

The dragonfly solution, in **fig. 25**, showed a strong gravitational binding over the range from 0.001 to 0.0021, after which the system failed to exist. The solution was, however, totally distorted by the introduction of a new mass. The original system was turned into various distorted configurations; one of them is **fig. 26**. The x and y coordinates of the distorted configuration are shown in **fig. 27** and **fig. 28**. This solution can be labeled as stable over the mentioned range.

**Goggles**

The Goggles solution was very sensitive to the introduction of a fourth body to the system, showing instant instability over the whole range of masses. The original system, shown in **fig. 29**, was changed to the system in **fig. 30**. To further investigate the behavior of the system, the x and y coordinates were graphed in **fig. 31** and **fig. 32**. In the end, the system formed a two-body system and two single-body systems. This solution can be marked as unstable over the designated range.

**Moth-I**

The Moth-I solution showed a very strange behavior. Interestingly, the system was able to sustain its overall shape, with bodies moving faster, leading to the condensed graph; see **fig. 33** and **fig. 34** for comparison between the original configuration and the sustained configuration. Among the stability range, there existed points at which the system was unbound and split in a two-body system and two single-body systems. Those points are [0.0035, 0.0036]. The configuration after the addition of a fourth body with mass 0.0036 is shown in **fig. 35**. The x and y coordinates of the two four-body systems are shown in **fig. 35** and **fig. 36**. Finally, the stability range for the solution was [0.001:0.0041], after which the system became unstable and kept splitting. This solution can be labeled as stable over the determined range with the exception of the two unstable points.

**Moth-II**

This solution was very sensitive to the introduction of any additional body to the system. The system failed to exist after the mass 0.0014. In this small range of four masses, the solution was able to keep itself gravitationally bound, with great deviation from the original configuration, for two masses (0.0011 and 0.0012) and split for the other two. The original configuration is shown in **fig. 38**, while the gravitationally bound and the distorted solution are shown in **fig. 39** and **fig.** **40** respectively.

# Summary

The following table summarizes the results and the categorization of the solutions

# Conclusion

In conclusion, it was found that all three-body solutions are significantly affected by the introduction of a fourth body with any additional mass, even if it is 0.1% of the other bodies. The majority of the solutions experienced a drastic change in the configuration of the original system while maintaining a gravitational binding among the bodies. However, this was only on a certain range of masses. In addition, some solutions had only short ranges (or certain masses) where the solution either showed less chaotic configurations (like the Figure-8 Solution) or failed at all (Bumblebee Solution) among wide ranges of stability or instability. In addition, changing the position of the fourth body or its initial velocity led to complete chaos in the system. By further studying the stability of these chaotic systems with more computational power, we can grasp a better understanding of these complex systems and further work on the expanded n-body problem and incorporate the effects on the gravitational waves.

## Appendix

# References

- Laskar, J. (2012, September 26). Is the Solar System Stable?? arXiv.org. https://arxiv.org/abs/1209.5996 [↩] [↩]
- Poincare, H. (2015). Les mÃ©thodes nouvelles de la mÃ©canique cÃ©leste. Tome 1. Lilliad – UniversitÃ© De Lille – Sciences Et Technologies. https://iris.univ-lille.fr/handle/1908/3851 [↩]
- Bistafa, S. R. (2021). Eulerâ€™s three-body problem. Euleriana, 1(2), 181â€“187. https://doi.org/10.56031/2693-9908.1017 [↩]
- Broucke, R. (1975). On relative periodic solutions of the planar general three-body problem. Celestial Mechanics, 12(4), 439â€“462. https://doi.org/10.1007/bf01595390 [↩]
- Moore, C. (1993). Braids in classical dynamics. Physical Review Letters, 70(24), 3675â€“3679. https://doi.org/10.1103/physrevlett.70.3675; Chenciner, A. (2000, November 1). A remarkable periodic solution of the three-body problem in the case of equal masses. arXiv.org. https://arxiv.org/abs/math/0011268 [↩]
- Suvakov, M., & DmitraÅ¡inovi?, V. (2013). Three classes of Newtonian Three-Body Planar Periodic Orbits. Physical Review Letters, 110(11). https://doi.org/10.1103/physrevlett.110.114301; Li, X., Jing, Y. P., & Liao, S. (2018). Over a thousand new periodic orbits of a planar three-body system with unequal masses. Publications of the Astronomical Society of Japan, 70(4). https://doi.org/10.1093/pasj/psy057 [↩]
- Gutzwiller, M. C., & JosÃ©, J. V. (1990). Chaos in classical and quantum mechanics. In Springer eBooks. https://doi.org/10.1007/978-1-4612-0983-6 [↩]
- Montgomery, R. (2007). The zero angular momentum, three-body problem: All but one solution has syzygies. Ergodic Theory and Dynamical Systems, 27(6), 1933â€“1946. https://doi.org/10.1017/s0143385707000338 [↩] [↩]
- Suvakov, M., & DmitraÅ¡inovi?, V. (2014). A guide to hunting periodic three-body orbits. American Journal of Physics, 82(6), 609â€“619. https://doi.org/10.1119/1.4867608 [↩]
- Suvakov, M., & DmitraÅ¡inovic, V. (2014). A guide to hunting periodic three-body orbits. American Journal of Physics, 82(6), 609â€“619. https://doi.org/10.1119/1.4867608 [↩] [↩]
- Jacobi Coordinates. R. Montgomery (2015), The American Mathematical Monthly, Vol. 122, No. 4, 306 (http://www.jstor.org/stable/10.4169/amer.math.monthly.122.04.299) [↩]
- SciPy community (2024). scipy.integrate.solve_ivp â€” SciPy v1.11.4 Manual. https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp [↩]
- Poleshchikov, S. M. (2018). Regularization of the perturbed spatial restricted Three-Body problem by L-Transformations. Cosmic Research, 56(2), 151â€“163. https://doi.org/10.1134/s0010952518020077 [↩]
- Suvakov, M., & DmitraÅ¡inovi? (2013, March). Three-body Gallery. http://three-body.ipb.ac.rs/ [↩]