N-Body Simulation of a Proto-Solar System with a Neptune-like Object


The widely-accepted explanation for the initial formation of the Solar System is known as the Nebular Hypothesis, that the Sun formed from a collapsing, dense region of an interstellar molecular cloud, with the remainder of material gravitationally accreting and flattening into an orbiting circular disk. This excess material eventually coalesced to form the small solar system bodies, including planets, moons, asteroids, and comets. After their initial formation, however, the orbital parameters of small bodies were not static. Due to complex gravitational interactions and interchanges of orbital energy between particles, planets, remaining small-bodies, and particles experienced 'migration' – periods during which their orbital parameters evolved substantially, particularly in their orbital radii (Semi-Major Axis). Migration is believed to account for the current configuration of planets in the solar system, including the gas giants. Additionally, it offers a likely explanation for the small orbital radii and, consequently, short orbital periods of “Hot Jupiter” planets which have been detected in extrasolar planetary systems, as it suggests that their orbits may have undergone an inward migration and eventual stabilization. Evidence suggests that migration could account for the current orbit of Neptune, since the orbit periods of a number of small solar system bodies appear to be in resonance, which must be the result of gravitational interactions with Neptune. A numerical simulation of a gas-less, solar-like proto-planetary disk was performed by Rodney S Gomes et. al. (Planetary Migration in a Planetesimal Disk: Why Did Neptune Stop at 30 AU? - 2004) to investigate Neptune's migration. The result was an outward migration of Neptune, which stabilized at a semi-major axis of approximately 30 AU from the central star, remarkably close to it's actual semi-major axis of 29.21 AU. Migration was found to be sensitive to the mass density of the proto-planetary disk. Outlying material, representing an initially close-in 'Kuiper Belt', quickly dispersed and lost a significant portion of its mass before Neptune stabilized into its final orbit, thus placing Neptune, effectively, on the edge of the disk after it ceased migrating. This result was obtained by setting the initial radius of the disk to ~35 AU, with a linear mass density of ~1.5 Earth masses per AU. When the mass density of the disk was increased (relatively massive disk), Neptune showed a runaway migration out to very large distances, due to interaction with particles that continually “fed” its migration. In a previous simulation of 10,000 particles, with a disk which encompassed 60 Earth masses of material between 20 and 45 AU, and and r^-(1.5) surface density profile, Gomes (2003) observed an outward migration of Neptune and stabilization at 45 AU.


Aside from small-scale simuations to to ensure that the code functioned as expected, I ran a single 30,000 particle simulation, which produced the main results presented in this project.The simulation was coded in FORTRAN 90 and paralellized to run on 16 cores. Here, I describe the initial parameters of the simulation, the numerical integration scheme, and the subsequent data analysis.

Initial Masses and State Vectors

Initial Coordinates

The initial 'disk' of particles was created by randomly generating x and y coordinates between 0 and 1, with the constraint that every random point must lie no further from the origin than 1. Viewed top-down, the result was a roughly circular collection of 30,000 points with a uniform number density per unit area. The disk was given some 'thickness' by assigning each coordinate a random z-component between 0 and 0.1. The first two such points represented the coordinates of the Sun and Neptune, respectively. The initial position of the Sun was set to [0,0,0] for simplicity. The initial position of Neptune was set to [0,-0.5,0], at an initial distance of 0.5 from the Sun. This radius was chosen so as to immerse Neptune deep within the protoplanetary disk in which it would have formed. Interior material represented that from which planets and asteroids interior to the orbit of Neptune might form. Outlying material represented an initially close-in 'Kuiper Belt'. Neptune's initial radius also determined the initial Neptune orbit period (assuming a circular orbit), which served as the time unit for the simulation. This was calculated with the following equation.

\begin {equation} \2 \pi \sqrt{a^3/M} \end {equation}

Initial Masses

Mass was allocated between the Sun, the disk particles, and Neptune such that the total mass of the system summed to 1. Though it actually constitutes over 99% of the total mass of the solar system, the mass of the sun was ar btrarily set to a value of 0.9. The mass of Neptune was set to a value of $4.5\times 10^-5$\textsuperscript{,t} so that the ratio of the mass of Neptune to the mass of the Sun was a realistic value (0.005\%). The disk particles, with a combined mass of 0.099955, constituted the remainder of the total mass of the system, which was equally distributed among them (mass of an individual disk particle = 0.099955/30,000 = $3.33 10^-6$.

Initial Velocities

Though the combined mass of disk particles and that of Neptune constituted a significant portion of the total mass of the system, the Sun was the dominant mass. Therefore, since it ies at the center of the system, the mass enclosed by the radius of each particle was assumed to be that of the Sun alone in calculating initial circular velocities. Due to the large number of particles and they way in which they were initially randomly distributed, the disk was assumed to have a uniform mass density in the xy-plane. Under these assumptions, every particle was assigned an initial velocity vector which would kick them, more or less, into a circular orbit about the Sun. This was accomplished with the following equation:

\begin {equation} V_c=\sqrt{\frac{GM}{R}} \end {equation}

Integration Scheme

A leap frog integration scheme was used to ensure conservation of total energy. Therefore, the trajectory of particles was calculated in discrete time steps. The time step duration was set to $dt = 0.02$, a number found to be small enough to adequately preserve total energy but large enough to run the code with a practical computation time (~ 2 days). Every time step, new accelerations, velocities, and positions were calculated for each particle, in that order, taking the previous state vectors of each particle and the gravitational influence of every other particle into account. Acceleration was calculated for every particle in a separate routine, and called upon in the main Do Loop. To set the velocity and position a half time step apart, as is required for the leap frog integration scheme, the initial velocity and position were redefined as the following:

\begin {equation} vi=vi+0.5 a dt \end {equation} \begin {equation} xi=xi+v_i dt \end {equation}

And in the main Do Loop, the integration proceeded in the following way:

$call$ $acceleration$ $subroutine$ \begin {equation} v=v+a dt \end {equation} \begin {equation} x=x+v dt \end {equation}

A force softening of $\epsilon = 0.001$ was included in the acceleration subroutine to prevent particles from being ejected from the system. This was adjusted until particle motion appeared smooth and total energy was preserved.

Every 1000 time steps, position data was output for every particle in the system, providing frequent 'snapshots' of the system over the course of the simulation.

Every 10 time steps, the program output Neptune's distance from the Sun as a function of time.

Every 50 time steps, the program output the kinetic, potential, and total energy of the system, as a function of time, calculated via a separate subroutine and called upon in the main Do Loop.

Data Analysis

To study the evolution of the proto-solar system model, I analyzed how the orbital parameters of the particles and Neptune changed over the course of the simulation. I particularly focused on the orbit parameters of Neptune, which could be studied quantitatively given the data output by the FORTRAN code. The FORTRAN code output calculated raw data, which was plotted and further analyzed using two IDL codes.

The first IDL code, orbit_parameters.pro (see attached code), took Neptune orbital radius and energy data as its input. This program identified the apoapsis and periapsis of every orbit by finding local minima/maxima in the radius data, and used these to calculate the semi-major axis $a$, orbital period $p$, and eccentricity $e$ as a function of time, via the following formulae:

\begin {equation} a=\frac{r{apo}+r{peri}} {2} \end {equation}

\begin {equation} p=2(t{apo}-t{peri}) \end {equation}

\begin {equation} e=\frac{r{apo}-r{peri}} {r{apo}+r{peri}} \end {equation}

This program also plotted the kinetic, potential, and orbital energy of the system as a function of time, using data directly from the FORTRAN code.

The second IDL code, plot_position.pro (see attached code), took position data for every particle, which was output by the Fortran code every 1000 time steps, and plotted these to create 'snapshots' of the actual system. This provides both an overhead view (xy-plane) and an edge-on view (xz-plane) of the system, at roughly 10 orbit intervals over the course of the 1000 orbit simulation (about 100 frames). Migration of disk particles other than Neptune, though not analyzed quantitatively, can be seen in the evolution of these 'snapshots' overtime, as well as a change in the overall shape of the disk.

Analytical Estimates

Migration occurs because of an interchange of angular momentum between disk particles and planets, such that the total angular momentum of the system remains constant. Depending on the relative angular momenta of these components, a planet can either undergo an outward or inward migration. Ida et al. (2000b) and Gomes (2004) offer an analytical explanation of this process, based on Rutherford Scattering. That is, during an encounter between two particles, the velocity vectors of the two particles instantaneously change. If the planet is placed on a circular orbit, then particles with a z-component of angular momentum larger than that of the planet will drive an outward migration. This is because, when crossing near the planet, their tangential velocities are larger than the circular velocity of the planet, thus accelerating it into a higher orbit. Particles with a z-component of angular momentum smaller than that of the planet, however, will drive an inward migration of the planet. This is because, when crossing near the planet, their tangential velocities are smaller than the circular velocity of the planet, thus negatively accelerating it into a lower orbit. Though the disk has some ’thickness’ due to the fact that particles were initially given non-zero z-components in their positions, and therefore begin on trajectories that are inclined, we will assume a simple model in which all particles follow perfectly circular orbits in the plane of the system (\(e,i = 0\)) In this case, all angular momentum is in the z-component, with a magnitude of \(a\times m\times v\) for each particle. In this particular simulation, the mass of Neptune was 4.5e-5, and the disk particle mass was approximately \(3.33\times 10^{-6}\) (about 7.4% the mass of Neptune). All bodies were initially placed on circular orbits ranging between a =0 and a = 1, with Neptune at a radius of a=0.5. Therefore, particles which were initially in proximity to Neptune had a angular momentum of approximately \(L=0.5\times 3.33\times 10^{-6}*v_n\) whereas Neptune had an angular momentum of approximately \(L=0.5\times 4.5\times 10^{-5}*v_n\) These differ only by a factor of the difference in their masses, so all particles which were initially near Neptune had an angular momentum smaller than that of Neptune. Therefore, we would expect to have observed an initial i