3 body star & planet orbits
Between 2018 and 2021, I simulated star orbits using Newtonian gravity, by numerically integrating the differential equation. The motivation for this was to be able to extend the orbits to three or more stars or planets, since the analytic solution doesn't exist for those cases. I have extended my star orbits simulation to compute the orbits of three stars, or two stars and a planet. Right now, the code has been almost entirely tested for the case of an inner binary and an outer planet. For the orbit shown below, the inner binary stars have masses of 100 in units where G=c=1. The outer planet has a mass of 1. I realize these aren't realistic masses for Newtonian gravity but I have made no relativistic approximation so the orbits should still produce reasonable results for this theory. The separation between the inner two stars is 10 and the separation between the center of the inner binary and the outer planet, initially, is 1000. The stars and planet begin in a colinear position along the x axis, with circular orbit initial conditions.
Newtonian Gravity Differential Equations
The differential equations solved by the simulation are given below. The left hand side represents the x, y, and z components of the acceleration and the right hand side represents the forces due to the gravitational attraction between each pair of particles in the x, y, and z directions. This is a full three body simulation, though we solve it in the limit where the mass ratio and distance ratio is large. The planet produces a force on the binary, which perturbs the binary's orbit, which in turn produces a force on the planet. We expect some energy transfer from this backreaction.
This image shows the nearly circular orbit of the outer planet (black) around the inner binary's orbit (red), which is too small to resolve because the planet's orbit is 100 times larger. The cyan line is the orbit the planet would have had if it had truly been on a circular orbit at a radius of 1000.
The deviations are not due to numerical error. Instead, they are for two reasons. One is that I set up the initial conditions with some imbalance of mass from the origin, but used velocities as if the center of mass were at the origin. This causes a slightly elliptical nature to the orbit, which can be seen in the figure below.
The second cause of the difference between the simulated orbit (black) and the naïve orbit (cyan) is that in a 3 body orbit, Kepler's laws are not obeyed, and the orbits are not actually precisely circular or elliptical.
The above plot shows the orbit of the inner binary. The circle in black shows one orbit of the binary, and the red curve shows the path the orbit traces as it goes through many cycles. There is some drift upward because there appears to be some linear momentum (peculiar velocity) associated with the initial conditions, due to the initial momentum of the planet (the initial momentums of the inner binary are equal and opposite). The spiraling motion is because the center of mass is not at the origin, leading to precession of the orbit about the center of mass as the planet's position changes. This is over a little more than one orbit of the outer planet.
These plots show the convergence of energy (E) and angular momentum (L) with increasing time step (dt). According to the theory of Newtonian gravity in physics, energy and angular momentum should be conserved, but we can see that they aren't, precisely, in these plots. Since the curves are similar but offset in scale, and decrease toward zero as timestep decreases in size, we can conclude that finite timestep resolution is what's responsible for this non conservation of angular momentum and energy. As timestep increases, the error in L and E increase due to approximations made in the simulation. At a smaller timestep of dt=0.03125, the relative error is 10^-8 even after a time of 14000. The fact that the error decreases by about 1.2 magnitudes on a log scale for every factor of 2 in dt means that the error scales as dt^4. This is consistent with a fourth order Runga Kutta Numerical integration method.
These plots show the linear momentum as a function of timestep size. Linear momentum is limited by roundoff error (machine precision) as step size is reduced. That explains the random nature of the curve.
This plot shows the evolution of the outer planet's energy as it passes through one whole orbit. The oscillation demonstrates that energy is transferred between the inner binary and outer planet. By showing that the timescale doesn't produce a ladder effect like the previous plot, and that there is no trend toward zero as the timescale gets smaller, we demonstrate that the energy oscillation is real and not a numerical artifact of the simulation.