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 p