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, much-elongated 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 N-body 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.5-2 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²) 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/galaxy-formation

Space which evolution we simulate is defined as an array of Points:

1

class Space(val points: Array[Point])

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

case class Point(
  position: EuclideanVector,
  velocity: EuclideanVector,
  mass: Double = 1.0 // [M_sun]
) {
  def distanceVector(anotherPoint: Point): EuclideanVector = {
    EuclideanVector(
      anotherPoint.position.x - this.position.x,
      anotherPoint.position.y - this.position.y
    )
  }
}

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 Scala-way 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

package com.ontheserverside.galaxy

import java.time.LocalDateTime

import com.google.common.util.concurrent.AtomicDoubleArray
import com.ontheserverside.galaxy.Constants.{G, kmPcRatio}

import scala.annotation.tailrec

class NBody(
  totalSteps: Int,
  stepTimespan: Double,
  onStepCompleted: (Space, Int) => Unit,
  softeningLength: Double
) {
  def execute(space: Space): Space = execute(space, 1)

  @tailrec
  private[this] def execute(space: Space, stepsSoFar: Int): Space = {
    println(s"[${LocalDateTime.now}] Executing step $stepsSoFar out of $totalSteps")

    val transformed = executeStep(space)
    onStepCompleted(transformed, stepsSoFar)

    if (stepsSoFar == totalSteps) {
      transformed
    } else {
      execute(transformed, stepsSoFar + 1)
    }
  }

  private[this] def executeStep(space: Space): Space = {
    val forceFactorsX = new AtomicDoubleArray(space.points.length)
    val forceFactorsY = new AtomicDoubleArray(space.points.length)

    for (i <- 0 until space.points.length) {
      for (j <- (i + 1 until space.points.length).par) {

        val distanceVector = space.points(i).distanceVector(space.points(j))
        val softenedDistance = Math.pow((Math.pow(distanceVector.magnitude, 2) + Math.pow(softeningLength, 2)), 3/2)

        // [F] = [pc/M_sun * (km/s)^2 * M_sun^2 * 1/pc^3 * pc] = [ M_sun * km/s^2 * (km/pc) ]
        // in order to eliminate (km/pc) piece, we need to multiply by Constants.kmPcRatio
        val forceScalar = (G * space.points(i).mass * space.points(j).mass / softenedDistance) * kmPcRatio

        val forceFactor = (
          distanceVector.x * forceScalar,
          distanceVector.y * forceScalar
        )

        forceFactorsX.addAndGet(i, forceFactor._1)
        forceFactorsY.addAndGet(i, forceFactor._2)

        forceFactorsX.addAndGet(j, -1 * forceFactor._1)
        forceFactorsY.addAndGet(j, -1 * forceFactor._2)
      }
    }

    new Space(space.points.view.zipWithIndex.map { case (point, pointIdx) =>
      val newPosition = point.position.copy(
        x = point.position.x + stepTimespan * Constants.kmToPc(point.velocity.x),
        y = point.position.y + stepTimespan * Constants.kmToPc(point.velocity.y)
      )

      val newVelocity = point.velocity.copy(
        x = point.velocity.x + stepTimespan * forceFactorsX.get(pointIdx) / point.mass,
        y = point.velocity.y + stepTimespan * forceFactorsY.get(pointIdx) / point.mass
      )

      point.copy(position = newPosition, velocity = newVelocity)
    }.toArray)
  }
}

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³⁰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

object Constants {
  /**
    *  Gravitational constant expressed in parsecs [pc] and Solar mass [M_sun]
    */
  val G = 4.3e-3 // pc/M_sun (km/s)^2

  /**
    * Transforms units: [km] to [pc]
    *
    * 1 [pc] = 3.086e+13 [km]
    */
  def kmToPc(km: Double): Double = km / 3.086e+13

  /**
    * 1 [km] = 3.24078e-14 [pc]
    */
  val kmPcRatio = 3.24078e-14

  /**
    * Transforms years to seconds
    *
    * 1 [year] = 3.154e+7 [sec]
    */
  def yearsToSeconds(years: Double): Double = years * 3.154e+7
}

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:

1. Random generation of polar coordinates for the entire initial space area,
2. Uniform distribution of points for the entire initial space area,
3. 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

new Space(Array.fill(pointsCount) {
  val random = ThreadLocalRandom.current()

  val ρ = random.nextDouble(0, 10)
  val φ = random.nextDouble(0, Math.PI * 2)

  Point.fromPolarCoordinates(ρ, φ, velocity)
})

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

trait UniformPointDistributionSupport {
  /**
    * Lets to generate uniform distribution of random points in a circular ring.
    *
    * @param rMin - radius of the smaller circle determining area of annulus
    * @param rMax - radius of the bigger circle determining area of annulus
    * @param velocityFn - function used to evaluate point's velocity; should accept a distance from pole
    *                   to `Point` as an argument.
    */
  def createPoint(rMin: Double, rMax: Double, velocityFn: Double => Double): Point = {
    val random = ThreadLocalRandom.current()

    val ρ = random.nextDouble(0, 1)
    val φ = random.nextDouble(0, Math.PI * 2)

    val R = Math.sqrt((Math.pow(rMax, 2) - Math.pow(rMin, 2)) * ρ + Math.pow(rMin, 2))

    Point.fromPolarCoordinates(R, φ, rMax, velocityFn)
  }
}

Finally, we generate points from the user-provided 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

trait DensityFunctionGenSupport extends UniformPointDistributionSupport {
  /**
    * Generates `Space` with density described by the provided `densityFn`.
    *
    * The algorithm generates a set of circular rings with a width equal to `step` parameter.
    * Each ring is filled with `Point`s according to the provided `densityFn`.
    * Ring's density is computed from the mean density at its inner and outer boundaries and `Point`s are uniformly
    * distributed inside the ring.
    *
    * The algorithm generates rings until the `rMax` distance from the pole.
    *
    * Bear in mind that this is just an approximation and `step` value must be individually
    * adjusted for each `densityFn` and `rMax`.
    *
    * Each generated `Point` has a velocity computed from the provided `velocityFn` accepting a distance from
    * pole to the `Point` as an argument.
    */
  def generateSpaceFromDensityFn(
    densityFn: Double => Double,
    velocityFn: Double => Double,
    rMax: Double,
    step: Double = 0.01
  ): Space = {
    val generatedPoints = (for (r <- 0.0 until rMax by step) yield {
      val r2 = r + step

      val density = if (r == 0) densityFn(r2) else ((densityFn(r) + densityFn(r2)) / 2)
      val area = (Math.pow(r2, 2) - Math.pow(r, 2)) * Math.PI
      val numberOfPoints = (density * area).toInt

      Seq.fill(numberOfPoints)(createPoint(r, r2, velocityFn))
    }).flatten.toArray

    new Space(generatedPoints)
  }
}

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

object Point {
  // (...)
  def fromPolarCoordinates(ρ: Double, φ: Double, velocity: Double): Point = {
    fromPolarCoordinates(ρ, φ, 1.0, _ => velocity)
  }

  // (...)
  def fromPolarCoordinates(ρ: Double, φ: Double, R: Double, velocityFn: Double => Double): Point = {
    fromPolarCoordinates(ρ, φ, R, R, velocityFn)
  }

  /**
    * Creates an instance of `Point` from the provided polar coordinates.
    *
    * Requires to provide an additional εx and εy scale factors used to scale up / scale down cartesian coordinates
    * of the `Point`. This technique is useful when generating elliptical distributions or uniformly distributed random
    * points in a circular ring.
    *
    * Uses `velocityFn` to evaluate velocity of `Point` passing as an argument distance from pole to the `Point`.
    *
    * Created `Point` has a mass of 1 Solar masses.
    *
    * @param ρ - radial coordinate (distance from pole)
    * @param φ - angle in radians
    * @param εx - scale factor of X point's cartesian coordinate
    * @param εy - scale factor of Y point's cartesian coordinate
    * @param velocityFn - function used to evaluate point's velocity; should accept a distance from pole
    *                   to `Point` as an argument.
    */
  def fromPolarCoordinates(ρ: Double, φ: Double, εx: Double, εy: Double, velocityFn: Double => Double): Point = {
    val positionVector = EuclideanVector(
      ρ * εx * Math.cos(φ),
      ρ * εy * Math.sin(φ)
    )

    val velocity = velocityFn(positionVector.magnitude)

    val velocityVector = EuclideanVector(
      -1 * velocity * Math.sin(φ),
      velocity * Math.cos(φ)
    )

    new Point(positionVector, velocityVector)
  }
}

Here are sample galaxies formed from various profiles we defined.

System generated from CartwheelProfile compared to ESO 350-40 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 NGC-1300 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 N-body algorithm implementation. As we could expect, the direct sum algorithm is far from being super-effective. Although it’s good enough for simple simulations with profiles of 10⁴ or 10⁵ Solar Masses, it would be really hard to use it to simulate a Milky Way evolution which is 5.8 x 10¹¹ 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 master-worker 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 Barnes-Hut 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 real-world 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 350-40, NGC-1300 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 Time-Step in N-body Simulations. https://arxiv.org/pdf/astro-ph/0504573.pdf

[2] http://narimanfarsad.blogspot.com/2012/11/uniformly-distributed-points-inside_9.html

[3] L. Hernquist: An analytical model for spherical galaxies and bulges. http://adsabs.harvard.edu/abs/1990ApJ…356..359H