If you gaze out into the dark sky not polluted by artificial city lights, you have a chance to see a part of the Milky Way or a diffuse, muchelongated disk of the Andromeda galaxy. They are made of stars, stellar remnants, dust, interstellar gas and dark matter all bound together by the gravity. This blog post aims at implementing Nbody algorithm for simulating a galaxy assembly.
Before we go any further we need to talk about physics and recall Newton’s law of universal gravitation: Every point mass attracts every single other point mass by a force pointing along the line intersecting both points. The force is directly proportional to the product of the two masses and inversely proportional to the square of the distance between the point masses:
If we represent point’s position in space in time t
with a position vector, that is, with a vector having tail fixed in origin and head attached to a moving point, we can express the gravitational force value at the time t
as:
where r_{i} is a position vector for point i
, r_{j} is a position vector for point j
and r_{ij} = r_{j} – r_{i}.
Of course, computing a force between two points is not enough to simulate a galaxy which is made of a billion of visible and invisible objects. We need to compute the gravity force between each one of them.
This equation is what we need to solve for every point of mass. For simplicity, we limit our simulation to two dimensions (although it can be easily extended to 3 dimensions), we assume galaxy is made of N mass points and we ignore magnetohydrodynamics and all other effects except the gravity. Our drawing algorithm renders every single point of mass but in a reality, some of them should be treated as a baryonic matter and some of them as an invisible dark matter.
The equation looks very nice but it still does not tell us how our galaxy evolves. There are two things we need to know: how the velocity and position of a point change in time. For this, we’ll need to use Newton’s second law of motion:
Since acceleration is defined as the derivative of velocity, V
, with respect to time t
and velocity is defined as the derivative of position, X
, with respect to time, acceleration can be thought of as the second derivative of X
with respect to t
:
Now, we may define how velocity and position change over the course of time:
There is still one important thing missing here that is not easy to spot. The position of a point in our galaxy is represented by two doubles. Our simulation is collisionless so it means two points can come very close to each other (they are limited only by the range of a number representation) and the force will go to infinity. There is a numerical trick called softening allowing to prevent this kind of situations. The procedure modifies gravitational potential of each point in order to alter the law of gravity at small distances:
where ϵ is the softening length determining the degree of potential smoothing. Studies [1] show that the softening length should be a factor of 1.52 smaller than the mean distance between particles in the densest regions to be resolved.
Before we jump to the code we may still think of one small optimization to our algorithm. If F_{ab} represents force of gravity between points a
and b
, then:
It means we can slightly optimize our O(n^{2}) algorithm by computing scalar of the gravitational force between two given points only once.
The simulation is implemented in Scala and entire code is available on GitHub: https://github.com/pdyraga/galaxyformation
Space
which evolution we simulate is defined as an array of Point
s:
1


and each Point
is described by position
and velocity
vectors as well as a mass
:
1 2 3 4 5 6 7 8 9 10 11 12 

Here is the code evaluating system evolution. We use Scala’s parallel collections to split computations between concurrent threads and AtomicDoubleArray
to ensure safe writes. We tried several approaches here including more Scalaway solutions like folding, views or even new strawman collections. Still, the ordinary loop with AtomicDoubleArray
proved to be the most performance effective, sometimes even several times faster.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 

The simulation code contains some mysterious [km][pc]
unit conversions. The reason for doing them is that astronomical objects are large and operating on units we use on Earth to describe celestial objects is unwieldy and sometimes even not possible due to the range limits of a computer representation of numbers.
For this reason, we use parsec as a unit of length and Solar mass as a unit of mass. One parsec is equal to about 3.26 light years. One Solar mass is equal to the mass of Sun, about 2 * 10^{30}kg. The gravitational constant is scaled to parsecs and Solar masses as well as we define a number of conversions that are applied when evaluating gravitational force and a new position of an object:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 

Finally, we have to discuss how to generate initial conditions for the simulation. We defined a set of initial profiles in the com.ontheserverside.galaxy.profile
package which combine the following techniques:
 Random generation of polar coordinates for the entire initial space area,
 Uniform distribution of points for the entire initial space area,
 Generation of points from a density function for the initial space area.
The first technique is the easiest one and is used to generate profiles for carthweel and spiral galaxies. Two numbers representing point’s position in polar coordinates are randomly generated. Then, polar coordinates are transformed into cartesian coordinates. It results in a distribution with particle density increasing towards the center of the galaxy which simulates the greater mass in its center.
1 2 3 4 5 6 7 8 

The second technique uses the joint probability density function to generate uniformly distributed random points inside a circular ring [2]:
Generating uniform distribution inside a circle is just a special case of this technique where the inner ring’s circle radius is equal 0.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 

Finally, we generate points from the userprovided density function. In this case, we divide the space into circular rings inside which we use the joint probability density function to generate a uniform distribution of random points. The expected density of the ring is calculated as a mean density for the inner and outer ring boundary. This technique is used to generate a spherical galaxy model [3] in the BulgeProfile
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 

Code transforming polar coordinates into cartesian coordinates is placed in the Point
companion object. Each generated point is given an initial velocity to avoid a collapse of the entire system to the center. The initial velocity vector is always tangent to the circle having a center in the center of the initially generated system.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 

Here are sample galaxies formed from various profiles we defined.
System generated from CartwheelProfile
compared to ESO 35040 cartwheel galaxy
You can see the system evolution on YouTube. Please switch to HD fullscreen mode for the optimal quality.
System generated from SpiralProfile
compared to NGC1300 spiral galaxy.
You can see the system evolution on YouTube. Please switch to HD fullscreen mode for the optimal quality.
System generated from BulgeProfile
compared to Messier 54 globular cluster.
You can see the evolution of a small cluster generated from UniformProfile
on YouTube. Please switch to HD fullscreen mode for the optimal quality.
There are few lessons we learned during the Nbody algorithm implementation. As we could expect, the direct sum algorithm is far from being supereffective. Although it’s good enough for simple simulations with profiles of 10^{4} or 10^{5} Solar Masses, it would be really hard to use it to simulate a Milky Way evolution which is 5.8 x 10^{11} Solar Masses and it would require some serious effort even from a supercomputer. Speaking of the supercomputer, we found it impossible to convince JVM to take advantage of the whole computing power it offered. No matter what we tried, just a small percentage of available cores was utilized.
Having this experience, we now think about a masterworker architecture with Scala master server and OpenMP C/C++ workers as one of the best possibilities for further implementation. Depending on the performance of this approach, we may also switch from a direct sum algorithm into the BarnesHut approximation having O(n log n) complexity. Instead of comparing each point of mass with all other points, the system is divided into cells treated as a single large particle having a center at the cell’s center of mass and only particles from nearby cells are treated individually. Another possible improvement could be an implementation of more complex initial profiles better reflecting realworld conditions by, for example, taking into account dark matter halo.
This work has been made possible with the help of Tobiasz Górecki, PhD candidate in astronomy at Jagiellonian University.
ESO 35040, NGC1300 and Messier 54 photos are a courtesy ESA/Hubble & NASA.
References
[1] S.A. Rodionov, N.Ya. Sotnikova: Optimal Choice of the Softening Length and TimeStep in Nbody Simulations. https://arxiv.org/pdf/astroph/0504573.pdf
[2] http://narimanfarsad.blogspot.com/2012/11/uniformlydistributedpointsinside_9.html
[3] L. Hernquist: An analytical model for spherical galaxies and bulges. http://adsabs.harvard.edu/abs/1990ApJ…356..359H