How to simulate stellar motion within a galaxy

You are working on problems that involve the gravitational forces between N = 2, or 3, or maybe 10, objects. As you are probably learning, the time required to perform all the calculations increases quickly with N.


   Q:  If a simulation of N = 10 objects takes 1 hour
       on a computer, how long will a similar simulation
       of N = 20 objects take?

   





   A:    4 hours

In fact, the number of direct calculations will rise as N².

Now, galaxies are big collections of stars. Really big. Really, really, really big.


   Q:  If a simulation of N = 10 objects takes 1 hour
       on a computer, how long will a similar simulation
       of N = 10^9 objects take?

   





   A:    10^16 hours  or approx 90,000 years

And yet, some scientists do claim to be able to simulate the motions of stars in galaxies.

How do they do it?

One approach: take a coarse approximation

One way to handle the gravitational interactions of very large numbers of objects is to focus on the large-scale forces and ignore -- or pay less attention -- to the small-scale forces.

Here's the zero-th order version of this method:

The mass of the Milky Way is about 10¹¹ solar masses, and it is arranged in a roughly spherical shape, with the highest density near the center. So, let's pretend that all the mass M is concentrated at a point in the very center.

Well, that's easy: the force on an object of mass m at a radial distance r is simply

You can compute the motions of tens of thousands of objects with a force law which is this simple! Of course, their motions don't look very much like those of real stars in a real galaxy.

So, let's look at the first-order approach, which is more sophisticated.

The mass of the Milky Way is about 10¹¹ solar masses, and it is arranged in a roughly spherical shape; however, the mass is not very concentrated to the center, so there is plenty of mass far from the center. Let's pretend that the density is spherically symmetric, and follows some simple relationship with radius.

So, for example, we could pretend that the density is

In order to compute the gravitational force on a star at some radial distance r from the center, we

integrate the density outward to r to find the mass M_int interior to r
modify Newton's Law of Universal Gravitation to include only this interior mass:

That would do a better job of reproducing the motions of stars in real galaxies -- depending on the exact values in the density relationship, we could get circular velocities which rise with radius out to some point.

What about non-spherical distributions of mass?

But one of the distinguishing characteristics of many galaxies, our own Milky Way among them, is their relative flatness:

Image courtesy of Bruce Hugo and Leslie Gaul and Adam Block

In order to deal with this real-life complication, we need to devise a more sophisticated approach. Consider the following multi-step technique:

1. make a model of the mass distribution: It helps to build the model out of very simple pieces; for the galaxy above, we might start with a sphere plus a flat circular disk
2. use the mass to compute the gravitational potential: Even when the mass distribution isn't spherical, we can use the Poisson equation
3. take the derivative of the potential to find the force at any point: Once we have a way to compute the gravitational potential at any location, we can easily compute a derivative numerically to find the force in any desired direction
4. move each particle according to the local force: You've already done this, right?
5. go to step 1: Since the particles have moved, we (may) need to re-compute the distribution of mass

This sort of technique will yield forces which are accurate on the large scales -- they will reflect well the forces on the Sun due to the conglomeration of billions of stars in the nucleus of a spiral galaxy -- but not on the small scales -- they will fail to account for the force of a rogue star which happens to fly through the solar system. Why? Because we generally use a very coarse grid to sample the density of matter.

For example, suppose we build a grid which has 1000 x 1000 x 1000 points in the X and Y and Z directions, and arrange it so that the entire Milky Way galaxy just fits inside. The diameter of the Milky Way is very roughly 100 kiloparsecs, so the width of each grid point will be

The distance to the nearest star system, Alpha Centauri, is about 1 parsec. Within a distance of 100 parsecs, there are roughly 200,000 stars (based on the incomplete count of 3803 stars in in this catalog of stars within 25 parsecs of the Sun. All these hundreds of thousands of real stars are smushed together into a single lump and placed at the center of the grid point.

Doesn't sound very accurate, does it? Well, sometimes, it's the best one can do.

For more information

Amara Graps' description of N-body integration techniques
Sverre Aarseth's paper "From NBODY1 to NBODY6: The Growth of an Industry" tells of the evolution of N-body codes.