Simulation of 3 Body Exo-Planetary System Orbits

To find out the general motion model of exo-planetary systems with one star and two planets, a computer program was used to carry out simulations and generate graphs showing the orbits of planets. When given the orbital periods and masses of the planets and stars, it is possible to predict the location of the planets over time and plot the shape of the orbit by considering the gravitational interactions between planets and the star, assuming that the planetary orbits are co-planar. I used the program to reproduce the result of transit timing variations (TTVs) of Kepler-46 system, I then investigated on the magnitude of transit timing variations on a 3-body system with various masses and periods. I also ran simulations to investigate the pattern of orbits for different periods of planets in order to get a systematic conclusion.


Introduction
At the time this paper is being written, 3053 stellar systems have already been discovered with at least one planet, and 664 of them have multiple planets (exoplanet.eu). Calculation of the planet orbits is impossible to be done analytically, but must be done numerically, which is challenging without the aid of a computer.
Currently, most researches regarding exo-planetary orbits focus on analyzing data from observation using methods like transit timing variation and transit duration variations to predict the orbits of planets (Pavol Gajdoˇs et al., 2019;Eric Agol et al., 2017;Nesvorn'y et al., 2012;Saad-Olivera et al., 2018). However, in this paper I chose to directly simulate the gravitational interactions between planets and stars using a program, which was a more useful tool for carrying out research on orbits and TTVs.
In this study, I focused on three body systems with one star and two planets. It is predicted that in a system like this the orbits would be less stable if the masses of planets were greater or the distance between planets was closer.
The algorithm was first deduced to calculate the motion of the planets from Kepler's laws and Newton's laws based on the mass of star and planets, and the orbital period of planets, since these quantities can be gathered directly by TTV and radial velocity measurements. Then a program was wrritten to carry out the algorithm and plot the transit timing variation of planets directly. I ran the program using the data of Kepler-46 system (Saad-Olivera et al., 2018), then compared the result of transit timing variation with existing data (Saad-Olivera et al., 2018) to check the accuracy of the program.
After the accuracy of the program was proved, I investigated the trend of magnitude of transit timing variation by running simulations in different periods and masses of planets.
Then I ran more simulations to find the more detailed pattern of change in periods.
Finally, I ran simulations which the periods of planets were in a fixed ratio but the mass was changing to investigate the impact of mass on the transit timing variation.

Formulas
To simplify the calculation three assumptions have been made: 1. The planets are orbiting at the same inclination angle-they're at the same plane 2. The star is stationary at the center of the model, (0,0) apr.ccsenet.org Applied Physics Research Vol. 12, No. 2;2020 3. Both planets start at location (r,0) The assumptions are, generally, close to the real situation since the 3-body planet systems I am modelling after generally contains a low inclination angle where 3° (e.g. Kepler-46 (Saad-Olivera et al., 2018), KOI-82 (Nesvorn'y et al.,2012)), hence assuming the planets to be co-planar can improve the performance of simulation without losing much accuracy. The second assumption is made as the star has a greater mass compared to the planet and can be consider as still relative to the planet system. Finally, the third assumption of both planets being aligned with the star is possible to happen, as long as two stars have a different velocity, at some point of time for any system.
In order to calculate the velocity of the planets, the semi-major axis of the planet based on Kepler's third law should be found first:

√
( 1 ) So it is possible to calculate the velocity of the planets by: ( 2 ) Also, based on Newton's second law of motion: ( And the Newton's law of gravitation: ( 4 ) The average acceleration of the planet due to star is: Where the represents the relative x coordinates of the planets, as the y coordinates between the planets. And is the distance between the planets at the time. Combining them get the total acceleration of the planet: The formula is similar with the second planet, the only key difference is the force between planets is in the opposite direction.
Then the coordinates of the planets could be found easily:

Program
A program was wrritten to carry out the calculation and plots the diagram of orbits.
Then the program handles the inputs, which include the mass of star and planets, as well as the orbiting period of planets; these data can be gathered directly by methods like TTV and RV. There's also a value that states the accuracy of simulation. The initial inputs are in natural units, where the masses are in terms of Sun mass and Earth mass and orbiting periods are in Earth days. In the next section, the program converts all the units to standard units: kg and year.
Utilizing the periods of the planets and the masses of the star and planets, it calculates the initial distances between the star and planets using the simulations above and records the locations of the planets on a 2-dimensional grid for every time step.
The program keeps looping until the time limit is reached, then the plugin can connect the points one by one from the arrays and draw the orbit of 2 planets on an image file.

Simulation of Observed Transit Timing Variations in Kepler-46 System
To check the accuracy of the program, the program was used to simulate the orbits for planets Kepler-46b and Kepler-46c in Kepler-46 system. First, the information of planets and star from paper (Saad-Olivera et al., 2018) was gathered, and the periods and masses was used to run the program. Part of the result of the simulation is shown in Figure 1, which is similar to the result given by Nesvorn'y et al. in "The Detection and Characterization of a Non-transiting Planet by Transit Timing Variations" through observation ( Figure 2). There are some slight differences, which may be due to the assumptions that the orbits are co-planar and both planets start at the same orbit phase. This shows the result of simulation is generally reliable as both graphs have a pattern with sharp turning angles, identical magnitudes and similar frequency. The magnitude of transit timing variation is nearly the same as Figure 2.

Investigation of Orbital Period Variation
The simulations were run to find the trend on three body orbits.
To do the investigation systematically, the program was updated and made to run through simulations when the mass of the outer planet is 100-180 Earth mass, and the period of outer planet is 75-175 days, with mass of inner planet fixed at 1 Earth mass and period at 50 days. Then I recorded the magnitude of the transit timing variation (different between the longest and shortest period of inner planet) of 25 simulations and plot it in a graph:  Vol. 12, No. 2;2020 It is assumed that this is because when the orbits of planets get closer or when the mass of the second planet gets bigger, the gravitational interaction between them becomes greater, so the orbits fluctuate more.
Then, I investigated on the more detailed change of period against orbits for different ratio of periods between planets when the masses of planets are constant.
To test this I ran the simulation multiple times, and changed the period of the inner planet.
The mass of inner planet is set to be 1 Earth mass, the mass of outer planet to be 1 Jupiter mass, and mass of the star to be 1 solar mass. The period of the outer planet is fixed at 200 days.
The orbits of planets and the change in periods are shown in Figure 4.
First panel and second panel: In the situation where the ratio of periods between planets is less than 1:2, the orbits of both stars remain in a near perfect circle with the periods keep constant;

Second
Third panel: When the ratio of periods between planets is 1:2, the inner planet's orbit has an obvious periodic variation, but the outer planet's orbit is still near constant, due to its high mass; Fourth panel: As the orbits of planets becoming closer the variation of the orbit of inner planet has a higher frequency, while the outer planet still has a near constant orbit.
Fifth panel: When the ratio of periods reaches 2:3, the pattern of variation changes dramatically: where it no longer changes periodically but rather randomly. The magnitude of change also increases greatly from around 5e5 seconds to 5e7 seconds.
The only variable changed during the simulations is the period of the inner planet, which translates to that the distances of the planets are getting closer. Considering The Newton's Law of Gravitation: We can deduce that as the proximity of the planet increases, the forces acting on both planets increase. However since the outer planet's mass is 318 times the mass of inner planet, the force has less impact on its acceleration, and hence, its orbit. While the inner planet's orbit is disturbed significantly due to the additional gravitational force acting on it from outer planet.
In addition, as the distance from the star increases, the centripetal force it receives from the star decreases, so the orbit becomes more eccentric as the period increases.
It is notable on the last experiment, when the period of the inner planet is 100 days, which is exactly half that of the outer planet. The period of the planet changes periodically. Figure 6. When the graph is magnified we can see that the period of the inner planet fluctuates up and down rapidly, and generally form a sine curve To investigate the general trend of period when the periods of 2 planets are in 1:2 ratio when the ratio of mass is changing, I did more simulation where the periods of planets were fixed at 100 and 200 days, and mass of the inner planet at 1 earth mass, then changed the mass of outer planet.
Figure 7. This graph shows the pattern of periods over orbits in different outer planet mass, with other condition the same. This shows that the pattern of change is always a periodic variation, which as the mass of planet 2 gets greater, the frequency of change of the inner planet's period increases It is assumed that this is because as the mass of outer planet getting greater, the gravitational interaction on the inner planet is greater so the period changes more significantly.

Conclusion
A program was wrritten to simulate the orbits of planets in a three body exo-planetary system based on the gravitational interactions between planets and star.
The simulations have shown that when the orbits of 2 planets are closer when the masses of planets are greater, the variations of planets' orbits are greater, possibly because of the gravitational interaction of planets are more significant based on the newton's law of gravitation.
Also, when the masses of planets are constant, when the ratio of periods is smaller than 1:2 the orbits of inner planets are generally constant. When the ratio is 1:2-2:3, the period changes in a uniform wave, however, when the ratio is more than 2:3, the period fluctuates in a greater magnitude randomly. Simultaneously, the orbit of outer planet is always constant, it is assumed that because of the mass of outer planet is much greater.
The features of the period of the inner planet when the ratio of periods is 1:2 is also investigated, which the change of period shows an exact sine wave. The result of simulations shows that when the ratio of masses between planets becomes smaller, the frequency becomes greater.
The results from simulations provided a method to assess the stability of planet system models. Hence it is possible to predict the possibility for certain kind of planetary system to exist.
The reliability of transit timing variation of Kepler-46b was also proved by running simulations to calculate transit timing variation, and compared the result of the simulation with the observed data from Kepler Mission. The results obtained by simulation and observation are similar, and both indicated there were more planet(s) existing in Kepler-46 system.
More research can be done based on the current investigation. For example, it would be possible to expand the program and enable it to simulate orbits for planetary systems with more planets, so the range of application would become wider. Also, although most planetary systems that are discovered contains planets that are co-planar or near co-planar, the accuracy of the simulation would be improved if the program considered the